Use VMC id's rather than TGeo id's
[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   TH2D* fDiffM;       // Histogram 
51   AliFMDEdepMap  fMap;
52   AliFMDFloatMap fEta;
53   AliFMDFloatMap fPhi;
54   AliFMDFloatMap fMult;
55   Bool_t fPrimary;
56 public:
57   //__________________________________________________________________
58   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
59   {
60     TArrayF bins(n+1);
61     Float_t dp   = n / TMath::Log10(max / min);
62     Float_t pmin = TMath::Log10(min);
63     bins[0]      = min;
64     for (Int_t i = 1; i < n+1; i++) {
65       Float_t p = pmin + i / dp;
66       bins[i]   = TMath::Power(10, p);
67     }
68     return bins;
69   }
70   //__________________________________________________________________
71   DrawHitsRecs(Bool_t primary=kFALSE,
72                Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
73                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5) 
74   { 
75     fPrimary = primary;
76     AddLoad(kRecPoints);
77     AddLoad(kHits);
78     AddLoad(kDigits);
79     if (fPrimary) AddLoad(kKinematics);
80     TArrayF eloss(MakeLogScale(n, emin, emax));
81     TArrayF mults(m+1);
82     mults[0] = mmin;
83     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
84
85     fHitEvsAdc  = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
86                            n, emin, emax, 1025, -.5, 1024.5);
87     fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
88     fHitEvsAdc->SetYTitle("ADC");
89     
90     fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}", 
91                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
92     fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
93     fHitEvsRecM->SetYTitle("M_{rec}");
94
95     fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}", 
96                             n, emin, emax, n, emin, emax);
97     fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
98     fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
99
100
101     fDiffE = new TH1D("diffE", 
102                       "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}", 
103                       1100, -1, 1.1);
104     fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
105     
106     Double_t omin = -.5;
107     Double_t omax = 7.5;
108     Int_t    o    = 8;
109     fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
110                            o, omin, omax, m, mmin, mmax);
111     fHitsVsRecM->SetXTitle("# of Hits");
112     fHitsVsRecM->SetYTitle("M_{rec}");
113
114     fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}", 
115                       41, -20.5, 20.5, 70, 1.5, 5);
116     // 36, -TMath::Pi(),TMath::Pi());
117     fDiffM->SetXTitle("M_{sim} - M_{rec}");
118     fDiffM->SetYTitle("|#eta|");
119     // fDiffM->SetYTitle("Detector");
120     
121   }
122   //__________________________________________________________________
123   /** Begining of event
124       @param ev Event number
125       @return @c false on error */
126   Bool_t Begin(Int_t ev) 
127   {
128     fMap.Reset();
129     return AliFMDInput::Begin(ev);
130   }
131   //__________________________________________________________________
132   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
133   {
134     // Cache the energy loss 
135     if (!hit) return kFALSE;
136     UShort_t det = hit->Detector();
137     Char_t   rng = hit->Ring();
138     UShort_t sec = hit->Sector();
139     UShort_t str = hit->Strip();
140     if (str > 511) {
141       AliWarning(Form("Bad strip number %d in hit", str));
142       return kTRUE;
143     }
144     if (fPrimary) {
145       TParticle* kine = fStack->Particle(hit->Track());
146       if (!kine) return kTRUE;
147       if (!kine->IsPrimary()) return kTRUE;
148     }
149     
150     fMap(det, rng, sec, str).fEdep += hit->Edep();
151     fMap(det, rng, sec, str).fN++;
152     return kTRUE;
153   }
154
155   //__________________________________________________________________
156   Bool_t ProcessDigit(AliFMDDigit* digit)
157   {
158     if (!digit) return kTRUE;
159
160     UShort_t det = digit->Detector();
161     Char_t   rng = digit->Ring();
162     UShort_t sec = digit->Sector();
163     UShort_t str = digit->Strip();
164     if (str > 511) {
165       AliWarning(Form("Bad strip number %d in digit", str));
166       return kTRUE;
167     }
168     Double_t edep = fMap(det, rng, sec, str).fEdep;
169     if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
170       
171     return kTRUE;
172   }
173
174   //__________________________________________________________________
175   Bool_t ProcessRecPoint(AliFMDRecPoint* single) 
176   {
177     if (!single) return kTRUE;
178     UShort_t det = single->Detector();
179     Char_t   rng = single->Ring();
180     UShort_t sec = single->Sector();
181     UShort_t str = single->Strip();
182     if (str > 511) {
183       AliWarning(Form("Bad strip number %d in single", str));
184       return kTRUE;
185     }
186     Double_t edep = fMap(det, rng, sec, str).fEdep;
187     Int_t    nhit = fMap(det, rng, sec, str).fN;
188     if (edep > 0) {
189       fHitEvsRecM->Fill(edep, single->Particles());
190       fHitEvsRecE->Fill(edep, single->Edep());
191       fDiffE->Fill((single->Edep() - edep) / edep);
192     }
193     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
194     fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
195     return kTRUE;
196   }
197   //__________________________________________________________________
198   Bool_t Finish()
199   {
200     gStyle->SetPalette(1);
201     gStyle->SetOptTitle(0);
202     gStyle->SetCanvasColor(0);
203     gStyle->SetCanvasBorderSize(0);
204     gStyle->SetPadColor(0);
205     gStyle->SetPadBorderSize(0);
206     TCanvas* c = 0;
207
208     c = new TCanvas("c0", fHitEvsAdc->GetTitle());
209     fHitEvsAdc->SetStats(kFALSE);
210     fHitEvsAdc->Draw("COLZ");
211
212     c = new TCanvas("c1", fHitEvsRecM->GetTitle());
213     fHitEvsRecM->SetStats(kFALSE);
214     fHitEvsRecM->Draw("COLZ");
215
216     c = new TCanvas("c2", fHitEvsRecE->GetTitle());
217     fHitEvsRecE->SetStats(kFALSE);
218     fHitEvsRecE->Draw("COLZ");
219
220     c = new TCanvas("c3", fDiffE->GetTitle());
221     c->SetLogz();
222     fDiffE->Draw();
223
224     c = new TCanvas("c4", fHitsVsRecM->GetTitle());
225     c->SetLogz();
226     fHitsVsRecM->SetStats(kFALSE);
227     fHitsVsRecM->Draw("COLZ");
228
229     c = new TCanvas("c5", fDiffM->GetTitle());
230     fDiffM->SetFillColor(2);
231     fDiffM->SetFillStyle(3001);
232     c->SetLogz();
233     fDiffM->Draw("colz");
234
235
236     return kTRUE;
237   }
238
239   ClassDef(DrawHitsRecs,0);
240 };
241
242 //____________________________________________________________________
243 //
244 // EOF
245 //