Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxCalib.cxx
1 #include "AliHighPtDeDxCalib.h"
2
3 #ifndef __IOSTREAM__
4 #include <iostream>
5 #endif
6
7 using namespace std;
8
9 ClassImp(AliHighPtDeDxCalib);
10
11 //
12 // AliHighPtDeDxCalib class
13 //
14 // This class contains the AliHighPtDeDxCalib information 
15 //
16
17 //_________________________________________________________
18 AliHighPtDeDxCalib::AliHighPtDeDxCalib():
19   AliHighPtDeDxBase(),
20   fStep(1),             // Step 1 = Eta calibration, step 2 = dE/dx calibration, step = 3 ncl calib
21   fInit(0),             // Step 1 = Eta calibration, step 2 = dE/dx calibration, step = 3 ncl calib
22   fPMIPMin(0.4),        // Min P for MIP pion
23   fPMIPMax(0.6),        // Max P for MIP pion
24   fDeDxMIPMin(40),      // Min dE/dx for MIP pion
25   fDeDxMIPMax(60),      // Max dE/dx for MIP pion
26   fDeltaBeta(0.1),      // delta beta cut for electrons 
27   fDeDxPi(0x0),         // dE/dx vs p for pions
28   fSigmaDeDx(0x0),      // sigma dE/dx vs ncl
29   fDeDxVsEtaNeg(0x0),   // eta < 0 dE/dx calib
30   fDeDxVsEtaPos(0x0),   // eta > 0 dE/dx calib
31   fDeDxVsNcl(0x0),      // ncl dE/dx calib
32   hSelection1(0x0),             // selected region in p and dE/dx for pion MIPs
33   hSelection2(0x0),             // selected region in p and dE/dx for pion MIPs
34   hSelection3(0x0),             // selected region in p and dE/dx for pion MIPs
35   hSelectionElectrons2(0x0),    // selected region in p and dE/dx for electrons
36   hSelectionElectrons3(0x0),    // selected region in p and dE/dx for electrons
37   hDeDxVsEta(0x0),             // dE/dx vs eta uncalibrated (pion MIPs)
38   hDeDxVsEtaElectrons(0x0),    // dE/dx vs eta uncalibrated (electrons)
39   hNclVsEta(0x0),              // Ncl vs eta (pion MIPs)
40   hNclVsEtaElectrons(0x0),     // Ncl vs eta (electrons)
41   hDeDxVsEtaCal(0x0),          // dE/dx vs eta calibrated (pion MIPs)
42   hDeDxVsEtaCalElectrons(0x0), // dE/dx vs eta calibrated (electrons)
43   hMeanEta(0x0),               // <eta> in the 4 eta interval (pion MIPs)
44   hMeanEtaElectrons(0x0),      // <eta> in the 4 eta interval (electrons)
45   hDeDx(0x0),                  // dE/dx no eta cut    (pion MIPs)
46   hDeDx1(0x0),                 // dE/dx 0.0<|eta|<0.2 (pion MIPs)
47   hDeDx2(0x0),                 // dE/dx 0.2<|eta|<0.4 (pion MIPs)
48   hDeDx3(0x0),                 // dE/dx 0.4<|eta|<0.6 (pion MIPs)
49   hDeDx4(0x0),                 // dE/dx 0.6<|eta|<0.8 (pion MIPs)
50   hDeDxElectrons(0x0),         // dE/dx no eta cut    (electrons)
51   hDeDxElectrons1(0x0),        // dE/dx 0.0<|eta|<0.2 (electrons)
52   hDeDxElectrons2(0x0),        // dE/dx 0.2<|eta|<0.4 (electrons)
53   hDeDxElectrons3(0x0),        // dE/dx 0.4<|eta|<0.6 (electrons)
54   hDeDxElectrons4(0x0),        // dE/dx 0.6<|eta|<0.8 (electrons)
55   hDeDxVsNclBefore(0x0),       // dE/dx vs ncl for step 2 calib (pion MIPs)
56   hDeDxVsNcl(0x0),             // dE/dx vs ncl no eta cut    (pion MIPs)
57   hDeDxVsNcl1(0x0),            // dE/dx vs ncl 0.0<|eta|<0.2 (pion MIPs)
58   hDeDxVsNcl2(0x0),            // dE/dx vs ncl 0.2<|eta|<0.4 (pion MIPs)
59   hDeDxVsNcl3(0x0),            // dE/dx vs ncl 0.4<|eta|<0.6 (pion MIPs)
60   hDeDxVsNcl4(0x0),            // dE/dx vs ncl 0.6<|eta|<0.8 (pion MIPs)
61   hDeDxVsNclElectrons(0x0),    // dE/dx vs ncl no eta cut    (electrons)
62   hDeDxVsNclElectrons1(0x0),   // dE/dx vs ncl 0.0<|eta|<0.2 (electrons)
63   hDeDxVsNclElectrons2(0x0),   // dE/dx vs ncl 0.2<|eta|<0.4 (electrons)
64   hDeDxVsNclElectrons3(0x0),   // dE/dx vs ncl 0.4<|eta|<0.6 (electrons)
65   hDeDxVsNclElectrons4(0x0),   // dE/dx vs ncl 0.6<|eta|<0.8 (electrons)
66   hMeanP(0x0),        // <p> vs p
67   hDeDxVsP(0x0),      // dE/dx vs p 
68   hDeDxVsPPiMc(0x0),  // dE/dx vs p for MC pions
69   hDeDxVsPKMc(0x0),   // dE/dx vs p for MC Kaons
70   hDeDxVsPPMc(0x0),   // dE/dx vs p for MC protons
71   hDeDxVsPEMc(0x0),   // dE/dx vs p for MC electrons
72   hDeDxVsPMuMc(0x0)   // dE/dx vs p for MC muons
73 {
74   // default constructor - do not use
75
76 }
77
78 //_________________________________________________________
79 AliHighPtDeDxCalib::AliHighPtDeDxCalib(const char* name, const char* title):
80   AliHighPtDeDxBase(name, title),
81   fStep(1),             // Step 1 = Eta calibration, step 2 = dE/dx calibration
82   fInit(0),             // Step 1 = Eta calibration, step 2 = dE/dx calibration, step = 3 ncl calib
83   fPMIPMin(0.4),        // Min P for MIP pion
84   fPMIPMax(0.6),        // Max P for MIP pion
85   fDeDxMIPMin(40),      // Min dE/dx for MIP pion
86   fDeDxMIPMax(60),      // Max dE/dx for MIP pion
87   fDeltaBeta(0.1),      // delta beta cut for electrons 
88   fDeDxPi(0x0),         // dE/dx vs p for pions
89   fSigmaDeDx(0x0),      // sigma dE/dx vs ncl
90   fDeDxVsEtaNeg(0x0),   // eta < 0 dE/dx calib
91   fDeDxVsEtaPos(0x0),   // eta > 0 dE/dx calib
92   fDeDxVsNcl(0x0),      // ncl dE/dx calib
93   hSelection1(0x0),             // selected region in p and dE/dx for pion MIPs
94   hSelection2(0x0),             // selected region in p and dE/dx for pion MIPs
95   hSelection3(0x0),             // selected region in p and dE/dx for pion MIPs
96   hSelectionElectrons2(0x0),    // selected region in p and dE/dx for electrons
97   hSelectionElectrons3(0x0),    // selected region in p and dE/dx for electrons
98   hDeDxVsEta(0x0),             // dE/dx vs eta uncalibrated (pion MIPs)
99   hDeDxVsEtaElectrons(0x0),    // dE/dx vs eta uncalibrated (electrons)
100   hNclVsEta(0x0),              // Ncl vs eta (pion MIPs)
101   hNclVsEtaElectrons(0x0),     // Ncl vs eta (electrons)
102   hDeDxVsEtaCal(0x0),          // dE/dx vs eta calibrated (pion MIPs)
103   hDeDxVsEtaCalElectrons(0x0), // dE/dx vs eta calibrated (electrons)
104   hMeanEta(0x0),               // <eta> in the 4 eta interval (pion MIPs)
105   hMeanEtaElectrons(0x0),      // <eta> in the 4 eta interval (electrons)
106   hDeDx(0x0),                  // dE/dx no eta cut    (pion MIPs)
107   hDeDx1(0x0),                 // dE/dx 0.0<|eta|<0.2 (pion MIPs)
108   hDeDx2(0x0),                 // dE/dx 0.2<|eta|<0.4 (pion MIPs)
109   hDeDx3(0x0),                 // dE/dx 0.4<|eta|<0.6 (pion MIPs)
110   hDeDx4(0x0),                 // dE/dx 0.6<|eta|<0.8 (pion MIPs)
111   hDeDxElectrons(0x0),         // dE/dx no eta cut    (electrons)
112   hDeDxElectrons1(0x0),        // dE/dx 0.0<|eta|<0.2 (electrons)
113   hDeDxElectrons2(0x0),        // dE/dx 0.2<|eta|<0.4 (electrons)
114   hDeDxElectrons3(0x0),        // dE/dx 0.4<|eta|<0.6 (electrons)
115   hDeDxElectrons4(0x0),        // dE/dx 0.6<|eta|<0.8 (electrons)
116   hDeDxVsNclBefore(0x0),       // dE/dx vs ncl no eta cut    (pion MIPs)
117   hDeDxVsNcl(0x0),             // dE/dx vs ncl no eta cut    (pion MIPs)
118   hDeDxVsNcl1(0x0),            // dE/dx vs ncl 0.0<|eta|<0.2 (pion MIPs)
119   hDeDxVsNcl2(0x0),            // dE/dx vs ncl 0.2<|eta|<0.4 (pion MIPs)
120   hDeDxVsNcl3(0x0),            // dE/dx vs ncl 0.4<|eta|<0.6 (pion MIPs)
121   hDeDxVsNcl4(0x0),            // dE/dx vs ncl 0.6<|eta|<0.8 (pion MIPs)
122   hDeDxVsNclElectrons(0x0),    // dE/dx vs ncl no eta cut    (electrons)
123   hDeDxVsNclElectrons1(0x0),   // dE/dx vs ncl 0.0<|eta|<0.2 (electrons)
124   hDeDxVsNclElectrons2(0x0),   // dE/dx vs ncl 0.2<|eta|<0.4 (electrons)
125   hDeDxVsNclElectrons3(0x0),   // dE/dx vs ncl 0.4<|eta|<0.6 (electrons)
126   hDeDxVsNclElectrons4(0x0),   // dE/dx vs ncl 0.6<|eta|<0.8 (electrons)
127   hMeanP(0x0),        // <p> vs p
128   hDeDxVsP(0x0),      // dE/dx vs p 
129   hDeDxVsPPiMc(0x0),  // dE/dx vs p for MC pions
130   hDeDxVsPKMc(0x0),   // dE/dx vs p for MC Kaons
131   hDeDxVsPPMc(0x0),   // dE/dx vs p for MC protons
132   hDeDxVsPEMc(0x0),   // dE/dx vs p for MC electrons
133   hDeDxVsPMuMc(0x0)   // dE/dx vs p for MC muons
134 {
135   // named constructor
136 }
137
138 //_________________________________________________________
139 AliHighPtDeDxCalib::~AliHighPtDeDxCalib()
140 {
141   delete fDeDxPi;
142   delete fSigmaDeDx;
143   delete fDeDxVsEtaNeg;
144   delete fDeDxVsEtaPos;
145   delete fDeDxVsNcl;
146   delete hSelection1;
147   delete hSelection2;
148   delete hSelection3;
149   delete hSelectionElectrons2;
150   delete hSelectionElectrons3;
151   delete hDeDxVsEta;
152   delete hDeDxVsEtaElectrons;
153   delete hNclVsEta;
154   delete hNclVsEtaElectrons;
155   delete hDeDxVsEtaCal;
156   delete hDeDxVsEtaCalElectrons;
157   delete hMeanEta;
158   delete hMeanEtaElectrons;
159   delete hDeDx;
160   delete hDeDx1;
161   delete hDeDx2;
162   delete hDeDx3;
163   delete hDeDx4;
164   delete hDeDxElectrons;
165   delete hDeDxElectrons1;
166   delete hDeDxElectrons2;
167   delete hDeDxElectrons3;
168   delete hDeDxElectrons4;
169   delete hDeDxVsNclBefore;
170   delete hDeDxVsNcl;
171   delete hDeDxVsNcl1;
172   delete hDeDxVsNcl2;
173   delete hDeDxVsNcl3;
174   delete hDeDxVsNcl4;
175   delete hDeDxVsNclElectrons;
176   delete hDeDxVsNclElectrons1;
177   delete hDeDxVsNclElectrons2;
178   delete hDeDxVsNclElectrons3;
179   delete hDeDxVsNclElectrons4;
180   delete hMeanP;
181   delete hDeDxVsP;
182   delete hDeDxVsPPiMc;
183   delete hDeDxVsPKMc;
184   delete hDeDxVsPPMc;
185   delete hDeDxVsPEMc;
186   delete hDeDxVsPMuMc;
187 }
188
189 //_________________________________________________________
190 void AliHighPtDeDxCalib::Init(Int_t nPtBins, Double_t* ptBins)
191 {
192   //
193   // Create histograms and functions
194   //
195
196   //
197   // init base class
198   //
199   AliHighPtDeDxBase::Init(nPtBins, ptBins);
200
201   //
202   // functions
203   //
204   fDeDxVsEtaNeg = new TF1("fDeDxVsEtaNeg", "pol3", -1,  0);
205   fDeDxVsEtaNeg->SetParameters(50, 0, 0, 0);
206   fDeDxVsEtaPos = new TF1("fDeDxVsEtaPos", "pol3",  0, +1);
207   fDeDxVsEtaPos->SetParameters(50, 0, 0, 0);
208
209   fDeDxVsNcl = new TF1("fDeDxVsNcl", "pol1",  0, 160);
210   fDeDxVsNcl->SetParameters(50, 0);
211 }
212
213 void AliHighPtDeDxCalib::Init(Int_t step, Int_t nPtBins, Double_t* ptBins)
214 {
215   if(step != fStep)
216     cout << "Warning! Innit called for step " << step << " but next step is " << fStep << endl;
217
218   if(fInit>=step) {
219     cout << "Step " << step << " has been initialized!" << endl;
220     return;
221   }
222
223   switch (step) {
224   case 1:
225     //
226     // step 1 histograms - nCl calibration
227     //
228     hSelection1 = new TH2D("hSelection1", "dE/dx vs P (pion MIPs) step 1; P [GeV/c]; dE/dx",
229                            100, 0.2, 0.8, 
230                            100, 0, 100);
231     hSelection1->SetDirectory(0);
232
233     hDeDxVsNclBefore = new TH2D("hDeDxVsNclBefore", "dE/dx vs Ncl (pion MIPs); Ncl; dE/dx",
234                                 18, 69.5, 159.5, 
235                                 fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
236     hDeDxVsNclBefore->Sumw2();
237     hDeDxVsNclBefore->SetDirectory(0);
238
239     fInit = 1;
240     
241     break;
242
243   case 2:
244     //
245     // step 2 histograms - eta calibration
246     //
247     
248     hSelection2 = new TH2D("hSelection1", "dE/dx vs P (pion MIPs) step 2; P [GeV/c]; dE/dx",
249                            100, 0.2, 0.8, 
250                            100, 0, 100);
251     hSelection2->SetDirectory(0);
252
253     hSelectionElectrons2 = new TH2D("hSelectionElectrons3", "dE/dx vs P (electrons) step 2; P [GeV/c]; dE/dx",
254                                    100, 0.2, 0.8, 
255                                    100, 0, 100);
256    hSelectionElectrons2->SetDirectory(0);
257
258     hDeDxVsEta = new TH2D("hDeDxVsEta", "dE/dx vs #eta (pion MIPs); #eta; dE/dx",
259                           50, -1.0, 1.0, 
260                           fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
261     hDeDxVsEta->Sumw2();
262     hDeDxVsEta->SetDirectory(0);
263     
264     hDeDxVsEtaElectrons = new TH2D("hDeDxVsEtaElectrons", "dE/dx vs #eta (electrons); #eta; dE/dx",
265                                  50, -1.0, 1.0, 
266                                    95-fDeDxMIPMax, fDeDxMIPMax, 95);
267     hDeDxVsEtaElectrons->Sumw2();
268     hDeDxVsEtaElectrons->SetDirectory(0);
269
270     hNclVsEta = new TH2D("hNclVsEta", "Ncl vs #eta (pion MIPs); #eta; Ncl",
271                          50, -1.0, 1.0, 
272                          18, 69.5, 159.5); 
273     hNclVsEta->Sumw2();
274     hNclVsEta->SetDirectory(0);
275     hNclVsEtaElectrons = new TH2D("hNclVsEtaElectrons", "Ncl vs #eta (electrons); #eta; Ncl",
276                                   50, -1.0, 1.0, 
277                                   18, 69.5, 159.5); 
278     hNclVsEtaElectrons->Sumw2();
279     hNclVsEtaElectrons->SetDirectory(0);
280     
281     fInit = 2;
282
283     break;
284
285   case 3:
286
287   //
288   // step 3 histograms 
289   //
290     hSelection3 = new TH2D("hSelection3", "dE/dx vs P (pion MIPs) step 3; P [GeV/c]; dE/dx",
291                            100, 0.2, 0.8, 
292                            100, 0, 100);
293     hSelection3->SetDirectory(0);
294
295     hSelectionElectrons3 = new TH2D("hSelectionElectrons3", "dE/dx vs P (electrons) step 3; P [GeV/c]; dE/dx",
296                                    100, 0.2, 0.8, 
297                                    100, 0, 100);
298     hSelectionElectrons3->SetDirectory(0);
299     
300
301     hMeanP = new TProfile("hMeanP", "mean p; p [GeV/c]; mean p",
302                           nPtBins, ptBins);
303     hMeanP->SetDirectory(0);
304     
305     hDeDxVsP = new TH2D("hDeDxVsP", "dE/dx vs P; p [GeV/c]; dE/dx",
306                         nPtBins, ptBins, 55, 40, 95);
307     hDeDxVsP->Sumw2();
308     hDeDxVsP->SetDirectory(0);
309     // dE/dx vs p 
310     if(fIsMc) {
311       
312       hDeDxVsPPiMc = new TH2D("hDeDxVsPPiMc", "dE/dx vs P; p [GeV/c]; dE/dx",
313                               nPtBins, ptBins, 55, 40, 95);
314       hDeDxVsPPiMc->Sumw2();
315       hDeDxVsPPiMc->SetDirectory(0);
316       
317       hDeDxVsPKMc = new TH2D("hDeDxVsPKMc", "dE/dx vs P; p [GeV/c]; dE/dx",
318                              nPtBins, ptBins, 55, 40, 95);
319       hDeDxVsPKMc->Sumw2();
320       hDeDxVsPKMc->SetDirectory(0);
321       
322       hDeDxVsPPMc = new TH2D("hDeDxVsPPMc", "dE/dx vs P; p [GeV/c]; dE/dx",
323                              nPtBins, ptBins, 55, 40, 95);
324       hDeDxVsPPMc->Sumw2();
325       hDeDxVsPPMc->SetDirectory(0);
326       
327       hDeDxVsPEMc = new TH2D("hDeDxVsPEMc", "dE/dx vs P; p [GeV/c]; dE/dx",
328                              nPtBins, ptBins, 55, 40, 95);
329       hDeDxVsPEMc->Sumw2();
330       hDeDxVsPEMc->SetDirectory(0);
331       
332       hDeDxVsPMuMc = new TH2D("hDeDxVsPMuMc", "dE/dx vs P; p [GeV/c]; dE/dx",
333                               nPtBins, ptBins, 55, 40, 95);
334       hDeDxVsPMuMc->Sumw2();
335       hDeDxVsPMuMc->SetDirectory(0);
336     }
337     
338     //
339     // step 3 histograms reated to the eta calibration
340     //
341     
342     hDeDxVsEtaCal = new TH2D("hDeDxVsEtaCal", "dE/dx vs #eta calibrated (pion MIPs); #eta; dE/dx",
343                              50, -1.0, 1.0, 
344                              fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
345     hDeDxVsEtaCal->Sumw2();
346     hDeDxVsEtaCal->SetDirectory(0);
347     
348     hDeDxVsEtaCalElectrons = new TH2D("hDeDxVsEtaCalElectrons", "dE/dx vs #eta calibrated (electrons); #eta; dE/dx",
349                                       50, -1.0, 1.0, 
350                                       95-fDeDxMIPMax, fDeDxMIPMax, 95);
351     hDeDxVsEtaCalElectrons->Sumw2();
352     hDeDxVsEtaCalElectrons->SetDirectory(0);
353   
354     //
355     // step 3 histograms related to resolution
356     //
357     hMeanEta = new TProfile("hMeanEta", "<|#eta|> in |#eta| intervals (pion MIPs); |#eta|; <|#eta|>",
358                             4, 0, 0.8);
359     hMeanEta->SetDirectory(0);
360
361     hMeanEtaElectrons = new TProfile("hMeanEtaElectrons", "<|#eta|> in |#eta| intervals (electrons); |#eta|; <|#eta|>",
362                                      4, 0, 0.8);
363     hMeanEtaElectrons->SetDirectory(0);
364     
365     hDeDx = new TH1D("hDeDx", "dE/dx (pion MIPs); dE/dx; Counts",
366                      fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
367     hDeDx->Sumw2();
368     hDeDx->SetDirectory(0);
369     hDeDx1 = new TH1D("hDeDx1", "dE/dx |#eta|<0.2 (pion MIPs); dE/dx; Counts",
370                       fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
371     hDeDx1->Sumw2();
372     hDeDx1->SetDirectory(0);
373     hDeDx2 = new TH1D("hDeDx2", "dE/dx 0.2<|#eta|<0.4 (pion MIPs); dE/dx; Counts",
374                       fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
375     hDeDx2->Sumw2();
376     hDeDx2->SetDirectory(0);
377     hDeDx3 = new TH1D("hDeDx3", "dE/dx 0.4<|#eta|<0.6 (pion MIPs); dE/dx; Counts",
378                       fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
379     hDeDx3->Sumw2();
380     hDeDx3->SetDirectory(0);
381     hDeDx4 = new TH1D("hDeDx4", "dE/dx 0.6<|#eta|<0.8 (pion MIPs); dE/dx; Counts",
382                       fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
383     hDeDx4->Sumw2();
384     hDeDx4->SetDirectory(0);
385     
386     hDeDxElectrons = new TH1D("hDeDxElectrons", "dE/dx (electrons); dE/dx; Counts",
387                               95-fDeDxMIPMax, fDeDxMIPMax, 95);
388     hDeDxElectrons->Sumw2();
389     hDeDxElectrons->SetDirectory(0);
390     hDeDxElectrons1 = new TH1D("hDeDxElectrons1", "dE/dx |#eta|<0.2 (electrons); dE/dx; Counts",
391                                95-fDeDxMIPMax, fDeDxMIPMax, 95);
392     hDeDxElectrons1->Sumw2();
393     hDeDxElectrons1->SetDirectory(0);
394     hDeDxElectrons2 = new TH1D("hDeDxElectrons2", "dE/dx 0.2<|#eta|<0.4 (electrons); dE/dx; Counts",
395                                95-fDeDxMIPMax, fDeDxMIPMax, 95);
396     hDeDxElectrons2->Sumw2();
397     hDeDxElectrons2->SetDirectory(0);
398     hDeDxElectrons3 = new TH1D("hDeDxElectrons3", "dE/dx 0.4<|#eta|<0.6 (electrons); dE/dx; Counts",
399                                95-fDeDxMIPMax, fDeDxMIPMax, 95);
400     hDeDxElectrons3->Sumw2();
401     hDeDxElectrons3->SetDirectory(0);
402     hDeDxElectrons4 = new TH1D("hDeDxElectrons4", "dE/dx 0.6<|#eta|<0.8 (electrons); dE/dx; Counts",
403                                95-fDeDxMIPMax, fDeDxMIPMax, 95);
404     hDeDxElectrons4->Sumw2();
405     hDeDxElectrons4->SetDirectory(0);
406     
407     //
408     // step 3 histograms for resolution
409     //
410     
411     
412     
413     hDeDxVsNcl = new TH2D("hDeDxVsNcl", "dE/dx vs Ncl (pion MIPs); Ncl; dE/dx",
414                           18, 69.5, 159.5, 
415                           fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
416     hDeDxVsNcl->Sumw2();
417     hDeDxVsNcl->SetDirectory(0);
418     hDeDxVsNcl1 = new TH2D("hDeDxVsNcl1", "dE/dx vs Ncl |#eta|<0.2 (pion MIPs); Ncl; dE/dx",
419                            18, 69.5, 159.5, 
420                            fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
421     hDeDxVsNcl1->Sumw2();
422     hDeDxVsNcl1->SetDirectory(0);
423     hDeDxVsNcl2 = new TH2D("hDeDxVsNcl2", "dE/dx vs Ncl 0.2<|#eta|<0.4 (pion MIPs); Ncl; dE/dx",
424                            18, 69.5, 159.5, 
425                            fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
426     hDeDxVsNcl2->Sumw2();
427     hDeDxVsNcl2->SetDirectory(0);
428     hDeDxVsNcl3 = new TH2D("hDeDxVsNcl3", "dE/dx vs Ncl 0.4<|#eta|<0.6 (pion MIPs); Ncl; dE/dx",
429                            18, 69.5, 159.5, 
430                            fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
431     hDeDxVsNcl3->Sumw2();
432     hDeDxVsNcl3->SetDirectory(0);
433     hDeDxVsNcl4 = new TH2D("hDeDxVsNcl4", "dE/dx vs Ncl 0.6<|#eta|<0.8 (pion MIPs); Ncl; dE/dx",
434                            18, 69.5, 159.5, 
435                            fDeDxMIPMax-fDeDxMIPMin, fDeDxMIPMin, fDeDxMIPMax);
436     hDeDxVsNcl4->Sumw2();
437     hDeDxVsNcl4->SetDirectory(0);
438     
439     hDeDxVsNclElectrons = new TH2D("hDeDxVsNclElectrons", "dE/dx vs Ncl (electrons); Ncl; dE/dx",
440                                    18, 69.5, 159.5, 
441                                    95-fDeDxMIPMax, fDeDxMIPMax, 95);
442     hDeDxVsNclElectrons->Sumw2();
443     hDeDxVsNclElectrons->SetDirectory(0);
444     hDeDxVsNclElectrons1 = new TH2D("hDeDxVsNclElectrons1", "dE/dx vs Ncl |#eta|<0.2 (electrons); Ncl; dE/dx",
445                                     18, 69.5, 159.5, 
446                                     95-fDeDxMIPMax, fDeDxMIPMax, 95);
447     hDeDxVsNclElectrons1->Sumw2();
448     hDeDxVsNclElectrons1->SetDirectory(0);
449     hDeDxVsNclElectrons2 = new TH2D("hDeDxVsNclElectrons2", "dE/dx vs Ncl 0.2<|#eta|<0.4 (electrons); Ncl; dE/dx",
450                          18, 69.5, 159.5, 
451                                     95-fDeDxMIPMax, fDeDxMIPMax, 95);
452     hDeDxVsNclElectrons2->Sumw2();
453     hDeDxVsNclElectrons2->SetDirectory(0);
454     hDeDxVsNclElectrons3 = new TH2D("hDeDxVsNclElectrons3", "dE/dx vs Ncl 0.4<|#eta|<0.6 (electrons); Ncl; dE/dx",
455                                     18, 69.5, 159.5, 
456                          95-fDeDxMIPMax, fDeDxMIPMax, 95);
457     hDeDxVsNclElectrons3->Sumw2();
458     hDeDxVsNclElectrons3->SetDirectory(0);
459     hDeDxVsNclElectrons4 = new TH2D("hDeDxVsNclElectrons4", "dE/dx vs Ncl 0.6<|#eta|<0.8 (electrons); Ncl; dE/dx",
460                                     18, 69.5, 159.5, 
461                                     95-fDeDxMIPMax, fDeDxMIPMax, 95);
462     hDeDxVsNclElectrons4->Sumw2();
463     hDeDxVsNclElectrons4->SetDirectory(0);
464
465     fInit = 3;
466
467     break;
468   default:
469     cout << "No init implemented for step: " << step << endl;
470   }
471 }
472
473 //_________________________________________________________
474 void AliHighPtDeDxCalib::FillTrackInfo(Float_t weight)
475 {
476   if(fStep > 2) { // calibate eta dependence
477     
478     if(fTrackEta < 0) 
479       fTrackDeDx *= 50.0 / fDeDxVsEtaNeg->Eval(fTrackEta);
480     else
481       fTrackDeDx *= 50.0 / fDeDxVsEtaPos->Eval(fTrackEta);
482   }
483
484   if(fStep > 1) { // calibate ncl dependence
485     
486     fTrackDeDx *= 50.0 / fDeDxVsNcl->Eval(fTrackNcl);
487   }
488   
489   //
490   // Fill information for MIP pions
491   //
492   if(IsMIP()) {
493     
494     if(fStep == 2) { // calibate eta
495       
496       hSelection2->Fill(fTrackP, fTrackDeDx);
497       hDeDxVsEta->Fill(fTrackEta, fTrackDeDx);
498       hNclVsEta->Fill(fTrackEta, fTrackNcl);
499     } else if (fStep==1) { // calibrate nCl dependence
500       
501       hSelection1->Fill(fTrackP, fTrackDeDx);
502       hDeDxVsNclBefore->Fill(fTrackNcl, fTrackDeDx);
503
504     } else if (fStep==3) { // calibrate <dE/dx>pi and check eta calibration
505       
506       AliHighPtDeDxBase::FillTrackInfo(weight);
507
508       hSelection3->Fill(fTrackP, fTrackDeDx);
509
510       hDeDxVsEtaCal->Fill(fTrackEta, fTrackDeDx);
511       
512       hMeanEta->Fill(TMath::Abs(fTrackEta), TMath::Abs(fTrackEta));
513       hDeDx->Fill(fTrackDeDx);
514       hDeDxVsNcl->Fill(fTrackNcl, fTrackDeDx);
515       if(TMath::Abs(fTrackEta)<0.2) {
516         hDeDx1->Fill(fTrackDeDx);
517         hDeDxVsNcl1->Fill(fTrackNcl, fTrackDeDx);
518       } else if(TMath::Abs(fTrackEta)>=0.2&&TMath::Abs(fTrackEta)<0.4) {
519         hDeDx2->Fill(fTrackDeDx);
520         hDeDxVsNcl2->Fill(fTrackNcl, fTrackDeDx);
521       } else if(TMath::Abs(fTrackEta)>=0.4&&TMath::Abs(fTrackEta)<0.6) {
522         hDeDx3->Fill(fTrackDeDx);
523         hDeDxVsNcl3->Fill(fTrackNcl, fTrackDeDx);
524       }else if(TMath::Abs(fTrackEta)>=0.6&&TMath::Abs(fTrackEta)<0.8) { 
525         hDeDx4->Fill(fTrackDeDx);
526         hDeDxVsNcl4->Fill(fTrackNcl, fTrackDeDx);
527       }
528     }
529   }
530
531   //
532   // Fill information for electrons with same p as MIP pions
533   //
534   // This is done to validate the calibration/assumptions about sigma for a
535   // group of tracks with same topology but much higher dE/dx (plateau)
536   //
537   if(fStep > 1 && IsElectron()) {
538
539     if(fStep == 2) { // calibate eta
540       
541       hSelectionElectrons2->Fill(fTrackP, fTrackDeDx);
542       hDeDxVsEtaElectrons->Fill(fTrackEta, fTrackDeDx);
543       hNclVsEtaElectrons->Fill(fTrackEta, fTrackNcl);
544     } else if (fStep==3) { // calibrate <dE/dx>pi and check eta calibration
545       
546       hSelectionElectrons3->Fill(fTrackP, fTrackDeDx);
547
548       hDeDxVsEtaCalElectrons->Fill(fTrackEta, fTrackDeDx);
549       
550       hMeanEtaElectrons->Fill(TMath::Abs(fTrackEta), TMath::Abs(fTrackEta));
551       hDeDxElectrons->Fill(fTrackDeDx);
552       hDeDxVsNclElectrons->Fill(fTrackNcl, fTrackDeDx);
553       if(TMath::Abs(fTrackEta)<0.2) {
554         hDeDxElectrons1->Fill(fTrackDeDx);
555         hDeDxVsNclElectrons1->Fill(fTrackNcl, fTrackDeDx);
556       } else if(TMath::Abs(fTrackEta)>=0.2&&TMath::Abs(fTrackEta)<0.4) {
557         hDeDxElectrons2->Fill(fTrackDeDx);
558         hDeDxVsNclElectrons2->Fill(fTrackNcl, fTrackDeDx);
559       } else if(TMath::Abs(fTrackEta)>=0.4&&TMath::Abs(fTrackEta)<0.6) {
560         hDeDxElectrons3->Fill(fTrackDeDx);
561         hDeDxVsNclElectrons3->Fill(fTrackNcl, fTrackDeDx);
562       }else if(TMath::Abs(fTrackEta)>=0.6&&TMath::Abs(fTrackEta)<0.8) { 
563         hDeDxElectrons4->Fill(fTrackDeDx);
564         hDeDxVsNclElectrons4->Fill(fTrackNcl, fTrackDeDx);
565       }
566     }
567   }
568
569   //
570   // Fill information for high pT tracks (dE/dx vs p)
571   //
572   if(fStep==3) { // Fill dE/dx vs 
573     
574     hDeDxVsP->Fill(fTrackP, fTrackDeDx, weight);
575     hMeanP->Fill(fTrackP, fTrackP);    
576
577     if(fIsMc) {
578       
579       switch (fTrackPidMc) {
580         
581       case 1: // pion
582           hDeDxVsPPiMc->Fill(fTrackP, fTrackDeDx, weight);
583         break;
584       case 2: // kaon
585           hDeDxVsPKMc ->Fill(fTrackP, fTrackDeDx, weight);
586         break;
587       case 3: // proton
588           hDeDxVsPPMc ->Fill(fTrackP, fTrackDeDx, weight);
589         break;
590       case 4: // electron
591           hDeDxVsPEMc ->Fill(fTrackP, fTrackDeDx, weight);
592         break;
593       case 5: // muon
594           hDeDxVsPMuMc ->Fill(fTrackP, fTrackDeDx, weight);
595         break;
596       default:
597         break;
598       }
599     }
600   }
601 }
602   
603
604 //___________________________________________________________________________
605 Bool_t AliHighPtDeDxCalib::IsMIP()
606 {
607   if(fTrackP > fPMIPMin && fTrackP<fPMIPMax &&
608      fTrackDeDx > fDeDxMIPMin && fTrackDeDx < fDeDxMIPMax)
609     return kTRUE;
610   
611   return kFALSE;
612 }
613
614 //___________________________________________________________________________
615 Bool_t AliHighPtDeDxCalib::IsElectron()
616 {
617   if(fTrackP > fPMIPMin && fTrackP<fPMIPMax &&
618      fTrackDeDx > (fDeDxMIPMax+5) && TMath::Abs(fTrackBeta-1)<fDeltaBeta)
619     return kTRUE;
620   
621   return kFALSE;
622 }
623
624
625 //___________________________________________________________________________
626 void AliHighPtDeDxCalib::PerformEtaCal()
627 {
628   TProfile* hDeDxVsEtaProf = hDeDxVsEta->ProfileX();
629   hDeDxVsEtaProf->SetMarkerStyle(29);
630   hDeDxVsEtaProf->Fit(fDeDxVsEtaNeg, "0", "", -1, 0);
631   hDeDxVsEtaProf->Fit(fDeDxVsEtaPos, "0", "",  0, 1);
632   delete hDeDxVsEtaProf;
633 }
634
635 //___________________________________________________________________________
636 void AliHighPtDeDxCalib::PerformNclCal()
637 {
638   TProfile* hDeDxVsNclProf = hDeDxVsNclBefore->ProfileX();
639   hDeDxVsNclProf->SetMarkerStyle(29);
640   hDeDxVsNclProf->Fit(fDeDxVsNcl, "0", "", 69.5, 159.5);
641   delete hDeDxVsNclProf;
642 }
643
644 //___________________________________________________________________________
645 TCanvas* AliHighPtDeDxCalib::DrawNclCal()
646 {
647   TCanvas* cNcl = new TCanvas("cNcl", "dE/dx vs Ncl for MIP pions", 600, 400);
648   cNcl->Clear();
649   cNcl->cd();
650   hDeDxVsNclBefore->DrawCopy("COL");
651   TProfile* hDeDxVsNclProf = hDeDxVsNclBefore->ProfileX();
652   hDeDxVsNclProf->SetMarkerStyle(29);
653   hDeDxVsNclProf->DrawCopy("SAME");
654   fDeDxVsNcl->Draw("SAME");
655   delete hDeDxVsNclProf;
656   return cNcl;
657 }
658
659 //___________________________________________________________________________
660 TCanvas* AliHighPtDeDxCalib::DrawEta(Bool_t forMIP)
661 {
662   // Draw dE/dx vs eta for step 1 (calibration)
663   if(forMIP) {
664
665     TCanvas* cEta = new TCanvas("cEta", "dE/dx vs Eta for MIP pions", 600, 400);
666     cEta->Clear();
667     cEta->cd();
668     hDeDxVsEta->DrawCopy("COL");
669
670     TProfile* hDeDxVsEtaProf = hDeDxVsEta->ProfileX();
671     hDeDxVsEtaProf->SetMarkerStyle(29);
672     hDeDxVsEtaProf->DrawCopy("SAME");
673     fDeDxVsEtaNeg->DrawCopy("SAME");
674     fDeDxVsEtaPos->DrawCopy("SAME");
675
676     return cEta;
677   }
678
679   TCanvas* cEta = new TCanvas("cEtaElectrons", "dE/dx vs Eta for electrons", 600, 400);
680   cEta->Clear();
681   cEta->cd();
682   hDeDxVsEtaElectrons->DrawCopy("COL");
683
684   TProfile* hDeDxVsEtaProfElectrons = hDeDxVsEtaElectrons->ProfileX();
685   hDeDxVsEtaProfElectrons->SetMarkerStyle(29);
686   hDeDxVsEtaProfElectrons->DrawCopy("SAME");
687
688   return cEta;
689 }     
690
691 //___________________________________________________________________________
692 TCanvas* AliHighPtDeDxCalib::DrawEtaCalibrated(Bool_t forMIP)
693 {
694   // Draw dE/dx vs eta for step 2 (after calibration)
695   if(forMIP) {
696
697     TCanvas* cEtaCal = new TCanvas("cEtaCal", "dE/dx vs Eta for MIP pions (calibrated)", 600, 400);
698     cEtaCal->Clear();
699     cEtaCal->cd();
700     hDeDxVsEtaCal->DrawCopy("COL");
701
702     TProfile* hDeDxVsEtaCalProf = hDeDxVsEtaCal->ProfileX();
703     hDeDxVsEtaCalProf->SetMarkerStyle(29);
704     hDeDxVsEtaCalProf->DrawCopy("SAME");
705     
706     return cEtaCal;
707   }
708
709   TCanvas* cEtaCalElectrons = new TCanvas("cEtaCalElectrons", "dE/dx vs Eta (calibrated electrons)", 600, 400);
710   cEtaCalElectrons->Clear();
711   cEtaCalElectrons->cd();
712   hDeDxVsEtaCalElectrons->DrawCopy("COL");
713   TProfile* hDeDxVsEtaCalElectronsProf = 
714     hDeDxVsEtaCalElectrons->ProfileX();
715   hDeDxVsEtaCalElectronsProf->SetMarkerStyle(29);
716   hDeDxVsEtaCalElectronsProf->DrawCopy("SAME");
717   
718   return cEtaCalElectrons;
719 }     
720
721 //___________________________________________________________________________
722 TCanvas* AliHighPtDeDxCalib::DrawSelectionHistograms(Int_t step)
723 {
724   TCanvas* cSelection = new TCanvas("cSelection", "dE/dx vs p selection", 800, 800);
725   cSelection->Clear();
726   cSelection->Divide(2,2);
727
728   cSelection->cd(1);
729   if(step==1)
730     hSelection1->DrawCopy("COL");
731   else if(step==2)
732     hSelection2->DrawCopy("COL");
733   else if(step==3)
734     hSelection3->DrawCopy("COL");
735
736   cSelection->cd(2);
737   if(step==2)
738     hSelectionElectrons2->DrawCopy("COL");
739   else if(step==3)
740     hSelectionElectrons3->DrawCopy("COL");
741   
742   cSelection->cd(3);
743   hNclVsEta->DrawCopy("COL");
744
745   cSelection->cd(4);
746   hNclVsEtaElectrons->DrawCopy("COL");
747   
748   return cSelection;
749 }
750
751 //___________________________________________________________________________
752 TH1D* AliHighPtDeDxCalib::GetHistDeDx(Bool_t forMIP, Int_t etaBin)
753 {
754   switch (etaBin) {
755         
756   case 0:
757     if(forMIP)
758       return hDeDx;
759     else
760       return hDeDxElectrons;
761     break;
762   case 1:
763     if(forMIP)
764       return hDeDx1;
765     else
766       return hDeDxElectrons1;
767     break;
768   case 2:
769     if(forMIP)
770       return hDeDx2;
771     else
772       return hDeDxElectrons2;
773     break;
774   case 3:
775     if(forMIP)
776       return hDeDx3;
777     else
778       return hDeDxElectrons3;
779     break;
780   case 4:
781     if(forMIP)
782       return hDeDx4;
783     else
784       return hDeDxElectrons4;
785     break;
786   default:
787     cout << "Eta bin: " << etaBin << " not found" << endl;
788     break;
789   }
790   return 0;
791 }
792
793 //___________________________________________________________________________
794 TH2D* AliHighPtDeDxCalib::GetHistDeDxVsNcl(Bool_t forMIP, Int_t etaBin)
795 {
796   switch (etaBin) {
797         
798   case 0:
799     if(forMIP)
800       return hDeDxVsNcl;
801     else
802       return hDeDxVsNclElectrons;
803     break;
804   case 1:
805     if(forMIP)
806       return hDeDxVsNcl1;
807     else
808       return hDeDxVsNclElectrons1;
809     break;
810   case 2:
811     if(forMIP)
812       return hDeDxVsNcl2;
813     else
814       return hDeDxVsNclElectrons2;
815     break;
816   case 3:
817     if(forMIP)
818       return hDeDxVsNcl3;
819     else
820       return hDeDxVsNclElectrons3;
821     break;
822   case 4:
823     if(forMIP)
824       return hDeDxVsNcl4;
825     else
826       return hDeDxVsNclElectrons4;
827     break;
828   default:
829     cout << "Eta bin: " << etaBin << " not found" << endl;
830     break;
831   }
832   return 0;
833 }
834
835 //___________________________________________________________________________
836 TH2D* AliHighPtDeDxCalib::GetHistDeDxVsP(Int_t pid)
837 {
838   switch (pid) {
839         
840   case 0:
841     return hDeDxVsP;
842     break;
843   case 1:
844     return hDeDxVsPPiMc;
845     break;
846   case 2:
847     return hDeDxVsPKMc;
848     break;
849   case 3:
850     return hDeDxVsPPMc;
851     break;
852   case 4:
853     return hDeDxVsPEMc;
854     break;
855   case 5:
856     return hDeDxVsPMuMc;
857     break;
858   default:
859     cout << "PID: " << pid << " not found" << endl;
860     break;
861   }
862   return 0;
863 }
864
865 /*
866
867       }
868
869   }
870
871 */