The patch fixes
[u/mrichter/AliRoot.git] / FMD / scripts / DrawHitsRecs.C
1 //____________________________________________________________________
2 //
3 // $Id$
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 <TH2D.h>
15 #include <AliFMDHit.h>
16 #include <AliFMDDigit.h>
17 #include <AliFMDRecPoint.h>
18 #include <AliFMDInput.h>
19 #include <AliFMDEdepMap.h>
20 #include <AliFMDFloatMap.h>
21 #include <iostream>
22 #include <TStyle.h>
23 #include <TArrayF.h>
24 #include <TCanvas.h>
25 #include <TParticle.h>
26 #include <TClonesArray.h>
27 #include <TTree.h>
28 #include <AliStack.h>
29 #include <AliLog.h>
30
31 //____________________________________________________________________
32 /** @class DrawHitsRecs
33     @brief Draw hit energy loss versus rec point mult
34     @code 
35     Root> .L Compile.C
36     Root> Compile("DrawHitsRecs.C")
37     Root> DrawHitsRecs c
38     Root> c.Run();
39     @endcode
40     @ingroup FMD_script
41  */
42 class DrawHitsRecs : public AliFMDInput
43 {
44 private:
45   TH2D* fHitEvsAdc;
46   TH2D* fHitEvsRecM;  // Histogram 
47   TH2D* fHitEvsRecE;  // Histogram 
48   TH1D* fDiffE;       // Histogram 
49   TH2D* fHitsVsRecM;  // Histogram 
50   AliFMDEdepMap  fMap;
51   AliFMDFloatMap fEta;
52   AliFMDFloatMap fPhi;
53   AliFMDFloatMap fMult;
54   Bool_t fPrimary;
55 public:
56   //__________________________________________________________________
57   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
58   {
59     TArrayF bins(n+1);
60     Float_t dp   = n / TMath::Log10(max / min);
61     Float_t pmin = TMath::Log10(min);
62     bins[0]      = 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);
66     }
67     return bins;
68   }
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) 
73   { 
74     fPrimary = primary;
75     AddLoad(kRecPoints);
76     AddLoad(kHits);
77     AddLoad(kDigits);
78     if (fPrimary) AddLoad(kKinematics);
79     TArrayF eloss(MakeLogScale(n, emin, emax));
80     TArrayF mults(m+1);
81     mults[0] = mmin;
82     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
83
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");
88     
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}");
93
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]");
98
99
100     fDiffE = new TH1D("diffE", 
101                       "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}", 
102                       1100, -1, 1.1);
103     fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
104     
105     Double_t omin = -.5;
106     Double_t omax = 7.5;
107     Int_t    o    = 8;
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}");
112   }
113   //__________________________________________________________________
114   /** Begining of event
115       @param ev Event number
116       @return @c false on error */
117   Bool_t Begin(Int_t ev) 
118   {
119     fMap.Reset();
120     return AliFMDInput::Begin(ev);
121   }
122   //__________________________________________________________________
123   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
124   {
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();
131     if (str > 511) {
132       AliWarning(Form("Bad strip number %d in hit", str));
133       return kTRUE;
134     }
135     if (fPrimary) {
136       TParticle* kine = fStack->Particle(hit->Track());
137       if (!kine) return kTRUE;
138       if (!kine->IsPrimary()) return kTRUE;
139     }
140     
141     fMap(det, rng, sec, str).fEdep += hit->Edep();
142     fMap(det, rng, sec, str).fN++;
143     return kTRUE;
144   }
145
146   //__________________________________________________________________
147   Bool_t ProcessDigit(AliFMDDigit* digit)
148   {
149     if (!digit) return kTRUE;
150
151     UShort_t det = digit->Detector();
152     Char_t   rng = digit->Ring();
153     UShort_t sec = digit->Sector();
154     UShort_t str = digit->Strip();
155     if (str > 511) {
156       AliWarning(Form("Bad strip number %d in digit", str));
157       return kTRUE;
158     }
159     Double_t edep = fMap(det, rng, sec, str).fEdep;
160     if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
161       
162     return kTRUE;
163   }
164
165   //__________________________________________________________________
166   Bool_t ProcessRecPoint(AliFMDRecPoint* single) 
167   {
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();
173     if (str > 511) {
174       AliWarning(Form("Bad strip number %d in single", str));
175       return kTRUE;
176     }
177     Double_t edep = fMap(det, rng, sec, str).fEdep;
178     Int_t    nhit = fMap(det, rng, sec, str).fN;
179     if (edep > 0) {
180       fHitEvsRecM->Fill(edep, single->Particles());
181       fHitEvsRecE->Fill(edep, single->Edep());
182       fDiffE->Fill((single->Edep() - edep) / edep);
183     }
184     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
185     return kTRUE;
186   }
187   //__________________________________________________________________
188   Bool_t Finish()
189   {
190     gStyle->SetPalette(1);
191     gStyle->SetOptTitle(0);
192     gStyle->SetCanvasColor(0);
193     gStyle->SetCanvasBorderSize(0);
194     gStyle->SetPadColor(0);
195     gStyle->SetPadBorderSize(0);
196
197     new TCanvas("c0", fHitEvsAdc->GetTitle());
198     fHitEvsAdc->SetStats(kFALSE);
199     fHitEvsAdc->Draw("COLZ");
200
201     new TCanvas("c1", fHitEvsRecM->GetTitle());
202     fHitEvsRecM->SetStats(kFALSE);
203     fHitEvsRecM->Draw("COLZ");
204
205     new TCanvas("c2", fHitEvsRecE->GetTitle());
206     fHitEvsRecE->SetStats(kFALSE);
207     fHitEvsRecE->Draw("COLZ");
208
209     new TCanvas("c3", fDiffE->GetTitle());
210     fDiffE->Draw();
211
212     new TCanvas("c4", fHitsVsRecM->GetTitle());
213     fHitsVsRecM->SetStats(kFALSE);
214     fHitsVsRecM->Draw("COLZ");
215
216     return kTRUE;
217   }
218
219   ClassDef(DrawHitsRecs,0);
220 };
221
222 //____________________________________________________________________
223 //
224 // EOF
225 //