FMD/scripts/DrawSDigits.C

   1 //____________________________________________________________________
   2 //
   3 // $Id: DrawSDigits.C 20907 2007-09-25 08:44:03Z cholm $
   4 //
   5 // Script that contains a class to draw eloss from hits, versus ADC
   6 // counts from digits, using the AliFMDInputHits class in the util library.
   7 //
   8 // It draws the energy loss versus the p/(mq^2).  It can be overlayed
   9 // with the Bethe-Bloc curve to show how the simulation behaves
  10 // relative to the expected.
  11 //
  12 // Use the script `Compile.C' to compile this class using ACLic.
  13 //
  14 #include <TH1D.h>
  15 #include <TH2D.h>
  16 #include <TProfile2D.h>
  17 #include <TCanvas.h>
  18 #include <TMath.h>
  19 #include <AliFMDSDigit.h>
  20 #include <AliFMDInput.h>
  21 #include <AliFMDEdepMap.h>
  22 #include <AliFMDGeometry.h>
  23 #include <iostream>
  24 #include <TStyle.h>
  25 #include <TLatex.h>
  26 #include <TArrayF.h>
  27 #include <AliLog.h>
  28
  29 /** @class DrawSDigits
  30     @brief Draw hit energy loss versus digit ADC
  31     @code
  32     Root> .L Compile.C
  33     Root> Compile("DrawSDigits.C")
  34     Root> DrawSDigits c
  35     Root> c.Run();
  36     @endcode
  37     @ingroup FMD_script
  38  */
  39 class DrawSDigits : public AliFMDInput
  40 {
  41 private:
  42   TH1D* fAdc; // Histogram
  43   TProfile2D* fPrimRatio[5]; // Histograms
  44 public:
  45   void Idx2Det(UShort_t idx, UShort_t& d, Char_t& r) const
  46   {
  47     d = 0;
  48     r = '\0';
  49     switch (idx) {
  50     case 0: d = 1; r = 'I'; break;
  51     case 1: d = 2; r = 'I'; break;
  52     case 2: d = 2; r = 'O'; break;
  53     case 3: d = 3; r = 'I'; break;
  54     case 4: d = 3; r = 'O'; break;
  55     }
  56   }
  57   Short_t Det2Idx(UShort_t d, Char_t r) const
  58   {
  59     Short_t idx = -1;
  60     switch (d) {
  61     case 1: idx = 0; break;
  62     case 2: idx = 1; break;
  63     case 3: idx = 3; break;
  64     }
  65     return (idx + ((r == 'O' || r == 'o') ? 1 : 0));
  66   }
  67
  68
  69   //__________________________________________________________________
  70   DrawSDigits(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1023.5)
  71     : AliFMDInput("galice.root")
  72   {
  73     AddLoad(kSDigits);
  74     AddLoad(kGeometry);
  75     fAdc = new TH1D("adc", "ADC", m, amin, amax);
  76     fAdc->SetXTitle("ADC value");
  77     fAdc->Sumw2();
  78
  79     Int_t   eN   = 50;
  80     Float_t eMin = -4;
  81     Float_t eMax = 6;
  82     Int_t   pN   = 40;
  83     Float_t pMin = -.5;
  84     Float_t pMax = 359.5;
  85     for (UShort_t i = 0; i < 5; i++) {
  86       UShort_t d;
  87       Char_t   r;
  88       Idx2Det(i, d, r);
  89       fPrimRatio[i] = new TProfile2D(Form("primRatio_fmd%d%c", d, r),
  90                                      Form("Primary/Total - FMD%d%c", d, r),
  91                                      eN, eMin, eMax, pN, pMin, pMax);
  92       fPrimRatio[i]->SetXTitle("#eta");
  93       fPrimRatio[i]->SetYTitle("#varphi");
  94       fPrimRatio[i]->SetZTitle("M_{ch,primary}/M_{ch,total}");
  95       fPrimRatio[i]->GetXaxis()->SetTitleFont(132);
  96       fPrimRatio[i]->GetYaxis()->SetTitleFont(132);
  97       fPrimRatio[i]->GetZaxis()->SetTitleFont(132);
  98       fPrimRatio[i]->GetXaxis()->SetLabelFont(132);
  99       fPrimRatio[i]->GetYaxis()->SetLabelFont(132);
 100       fPrimRatio[i]->GetZaxis()->SetLabelFont(132);
 101       fPrimRatio[i]->GetXaxis()->SetTitleSize(.06);
 102       fPrimRatio[i]->GetYaxis()->SetTitleSize(.06);
 103       fPrimRatio[i]->GetZaxis()->SetTitleSize(.06);
 104       fPrimRatio[i]->GetXaxis()->SetLabelSize(.06);
 105       fPrimRatio[i]->GetYaxis()->SetLabelSize(.06);
 106       fPrimRatio[i]->GetZaxis()->SetLabelSize(.06);
 107       fPrimRatio[i]->GetXaxis()->SetNdivisions(10);
 108       fPrimRatio[i]->GetYaxis()->SetNdivisions(10);
 109       fPrimRatio[i]->GetZaxis()->SetNdivisions(10);
 110     }
 111   }
 112   Bool_t Init()
 113   {
 114     Bool_t ret = AliFMDInput::Init();
 115     AliFMDGeometry::Instance()->Init();
 116     AliFMDGeometry::Instance()->InitTransformations();
 117     return ret;
 118   }
 119
 120   //__________________________________________________________________
 121   Bool_t ProcessSDigit(AliFMDSDigit* digit)
 122   {
 123     if (!digit) return kTRUE;
 124     fAdc->Fill(digit->Counts());
 125     digit->Print("lp");
 126     if (digit->NParticles() == 0) return kTRUE;
 127
 128
 129     Short_t primIdx = Det2Idx(digit->Detector(), digit->Ring());
 130     if (primIdx < 0) return kTRUE;
 131
 132     AliFMDGeometry* geom = AliFMDGeometry::Instance();
 133     Double_t x, y, z;
 134     geom->Detector2XYZ(digit->Detector(),
 135                        digit->Ring(),
 136                        digit->Sector(),
 137                        digit->Strip(),
 138                        x, y, z);
 139     // We should get the primary vertex and do
 140     // z     += fCurrentVertex;
 141     Double_t phi   =  TMath::ATan2(y, x) * 180. / TMath::Pi();
 142     Double_t r     =  TMath::Sqrt(y * y + x * x);
 143     Double_t theta =  TMath::ATan2(r, z);
 144     Double_t eta   = -TMath::Log(TMath::Tan(theta / 2));
 145     Double_t ratio = digit->NPrimaries() / digit->NParticles();
 146     if (phi < 0) phi += 360;
 147     fPrimRatio[primIdx]->Fill(eta, phi, ratio);
 148     return kTRUE;
 149   }
 150   //__________________________________________________________________
 151   Bool_t Finish()
 152   {
 153     gStyle->SetPalette(1);
 154     gStyle->SetOptTitle(0);
 155     gStyle->SetCanvasColor(0);
 156     gStyle->SetCanvasBorderSize(0);
 157     gStyle->SetPadColor(0);
 158     gStyle->SetPadBorderSize(0);
 159
 160     TCanvas* c1 = new TCanvas("fmdSdigitSpectra", "FMD SDigit spectra");
 161     fAdc->SetStats(kFALSE);
 162     if (fAdc->GetEntries() != 0)
 163       fAdc->Scale(1. / fAdc->GetEntries());
 164     fAdc->Draw();
 165     c1->Modified();
 166     c1->Update();
 167     c1->cd();
 168
 169     TCanvas* c2 = new TCanvas("fmdSdigitPrims", "FMD SDigit # prim/# ntotal",
 170                               800, 800);
 171     c2->SetTopMargin(0);
 172     c2->SetRightMargin(0);
 173     c2->Divide(2, 3, 0, 0);
 174     for (UShort_t i = 0; i < 5; i++) {
 175       Int_t        np = i+1+(i == 0 ? 0 : 1);
 176       TVirtualPad* p   = c2->cd(np);
 177       p->SetFillColor(0);
 178       p->SetTopMargin((np <= 2) ? 0.05 : 0);
 179       p->SetLeftMargin((np % 2 == 1) ? 0.20 : 0);
 180       p->SetRightMargin((np % 2 == 1) ? 0 : 0.20);
 181       p->SetBottomMargin((np >= 5) ? 0.15 : 0);
 182       fPrimRatio[i]->SetStats(kFALSE);
 183       fPrimRatio[i]->Draw(Form("COL%c", (np % 2 == 1) ? ' ' : 'Z'));
 184       UShort_t d;
 185       Char_t   r;
 186       Idx2Det(i, d, r);
 187       TLatex* text = new TLatex(0, 180, Form("FMD%d%c", d, r));
 188       text->SetTextFont(132);
 189       text->SetTextAlign(22);
 190       text->SetTextSize(0.08);
 191       text->Draw();
 192     }
 193     c2->Modified();
 194     c2->Update();
 195     c2->cd();
 196
 197     return kTRUE;
 198   }
 199
 200   ClassDef(DrawSDigits,0);
 201 };
 202
 203 //____________________________________________________________________
 204 //
 205 // EOF
 206 //