f12cc89c96368c60e4d7b348ae4b2de78800b781
[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 <AliFMDMultStrip.h>
17 #include <AliFMDMultRegion.h>
18 #include <AliFMDDigit.h>
19 #include <AliFMDInput.h>
20 #include <AliFMDEdepMap.h>
21 #include <AliFMDFloatMap.h>
22 #include <iostream>
23 #include <TStyle.h>
24 #include <TArrayF.h>
25 #include <TCanvas.h>
26 #include <TParticle.h>
27 #include <TClonesArray.h>
28 #include <TTree.h>
29 #include <AliStack.h>
30 #include <AliLog.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 AliFMDInputHits
44 {
45 private:
46   TH2D* fElossVsMult; // Histogram 
47   TH2D* fHitsVsStrip;  // Histogram 
48   TH2D* fHitsVsRegion;  // Histogram 
49   AliFMDEdepMap fMap;
50   AliFMDFloatMap fEta;
51   AliFMDFloatMap fPhi;
52   AliFMDFloatMap fMult;
53   Bool_t fPrimary;
54 public:
55   //__________________________________________________________________
56   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
57   {
58     TArrayF bins(n+1);
59     Float_t dp   = n / TMath::Log10(max / min);
60     Float_t pmin = TMath::Log10(min);
61     bins[0]      = min;
62     for (Int_t i = 1; i < n+1; i++) {
63       Float_t p = pmin + i / dp;
64       bins[i]   = TMath::Power(10, p);
65     }
66     return bins;
67   }
68   //__________________________________________________________________
69   DrawHitsRecs(Bool_t primary=kFALSE,
70                Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
71                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5) 
72   { 
73     fPrimary = primary;
74     AddLoad(kRecPoints);
75     AddLoad(kHits);
76     if (fPrimary) AddLoad(kKinematics);
77     TArrayF eloss(MakeLogScale(n, emin, emax));
78     TArrayF mults(m+1);
79     mults[0] = mmin;
80     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
81     fElossVsMult = new TH2D("elossVsMult", 
82                             "#Delta E vs. Multiplicity (single)", 
83                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
84     fElossVsMult->SetXTitle("#Delta E/#Delta x [MeV/cm]");
85     fElossVsMult->SetYTitle("Strip Multiplicity");
86
87     Double_t omin = -.5;
88     Double_t omax = 7.5;
89     Int_t    o    = 8;
90     fHitsVsStrip = new TH2D("hitsVsStrip", 
91                             "# of Hits vs. Multiplicity (strip)",
92                            o, omin, omax, m, mmin, mmax);
93     fHitsVsStrip->SetXTitle("# of Hits");
94     fHitsVsStrip->SetYTitle("Strip Multiplicity");
95   }
96   //__________________________________________________________________
97   /** Begining of event
98       @param ev Event number
99       @return @c false on error */
100   Bool_t Begin(Int_t ev) 
101   {
102     fMap.Reset();
103     return AliFMDInputHits::Begin(ev);
104   }
105   //__________________________________________________________________
106   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
107   {
108     // Cache the energy loss 
109     if (!hit) return kFALSE;
110     UShort_t det = hit->Detector();
111     Char_t   rng = hit->Ring();
112     UShort_t sec = hit->Sector();
113     UShort_t str = hit->Strip();
114     if (str > 511) {
115       AliWarning(Form("Bad strip number %d in hit", str));
116       return kTRUE;
117     }
118     if (fPrimary) {
119       TParticle* kine = fStack->Particle(hit->Track());
120       if (!kine) return kTRUE;
121       if (!kine->IsPrimary()) return kTRUE;
122     }
123     
124     fMap(det, rng, sec, str).fEdep += hit->Edep();
125     fMap(det, rng, sec, str).fN++;
126     return kTRUE;
127   }
128   //__________________________________________________________________
129   Bool_t ProcessRecPoint(AliFMDRecPoint* single) 
130   {
131     if (!single) continue;
132     UShort_t det = single->Detector();
133     Char_t   rng = single->Ring();
134     UShort_t sec = single->Sector();
135     UShort_t str = single->Strip();
136     if (str > 511) {
137       AliWarning(Form("Bad strip number %d in single", str));
138       continue;
139     }
140     if (fMap(det, rng, sec, str).fEdep > 0) 
141       fElossVsMult->Fill(fMap(det, rng, sec, str).fEdep, 
142                          single->Particles());
143     if (fMap(det, rng, sec, str).fN > 0) 
144       fHitsVsStrip->Fill(fMap(det, rng, sec, str).fN, 
145                          single->Particles());
146     return kTRUE;
147   }
148   //__________________________________________________________________
149   Bool_t Finish()
150   {
151     gStyle->SetPalette(1);
152     gStyle->SetOptTitle(0);
153     gStyle->SetCanvasColor(0);
154     gStyle->SetCanvasBorderSize(0);
155     gStyle->SetPadColor(0);
156     gStyle->SetPadBorderSize(0);
157
158     new TCanvas("c1", "Energy loss vs. Strip Multiplicity");
159     fElossVsMult->SetStats(kFALSE);
160     fElossVsMult->Draw("COLZ");
161
162     new TCanvas("c2", "# of Hits vs. Strip Multiplicity");
163     fHitsVsStrip->SetStats(kFALSE);
164     fHitsVsStrip->Draw("COLZ");
165
166     return kTRUE;
167   }
168
169   ClassDef(DrawHitsRecs,0);
170 };
171
172 //____________________________________________________________________
173 //
174 // EOF
175 //