1 //____________________________________________________________________
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.
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.
12 // Use the script `Compile.C' to compile this class using ACLic.
15 #include <AliFMDHit.h>
16 #include <AliFMDDigit.h>
17 #include <AliFMDRecPoint.h>
18 #include <AliFMDInput.h>
19 #include <AliFMDEdepMap.h>
20 #include <AliFMDFloatMap.h>
25 #include <TParticle.h>
26 #include <TClonesArray.h>
31 //____________________________________________________________________
32 /** @class DrawHitsRecs
33 @brief Draw hit energy loss versus rec point mult
36 Root> Compile("DrawHitsRecs.C")
42 class DrawHitsRecs : public AliFMDInput
46 TH2D* fHitEvsRecM; // Histogram
47 TH2D* fHitEvsRecE; // Histogram
48 TH1D* fDiffE; // Histogram
49 TH2D* fHitsVsRecM; // Histogram
56 //__________________________________________________________________
57 TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
60 Float_t dp = n / TMath::Log10(max / min);
61 Float_t pmin = TMath::Log10(min);
63 for (Int_t i = 1; i < n+1; i++) {
64 Float_t p = pmin + i / dp;
65 bins[i] = TMath::Power(10, p);
69 //__________________________________________________________________
70 DrawHitsRecs(Bool_t primary=kFALSE,
71 Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
72 Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
78 if (fPrimary) AddLoad(kKinematics);
79 TArrayF eloss(MakeLogScale(n, emin, emax));
82 for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
84 fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
85 n, emin, emax, 1025, -.5, 1024.5);
86 fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
87 fHitEvsAdc->SetYTitle("ADC");
89 fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
90 eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
91 fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
92 fHitEvsRecM->SetYTitle("M_{rec}");
94 fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
95 n, emin, emax, n, emin, emax);
96 fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
97 fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
100 fDiffE = new TH1D("diffE",
101 "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
103 fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
108 fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
109 o, omin, omax, m, mmin, mmax);
110 fHitsVsRecM->SetXTitle("# of Hits");
111 fHitsVsRecM->SetYTitle("M_{rec}");
113 //__________________________________________________________________
114 /** Begining of event
115 @param ev Event number
116 @return @c false on error */
117 Bool_t Begin(Int_t ev)
120 return AliFMDInput::Begin(ev);
122 //__________________________________________________________________
123 Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
125 // Cache the energy loss
126 if (!hit) return kFALSE;
127 UShort_t det = hit->Detector();
128 Char_t rng = hit->Ring();
129 UShort_t sec = hit->Sector();
130 UShort_t str = hit->Strip();
132 AliWarning(Form("Bad strip number %d in hit", str));
136 TParticle* kine = fStack->Particle(hit->Track());
137 if (!kine) return kTRUE;
138 if (!kine->IsPrimary()) return kTRUE;
141 fMap(det, rng, sec, str).fEdep += hit->Edep();
142 fMap(det, rng, sec, str).fN++;
146 //__________________________________________________________________
147 Bool_t ProcessDigit(AliFMDDigit* digit)
149 if (!digit) return kTRUE;
151 UShort_t det = digit->Detector();
152 Char_t rng = digit->Ring();
153 UShort_t sec = digit->Sector();
154 UShort_t str = digit->Strip();
156 AliWarning(Form("Bad strip number %d in digit", str));
159 Double_t edep = fMap(det, rng, sec, str).fEdep;
160 if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
165 //__________________________________________________________________
166 Bool_t ProcessRecPoint(AliFMDRecPoint* single)
168 if (!single) return kTRUE;
169 UShort_t det = single->Detector();
170 Char_t rng = single->Ring();
171 UShort_t sec = single->Sector();
172 UShort_t str = single->Strip();
174 AliWarning(Form("Bad strip number %d in single", str));
177 Double_t edep = fMap(det, rng, sec, str).fEdep;
178 Int_t nhit = fMap(det, rng, sec, str).fN;
180 fHitEvsRecM->Fill(edep, single->Particles());
181 fHitEvsRecE->Fill(edep, single->Edep());
182 fDiffE->Fill((single->Edep() - edep) / edep);
184 if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
187 //__________________________________________________________________
190 gStyle->SetPalette(1);
191 gStyle->SetOptTitle(0);
192 gStyle->SetCanvasColor(0);
193 gStyle->SetCanvasBorderSize(0);
194 gStyle->SetPadColor(0);
195 gStyle->SetPadBorderSize(0);
197 new TCanvas("c0", fHitEvsAdc->GetTitle());
198 fHitEvsAdc->SetStats(kFALSE);
199 fHitEvsAdc->Draw("COLZ");
201 new TCanvas("c1", fHitEvsRecM->GetTitle());
202 fHitEvsRecM->SetStats(kFALSE);
203 fHitEvsRecM->Draw("COLZ");
205 new TCanvas("c2", fHitEvsRecE->GetTitle());
206 fHitEvsRecE->SetStats(kFALSE);
207 fHitEvsRecE->Draw("COLZ");
209 new TCanvas("c3", fDiffE->GetTitle());
212 new TCanvas("c4", fHitsVsRecM->GetTitle());
213 fHitsVsRecM->SetStats(kFALSE);
214 fHitsVsRecM->Draw("COLZ");
219 ClassDef(DrawHitsRecs,0);
222 //____________________________________________________________________