Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
[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 #include <TF1.h>
31
32 //____________________________________________________________________
33 /** @class DrawHitsRecs
34     @brief Draw hit energy loss versus rec point mult
35     @code 
36     Root> .L Compile.C
37     Root> Compile("DrawHitsRecs.C")
38     Root> DrawHitsRecs c
39     Root> c.Run();
40     @endcode
41     @ingroup FMD_script
42  */
43 class DrawHitsRecs : public AliFMDInput
44 {
45 private:
46   TH2D* fHitEvsAdc;
47   TH2D* fHitEvsRecM;  // Histogram 
48   TH2D* fHitEvsRecE;  // Histogram 
49   TH1D* fDiffE;       // Histogram 
50   TH2D* fHitsVsRecM;  // Histogram 
51   TH2D* fDiffM;       // Histogram 
52   TH1*  fHitEloss;
53   TH1*  fRecEloss;
54   AliFMDEdepMap  fMap;
55   AliFMDFloatMap fEta;
56   AliFMDFloatMap fPhi;
57   AliFMDFloatMap fMult;
58   Bool_t fPrimary;
59 public:
60   //__________________________________________________________________
61   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
62   {
63     TArrayF bins(n+1);
64     Float_t dp   = n / TMath::Log10(max / min);
65     Float_t pmin = TMath::Log10(min);
66     bins[0]      = min;
67     for (Int_t i = 1; i < n+1; i++) {
68       Float_t p = pmin + i / dp;
69       bins[i]   = TMath::Power(10, p);
70     }
71     return bins;
72   }
73   //__________________________________________________________________
74   DrawHitsRecs(Bool_t primary=kFALSE,
75                Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
76                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5) 
77   { 
78     fPrimary = primary;
79     AddLoad(kRecPoints);
80     AddLoad(kHits);
81     AddLoad(kDigits);
82     if (fPrimary) AddLoad(kKinematics);
83     TArrayF eloss(MakeLogScale(n, emin, emax));
84     TArrayF mults(m+1);
85     mults[0] = mmin;
86     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
87
88     fHitEvsAdc  = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
89                            n, emin, emax, 1025, -.5, 1024.5);
90     fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
91     fHitEvsAdc->SetYTitle("ADC");
92     
93     fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}", 
94                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
95     fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
96     fHitEvsRecM->SetYTitle("M_{rec}");
97
98     fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}", 
99                             n, emin, emax, n, emin, emax);
100     fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
101     fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
102
103
104     fDiffE = new TH1D("diffE", 
105                       "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}", 
106                       1100, -1, 1.1);
107     fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
108     
109     Double_t omin = -.5;
110     Double_t omax = 7.5;
111     Int_t    o    = 8;
112     fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
113                            o, omin, omax, m, mmin, mmax);
114     fHitsVsRecM->SetXTitle("# of Hits");
115     fHitsVsRecM->SetYTitle("M_{rec}");
116
117     fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}", 
118                       41, -20.5, 20.5, 70, 1.5, 5);
119     // 36, -TMath::Pi(),TMath::Pi());
120     fDiffM->SetXTitle("M_{sim} - M_{rec}");
121     fDiffM->SetYTitle("|#eta|");
122     // fDiffM->SetYTitle("Detector");
123
124     fHitEloss = new TH1D("hitEloss", "#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)", 
125                          100, 0, 10);
126     fHitEloss->SetFillColor(2);
127     fHitEloss->SetFillStyle(3001);
128     
129     fRecEloss = new TH1D("recEloss", "#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)", 
130                          100, 0, 10);
131     fRecEloss->SetFillColor(4);
132     fRecEloss->SetFillStyle(3001);
133   }
134   //__________________________________________________________________
135   /** Begining of event
136       @param ev Event number
137       @return @c false on error */
138   Bool_t Begin(Int_t ev) 
139   {
140     fMap.Reset();
141     return AliFMDInput::Begin(ev);
142   }
143   //__________________________________________________________________
144   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
145   {
146     // Cache the energy loss 
147     if (!hit) return kFALSE;
148     UShort_t det = hit->Detector();
149     Char_t   rng = hit->Ring();
150     UShort_t sec = hit->Sector();
151     UShort_t str = hit->Strip();
152     if (str > 511) {
153       AliWarning(Form("Bad strip number %d in hit", str));
154       return kTRUE;
155     }
156     if (fPrimary) {
157       TParticle* kine = fStack->Particle(hit->Track());
158       if (!kine) return kTRUE;
159       if (!kine->IsPrimary()) return kTRUE;
160     }
161     
162     if (hit->Edep()/hit->Length() > 0.1) fHitEloss->Fill(hit->Edep() / hit->Length());
163     fMap(det, rng, sec, str).fEdep += hit->Edep();
164     fMap(det, rng, sec, str).fN++;
165     return kTRUE;
166   }
167
168   //__________________________________________________________________
169   Bool_t ProcessDigit(AliFMDDigit* digit)
170   {
171     if (!digit) return kTRUE;
172
173     UShort_t det = digit->Detector();
174     Char_t   rng = digit->Ring();
175     UShort_t sec = digit->Sector();
176     UShort_t str = digit->Strip();
177     if (str > 511) {
178       AliWarning(Form("Bad strip number %d in digit", str));
179       return kTRUE;
180     }
181     Double_t edep = fMap(det, rng, sec, str).fEdep;
182     if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
183       
184     return kTRUE;
185   }
186
187   //__________________________________________________________________
188   Bool_t ProcessRecPoint(AliFMDRecPoint* single) 
189   {
190     if (!single) return kTRUE;
191     UShort_t det = single->Detector();
192     Char_t   rng = single->Ring();
193     UShort_t sec = single->Sector();
194     UShort_t str = single->Strip();
195     if (str > 511) {
196       AliWarning(Form("Bad strip number %d in single", str));
197       return kTRUE;
198     }
199     Double_t edep = fMap(det, rng, sec, str).fEdep;
200     Int_t    nhit = fMap(det, rng, sec, str).fN;
201     if (edep > 0) {
202       fHitEvsRecM->Fill(edep, single->Particles());
203       fHitEvsRecE->Fill(edep, single->Edep());
204       fDiffE->Fill((single->Edep() - edep) / edep);
205     }
206     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
207     fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
208     if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
209     return kTRUE;
210   }
211   //__________________________________________________________________
212   Bool_t Finish()
213   {
214     gStyle->SetPalette(1);
215     gStyle->SetOptTitle(0);
216     gStyle->SetCanvasColor(0);
217     gStyle->SetCanvasBorderSize(0);
218     gStyle->SetPadColor(0);
219     gStyle->SetPadBorderSize(0);
220     TCanvas* c = 0;
221
222     c = new TCanvas("c0", fHitEvsAdc->GetTitle());
223     fHitEvsAdc->SetStats(kFALSE);
224     fHitEvsAdc->Draw("COLZ");
225
226     c = new TCanvas("c1", fHitEvsRecM->GetTitle());
227     fHitEvsRecM->SetStats(kFALSE);
228     fHitEvsRecM->Draw("COLZ");
229
230     c = new TCanvas("c2", fHitEvsRecE->GetTitle());
231     fHitEvsRecE->SetStats(kFALSE);
232     fHitEvsRecE->Draw("COLZ");
233
234     c = new TCanvas("c3", fDiffE->GetTitle());
235     c->SetLogz();
236     fDiffE->Draw();
237
238     c = new TCanvas("c4", fHitsVsRecM->GetTitle());
239     c->SetLogz();
240     fHitsVsRecM->SetStats(kFALSE);
241     fHitsVsRecM->Draw("COLZ");
242
243     c = new TCanvas("c5", fDiffM->GetTitle());
244     fDiffM->SetFillColor(2);
245     fDiffM->SetFillStyle(3001);
246     c->SetLogz();
247     fDiffM->Draw("colz");
248
249     c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
250     fRecEloss->Scale(1./fRecEloss->GetMaximum());
251     fRecEloss->Draw();
252     fRecEloss->Fit("landau", "", "SAME", 2, 4);
253     TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
254     fHitEloss->Scale(1./fHitEloss->GetMaximum());
255     fHitEloss->Draw("same");
256     fHitEloss->Fit("landau", "", "SAME", 2, 10);
257     TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
258     std::cout << "Hit MPV,width: " << hitResp->GetParameter(1) << ","
259               << hitResp->GetParameter(2) << "\n"
260               << "Rec MPV,width: " << recResp->GetParameter(1) << ","
261               << recResp->GetParameter(2) << std::endl;
262     c->SetLogy();
263
264     return kTRUE;
265   }
266
267   ClassDef(DrawHitsRecs,0);
268 };
269
270 //____________________________________________________________________
271 //
272 // EOF
273 //