New book-keeping class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskPWG4PidDetEx.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 #include <TROOT.h>
17 #include <TSystem.h>
18 #include <TChain.h>
19 #include <TTree.h>
20 #include <TFile.h>
21 #include <TList.h>
22 #include <TStyle.h>
23 #include <TCanvas.h>
24 #include <TMath.h>
25 #include <TH1.h>
26 #include <TH2.h>
27 #include <TF1.h>
28 #include <TProfile.h>
29 #include <TLegend.h>
30 #include <TDatabasePDG.h>
31 #include <TParticle.h>
32
33 #include "AliAnalysisManager.h"
34 #include "AliAnalysisFilter.h"
35 #include "AliESDInputHandler.h"
36 #include "AliESDEvent.h"
37 #include "AliESDVertex.h"
38 #include "AliMCEventHandler.h"
39 #include "AliMCEvent.h"
40 #include "AliStack.h"
41 #include "AliAODInputHandler.h"
42 #include "AliAODEvent.h"
43 #include "AliAODVertex.h"
44 #include "AliAODMCParticle.h"
45 #include "AliLog.h"
46
47 #include "AliAnalysisTaskPWG4PidDetEx.h"
48
49 // STL includes
50 #include <iostream>
51 using namespace std;
52
53 //
54 // Analysis class example for PID using detector signals.
55 // Works on ESD and AOD; has a flag for MC analysis 
56 //
57 //    Alexandru.Dobrin@hep.lu.se
58 //    Peter.Christiansen@hep.lu.se
59 // 
60
61 ClassImp(AliAnalysisTaskPWG4PidDetEx)
62
63 const Double_t AliAnalysisTaskPWG4PidDetEx::fgkCtau = 370;                //  distance for kaon decay
64 const Double_t AliAnalysisTaskPWG4PidDetEx::fgkPionMass = 1.39570000000000000e-01;
65 const Double_t AliAnalysisTaskPWG4PidDetEx::fgkKaonMass = 4.93599999999999983e-01; 
66 const Double_t AliAnalysisTaskPWG4PidDetEx::fgkProtonMass = 9.38270000000000048e-01;
67 const Double_t AliAnalysisTaskPWG4PidDetEx::fgkC = 2.99792458e-2;
68
69
70
71
72
73 //_____________________________________________________________________________
74 AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx():
75   AliAnalysisTaskSE(),
76   fESD(0),
77   fAOD(0), 
78   fListOfHists(0), 
79   fEtaCut(0.9), 
80   fPtCut(0.5), 
81   fXbins(20), 
82   fXmin(0.),  
83   fTOFCutP(2.), 
84   fTrackFilter(0x0),
85   fAnalysisMC(kTRUE),
86   fAnalysisType("ESD"),
87   fTriggerMode(kMB2),
88   fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0),
89   fdNdPt(0), fMassAll(0),
90   fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0),
91   fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0),
92   fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0),
93   fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0) 
94 {
95   // Default constructor (should not be used)
96 }
97
98 //______________________________________________________________________________
99 AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx(const char *name):
100   AliAnalysisTaskSE(name),
101   fESD(0),
102   fAOD(0), 
103   fListOfHists(0), 
104   fEtaCut(0.9), 
105   fPtCut(0.5), 
106   fXbins(20), 
107   fXmin(0.), 
108   fTOFCutP(2.), 
109   fTrackFilter(0x0),
110   fAnalysisMC(kTRUE),
111   fAnalysisType("ESD"),
112   fTriggerMode(kMB2),
113   fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0),
114   fdNdPt(0), fMassAll(0),
115   fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0),
116   fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0),
117   fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0),
118   fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0) 
119 {
120   // Output slot #0 writes into a TList
121   DefineOutput(1, TList::Class());
122
123 }
124
125 //_____________________________________________________________________________
126 AliAnalysisTaskPWG4PidDetEx::~AliAnalysisTaskPWG4PidDetEx()
127 {
128   // Destructor
129 }
130
131 // //______________________________________________________________________________
132 // void AliAnalysisTaskPWG4PidDetEx::ConnectInputData(Option_t *)
133 // {
134 //   // Connect AOD here
135 //   // Called once
136 //   if (fDebug > 1)  AliInfo("ConnectInputData() \n");
137   
138 //   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
139 //   if (!tree) {
140 //     Printf("ERROR: Could not read chain from input slot 0");
141 //   } 
142 //   else {
143 //     if(fAnalysisType == "ESD") {
144 //       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
145 //       if (!esdH) {
146 //      Printf("ERROR: Could not get ESDInputHandler");
147 //       } else
148 //      fESD = esdH->GetEvent();
149 //     }
150 //     else if(fAnalysisType == "AOD") {
151 //       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
152 //       if (!aodH) {
153 //      Printf("ERROR: Could not get AODInputHandler");
154 //       } else
155 //      fAOD = aodH->GetEvent();
156 //     }
157 //   }
158
159 // }
160
161 //______________________________________________________________________________
162 void AliAnalysisTaskPWG4PidDetEx::UserCreateOutputObjects()
163
164   OpenFile(1);
165   fListOfHists = new TList();
166
167   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1);
168   fListOfHists->Add(fEvents);
169
170   fEffTot = new TH1F ("fEffTot","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
171   fEffTot->Sumw2();
172   fListOfHists->Add(fEffTot);
173
174   fEffPID = new TH1F ("fEffPID","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
175   fEffPID->Sumw2();
176   fListOfHists->Add(fEffPID);
177
178   fAccP = new TH1F ("fAccP","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
179   fAccP->Sumw2();
180   fListOfHists->Add(fAccP);
181
182   fAccPt = new TH1F ("fAccPt","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
183   fAccPt->Sumw2();
184   fListOfHists->Add(fAccPt);
185
186   fKaonDecayCorr = new TProfile("fKaonDecayCorr","Kaon decay correction vs p_{T}; p_{T} [GeV/c]; Kaon decay correction", fXbins, fXmin, fTOFCutP);
187   fListOfHists->Add(fKaonDecayCorr);
188
189   fdNdPt = new TH1F("fdNdPt","p_{T} distribution; p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100 , 0, 10);
190   fdNdPt->Sumw2();
191   fListOfHists->Add(fdNdPt);
192
193   fMassAll = new TH2F("fMassAll","Mass^{2} vs momentum (ToF); p [GeV/c]; M^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
194   fListOfHists->Add(fMassAll);
195
196   //Pion
197   fdNdPtPion = new TH1F("fdNdPtPion","p_{T} distribution (Pions); p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
198   fdNdPtPion->Sumw2();
199   fListOfHists->Add(fdNdPtPion);
200
201   fMassPion = new TH2F("fMassPion","Mass squared for pions after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
202   fListOfHists->Add(fMassPion);
203
204   fdEdxTPCPion = new TH2F("fdEdxTPCPion","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
205   fListOfHists->Add(fdEdxTPCPion);
206
207   fbgTPCPion = new TH2F("fbgTPCPion", "Energy loss vs #beta#gamma (Pions);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
208   fListOfHists->Add(fbgTPCPion);
209
210   //Kaon
211   fdNdPtKaon = new TH1F("fdNdPtKaon","p_{T} distribution (Kaons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
212   fdNdPtKaon->Sumw2();
213   fListOfHists->Add(fdNdPtKaon);
214
215   fMassKaon = new TH2F("fMassKaon","Mass squared for kaons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
216   fListOfHists->Add(fMassKaon);
217
218   fdEdxTPCKaon = new TH2F("fdEdxTPCKaon","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
219   fListOfHists->Add(fdEdxTPCKaon);
220
221   fbgTPCKaon = new TH2F("fbgTPCKaon", "Energy loss vs #beta#gamma (Kaons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
222   fListOfHists->Add(fbgTPCKaon);
223
224   //Proton
225   fdNdPtProton = new TH1F("fdNdPtProton","p_{T} distribution (Protons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
226   fdNdPtProton->Sumw2();
227   fListOfHists->Add(fdNdPtProton); 
228   
229   fMassProton = new TH2F("fMassProton","Mass squared for protons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
230   fListOfHists->Add(fMassProton);
231
232   fdEdxTPCProton = new TH2F("fdEdxTPCProton","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
233   fListOfHists->Add(fdEdxTPCProton);
234
235   fbgTPCProton = new TH2F("fbgTPCProton", "Energy loss vs #beta#gamma (Protons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
236   fListOfHists->Add(fbgTPCProton);
237
238
239   //MC
240   if(fAnalysisMC){
241     fdNdPtMC = new TH1F("fdNdPtMC","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100, 0, 10);
242     fdNdPtMC->SetLineColor(2);
243     fdNdPtMC->Sumw2();
244     fListOfHists->Add(fdNdPtMC);
245
246     fdNdPtMCPion = new TH1F("fdNdPtMCPion","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
247     fdNdPtMCPion->SetLineColor(2);
248     fdNdPtMCPion->Sumw2();
249     fListOfHists->Add(fdNdPtMCPion);
250   
251     fdNdPtMCKaon = new TH1F("fdNdPtMCKaon","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
252     fdNdPtMCKaon->SetLineColor(2);
253     fdNdPtMCKaon->Sumw2();
254     fListOfHists->Add(fdNdPtMCKaon);
255
256     fdNdPtMCProton = new TH1F("fdNdPtMCProton","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
257     fdNdPtMCProton->SetLineColor(2);
258     fdNdPtMCProton->Sumw2();
259     fListOfHists->Add(fdNdPtMCProton);
260   }
261
262 }
263
264 //______________________________________________________________________________
265 void AliAnalysisTaskPWG4PidDetEx::UserExec(Option_t *) 
266 {
267
268   // Create histograms
269   // Called once
270   if (fAnalysisType == "AOD") {
271     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
272     if(!fAOD){
273       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
274       return;
275     }
276     else{
277       //  assume that the AOD is in the general output...
278       //  fAOD  = AODEvent();
279       //  if(!fAOD){
280       //        Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
281     }
282   }
283   else if  (fAnalysisType == "ESD"){
284     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
285     if(!fESD){
286       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
287       this->Dump();
288       return;
289     }
290   }
291
292   // Main loop
293   // Called for each event
294   if (fDebug > 1) AliInfo("Exec() \n" );
295
296   if(fAnalysisType == "ESD") {
297     if (!fESD) {
298       Printf("ERROR: fESD not available");
299       return;
300     }
301     if(IsEventTriggered(fESD, fTriggerMode)) {
302       Printf("PWG4 PID ESD analysis");
303       AnalyzeESD(fESD);
304     }//triggered event
305   }//ESD analysis  
306               
307   else if(fAnalysisType == "AOD") {
308     if (!fAOD) {
309       Printf("ERROR: fAOD not available");
310       return;
311     } 
312     if(IsEventTriggered(fAOD, fTriggerMode)) {
313       Printf("PWG4 PID AOD analysis");
314       AnalyzeAOD(fAOD);
315     }//triggered event
316   }//AOD analysis  
317
318   // Post output data.
319   PostData(1, fListOfHists);
320 }
321
322 //________________________________________________________________________
323 void AliAnalysisTaskPWG4PidDetEx::AnalyzeESD(AliESDEvent* esd)
324 {
325   // Get vertex 
326   const AliESDVertex* primaryVertex = esd->GetPrimaryVertex();   
327   const Double_t vertexResZ = primaryVertex->GetZRes();
328
329   // Select only events with a reconstructed vertex
330   if(vertexResZ<5.0) {
331     const Int_t nESDTracks = esd->GetNumberOfTracks();
332     for(Int_t i = 0; i < nESDTracks; i++) {
333       AliESDtrack* trackESD = esd->GetTrack(i);
334
335       //Cuts
336       UInt_t selectInfo = 0;
337       if (fTrackFilter) {
338         selectInfo = fTrackFilter->IsSelected(trackESD);
339         if (!selectInfo) continue;
340       }
341
342       if ((trackESD->Pt() < fPtCut) || (TMath::Abs(trackESD->Eta()) > fEtaCut ))
343         continue;
344
345       //Corrections
346       fEffTot->Fill(trackESD->Pt());
347       if (trackESD->GetTOFsignal() !=0)
348         fEffPID->Fill(trackESD->Pt());
349
350       if (trackESD->P() < fTOFCutP) fAccP->Fill(trackESD->Pt());
351       if (trackESD->Pt() < fTOFCutP) fAccPt->Fill(trackESD->Pt());
352
353       //Analysis
354       fdNdPt->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
355
356       //TOF 
357       if ((trackESD->GetIntegratedLength() == 0) || (trackESD->GetTOFsignal() == 0))
358         continue;
359
360       fMassAll->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
361
362       if ((MassSquared(trackESD) < 0.15) && (MassSquared(trackESD) > -0.15) && (trackESD->P() < fTOFCutP)){
363         fdNdPtPion->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
364         fMassPion->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
365         fdEdxTPCPion->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal());
366         fbgTPCPion->Fill(trackESD->P()/fgkPionMass,trackESD->GetTPCsignal());
367       }
368
369       if ((MassSquared(trackESD) > 0.15) && (MassSquared(trackESD) < 0.45) && (trackESD->P() < fTOFCutP)){
370         fdNdPtKaon->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
371         fMassKaon->Fill(trackESD->Charge()*trackESD->P(),  MassSquared(trackESD));
372         fdEdxTPCKaon->Fill(trackESD->Charge()*trackESD->P(), trackESD->GetTPCsignal());
373         fbgTPCKaon->Fill(trackESD->P()/fgkKaonMass, trackESD->GetTPCsignal());
374         //Kaon decay correction
375         fKaonDecayCorr->Fill(trackESD->Pt(), TMath::Exp(KaonDecay(trackESD)));
376       }
377
378       if ((MassSquared(trackESD) > 0.75) && (MassSquared(trackESD) < 1.05) && (trackESD->P() < fTOFCutP)){
379         fdNdPtProton->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
380         fMassProton->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
381         fdEdxTPCProton->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal());
382         fbgTPCProton->Fill(trackESD->P()/fgkProtonMass,trackESD->GetTPCsignal());
383       }
384     }//ESD track loop
385
386     //MC
387     if(fAnalysisMC){      
388       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
389       if (!mcHandler) {
390         Printf("ERROR: Could not retrieve MC event handler");
391         return;
392       }
393
394       AliMCEvent* mcEvent = mcHandler->MCEvent();
395       if (!mcEvent) {
396         Printf("ERROR: Could not retrieve MC event");
397         return;
398       }
399
400       AliStack* mcStack = mcEvent->Stack();
401       if (!mcStack) {
402         Printf("ERROR: Could not retrieve MC stack");
403         return;
404       }
405
406       const Int_t nTracksMC = mcStack->GetNtrack();      
407       for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
408         //Cuts
409         if(!(mcStack->IsPhysicalPrimary(iTracks)))
410           continue;
411
412         TParticle* mcTrack = mcStack->Particle(iTracks);
413      
414         Double_t charge = mcTrack->GetPDG()->Charge();
415         if (charge == 0)
416           continue;
417         
418         if ((mcTrack->Pt() < fPtCut) || (TMath::Abs(mcTrack->Eta()) > fEtaCut ))
419           continue;
420
421         //Analysis
422         fdNdPtMC->Fill(mcTrack->Pt(), 1.0/mcTrack->Pt());
423
424         if ((mcTrack->GetPdgCode() == 211) || (mcTrack->GetPdgCode() == -211))
425           fdNdPtMCPion->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
426
427         if ((mcTrack->GetPdgCode() == 321) || (mcTrack->GetPdgCode() == -321))
428           fdNdPtMCKaon->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
429
430         if ((mcTrack->GetPdgCode() == 2212) || (mcTrack->GetPdgCode() == -2212))
431           fdNdPtMCProton->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
432       }//MC track loop 
433     }//if MC
434     fEvents->Fill(0);
435   }//if vertex
436
437 }
438
439 //________________________________________________________________________
440 void AliAnalysisTaskPWG4PidDetEx::AnalyzeAOD(AliAODEvent* aod)
441 {   
442   // Get vertex 
443   AliAODVertex* primaryVertex = aod->GetPrimaryVertex();    
444   const Double_t vertexResZ = TMath::Sqrt(primaryVertex->RotatedCovMatrixZZ());
445
446   // Select only events with a reconstructed vertex
447   if (vertexResZ < 5) {  
448     //AOD
449     Int_t nGoodTracks = aod->GetNumberOfTracks();
450
451     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {     
452       AliAODTrack* trackAOD = aod->GetTrack(iTracks);
453       //Cuts
454       if (!(trackAOD->IsPrimaryCandidate()))
455         continue;
456
457       if ((trackAOD->Pt() < fPtCut) || (TMath::Abs(trackAOD->Eta()) > fEtaCut ))
458         continue;
459
460       //Corrections
461       fEffTot->Fill(trackAOD->Pt());
462       if (trackAOD->GetDetPid()->GetTOFsignal() !=0 )
463         fEffPID->Fill(trackAOD->Pt());
464
465       if (trackAOD->P() < fTOFCutP) fAccP->Fill(trackAOD->Pt());
466       if (trackAOD->Pt() < fTOFCutP) fAccPt->Fill(trackAOD->Pt());
467
468       //Analysis
469       fdNdPt->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
470
471       //TOF
472       if ((IntegratedLength(trackAOD) == 0) || (trackAOD->GetDetPid()->GetTOFsignal() == 0))
473         continue;
474
475       fMassAll->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
476
477       if ((MassSquared(trackAOD) < 0.15) && (MassSquared(trackAOD) > -0.15) && (trackAOD->P() < fTOFCutP)){
478         fdNdPtPion->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
479         fMassPion->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
480         fdEdxTPCPion->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
481         fbgTPCPion->Fill(trackAOD->P()/fgkPionMass,trackAOD->GetDetPid()->GetTPCsignal());
482       }
483
484       if ((MassSquared(trackAOD) > 0.15) && (MassSquared(trackAOD) < 0.45) && (trackAOD->P() < fTOFCutP)){
485         fdNdPtKaon->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
486         fMassKaon->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
487         fdEdxTPCKaon->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
488         fbgTPCKaon->Fill(trackAOD->P()/fgkKaonMass,trackAOD->GetDetPid()->GetTPCsignal());
489         //Kaon decay correction
490         fKaonDecayCorr->Fill(trackAOD->Pt(), TMath::Exp(KaonDecay(trackAOD)));
491       }
492
493       if ((MassSquared(trackAOD) > 0.75) && (MassSquared(trackAOD) < 1.05) && (trackAOD->P() < fTOFCutP)){
494         fdNdPtProton->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
495         fMassProton->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
496         fdEdxTPCProton->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
497         fbgTPCProton->Fill(trackAOD->P()/fgkProtonMass,trackAOD->GetDetPid()->GetTPCsignal());
498       }
499     }//AOD track loop
500
501     //MC
502     if(fAnalysisMC){
503       TClonesArray* farray = (TClonesArray*)aod->FindListObject("mcparticles");
504       Int_t ntrks = farray->GetEntries();
505       for(Int_t i =0 ; i < ntrks; i++){   
506         AliAODMCParticle* trk = (AliAODMCParticle*)farray->At(i);
507         //Cuts
508         if (!(trk->IsPhysicalPrimary()))
509           continue;
510
511         if (trk->Charge() == 0)
512           continue;
513
514         if ((trk->Pt() < fPtCut) || (TMath::Abs(trk->Eta()) > fEtaCut ))
515           continue;
516         
517         //Analysis
518         fdNdPtMC->Fill(trk->Pt(), 1.0/trk->Pt());
519
520         if ((trk->GetPdgCode() == 211) || (trk->GetPdgCode() == -211))
521           fdNdPtMCPion->Fill(trk->Pt(),1.0/trk->Pt());
522
523         if ((trk->GetPdgCode() == 321) || (trk->GetPdgCode() == -321))
524           fdNdPtMCKaon->Fill(trk->Pt(),1.0/trk->Pt());
525
526         if ((trk->GetPdgCode() == 2212) || (trk->GetPdgCode() == -2212))
527           fdNdPtMCProton->Fill(trk->Pt(),1.0/trk->Pt());
528       }//MC track loop 
529     }//if MC 
530     fEvents->Fill(0);   
531   }//if vertex resolution
532
533 }
534
535 //________________________________________________________________________
536 Bool_t AliAnalysisTaskPWG4PidDetEx::IsEventTriggered(AliVEvent* ev, TriggerMode trigger) 
537 {
538   //adapted from PWG2 (AliAnalysisTaskProton)
539
540   ULong64_t triggerMask = ev->GetTriggerMask();
541   
542   // definitions from p-p.cfg
543   ULong64_t spdFO = (1 << 14);
544   ULong64_t v0left = (1 << 11);
545   ULong64_t v0right = (1 << 12);
546
547   switch (trigger) {
548   case kMB1: 
549     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
550       return kTRUE;
551     break;
552   
553   case kMB2: 
554     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
555       return kTRUE;
556     break;
557   
558   case kSPDFASTOR: 
559     if (triggerMask & spdFO)
560       return kTRUE;
561     break;
562   
563   }//switch
564
565   return kFALSE;
566 }
567
568 //_____________________________________________________________________________
569 Double_t AliAnalysisTaskPWG4PidDetEx::IntegratedLength(AliVTrack* track) const
570 {
571   Double_t intTime [5];
572   for (Int_t i = 0; i < 5; i++) intTime[i] = -100.;
573   Double_t timeElectron = 0, intLength = 0;
574
575   AliAODTrack* trackAOD = (AliAODTrack*)track;
576   trackAOD->GetDetPid()->GetIntegratedTimes(intTime);
577   timeElectron = intTime[0];
578   intLength = fgkC*timeElectron;
579
580   return intLength;
581 }
582
583 //_____________________________________________________________________________
584 Double_t AliAnalysisTaskPWG4PidDetEx::MassSquared(AliVTrack* track) const
585 {
586   Double_t beta = -10, mass = -10;
587
588   if(fAnalysisType == "ESD"){
589     AliESDtrack* trackESD = (AliESDtrack*)track;
590     beta = trackESD->GetIntegratedLength()/trackESD->GetTOFsignal()/fgkC;
591     mass = trackESD->P()*trackESD->P()*(1./(beta*beta) - 1.0);
592   }
593
594   if(fAnalysisType == "AOD"){
595     AliAODTrack* trackAOD = (AliAODTrack*)track;
596     beta =IntegratedLength(trackAOD)/trackAOD->GetDetPid()->GetTOFsignal()/fgkC;
597     mass = trackAOD->P()*trackAOD->P()*(1./(beta*beta) - 1.0);
598   }
599   
600   return mass;
601 }
602
603 //_____________________________________________________________________________
604 Double_t AliAnalysisTaskPWG4PidDetEx::KaonDecay(AliVTrack* track) const
605 {
606   Double_t decay = -10;
607
608   if(fAnalysisType == "ESD"){
609     AliESDtrack* trackESD = (AliESDtrack*)track;
610     decay = trackESD->GetIntegratedLength()*fgkKaonMass/fgkCtau/trackESD->P();
611   }
612
613   if (fAnalysisType == "AOD"){
614     AliAODTrack* trackAOD = (AliAODTrack*)track;
615     decay = IntegratedLength(trackAOD)*fgkKaonMass/fgkCtau/trackAOD->P();
616   }
617
618   return decay;
619 }
620
621 //_____________________________________________________________________________
622 void AliAnalysisTaskPWG4PidDetEx::Terminate(Option_t *)
623
624   // Terminate loop
625   if (fDebug > 1) Printf("Terminate()");
626  
627   //
628   // The followig code has now been moved to drawPid.C
629   //
630
631 //   fListOfHists = dynamic_cast<TList*> (GetOutputData(0));
632 //   if (!fListOfHists) {
633 //     Printf("ERROR: fListOfHists not available");
634 //     return;
635 //   }
636  
637 //   TH1I* hevents = dynamic_cast<TH1I*> (fListOfHists->FindObject("fEvents"));
638 //   Int_t nEvents = (Int_t) hevents->GetBinContent(1);
639
640 //   const Float_t normalization = 2.0*fEtaCut*2.0*TMath::Pi();
641
642 //   //-----------------
643 //   TCanvas* c1 = new TCanvas("c1", "c1");
644 //   c1->cd();
645
646 //   TH2F* hMassAll =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassAll"));
647 //   hMassAll->SetLineColor(1);
648 //   hMassAll->SetMarkerColor(1);
649 //   hMassAll->DrawCopy();
650
651 //   TH2F* hMassPion =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassPion"));
652 //   hMassPion->SetLineColor(2);
653 //   hMassPion->SetMarkerColor(2);
654 //   hMassPion->SetMarkerStyle(7);
655 //   hMassPion->DrawCopy("same");
656
657 //   TH2F* hMassKaon =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassKaon"));
658 //   hMassKaon->SetLineColor(3);
659 //   hMassKaon->SetMarkerColor(3);
660 //   hMassKaon->SetMarkerStyle(7);
661 //   hMassKaon->DrawCopy("same");
662
663 //   TH2F* hMassProton =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassProton"));
664 //   hMassProton->SetLineColor(4);
665 //   hMassProton->SetMarkerColor(4);
666 //   hMassProton->SetMarkerStyle(7);
667 //   hMassProton->DrawCopy("same");
668
669 //   TLegend* legend1 = new TLegend(0.8, 0.8, 1, 1);    
670 //   legend1->SetBorderSize(0);
671 //   legend1->SetFillColor(0);
672 //   legend1->AddEntry(hMassAll, "All","LP");
673 //   legend1->AddEntry(hMassPion, "Pions","LP");
674 //   legend1->AddEntry(hMassKaon, "Kaons","LP");
675 //   legend1->AddEntry(hMassProton, "Protons","LP");  
676 //   legend1->Draw();
677
678 //   c1->Update();
679
680 //   //-----------------
681 //   TCanvas* c2 = new TCanvas("c2", "c2");
682 //   c2->cd();
683
684 //   TH2F* hdEdxTPCPion =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCPion"));
685 //   hdEdxTPCPion->SetTitle("dE/dx vs p (TPC)");
686 //   hdEdxTPCPion->SetLineColor(2);
687 //   hdEdxTPCPion->SetMarkerColor(2);
688 //   hdEdxTPCPion->SetMarkerStyle(7);
689 //   hdEdxTPCPion->DrawCopy();
690
691 //   TH2F* hdEdxTPCKaon =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCKaon"));
692 //   hdEdxTPCKaon->SetLineColor(3);
693 //   hdEdxTPCKaon->SetMarkerColor(3);
694 //   hdEdxTPCKaon->SetMarkerStyle(7);
695 //   hdEdxTPCKaon->DrawCopy("same");
696
697 //   TH2F* hdEdxTPCProton =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCProton"));
698 //   hdEdxTPCProton->SetLineColor(4);
699 //   hdEdxTPCProton->SetMarkerColor(4);
700 //   hdEdxTPCProton->SetMarkerStyle(7);
701 //   hdEdxTPCProton->DrawCopy("same");
702
703 //   TLegend* legend2 = new TLegend(0.66, 0.66, 0.88, 0.88);    
704 //   legend2->SetBorderSize(0);
705 //   legend2->SetFillColor(0);
706 //   legend2->AddEntry(hdEdxTPCPion, "Pions","LP");
707 //   legend2->AddEntry(hdEdxTPCKaon, "Kaons","LP");
708 //   legend2->AddEntry(hdEdxTPCProton, "Protons","LP");
709 //   legend2->Draw();
710
711 //   c2->Update();
712
713 //   //-----------------
714 //   TCanvas* c3 = new TCanvas("c3", "c3");
715 //   c3->cd();
716
717 //   TH2F* hdEdxTPCbgPion =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCPion"));
718 //   hdEdxTPCbgPion->SetTitle("dE/dx vs #beta#gamma (TPC)");
719 //   hdEdxTPCbgPion->SetLineColor(2);
720 //   hdEdxTPCbgPion->SetMarkerColor(2);
721 //   hdEdxTPCbgPion->SetMarkerStyle(7);
722 //   hdEdxTPCbgPion->DrawCopy();
723
724 //   TH2F* hdEdxTPCbgKaon =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCKaon"));
725 //   hdEdxTPCbgKaon->SetLineColor(3);
726 //   hdEdxTPCbgKaon->SetMarkerColor(3);
727 //   hdEdxTPCbgKaon->SetMarkerStyle(7);
728 //   hdEdxTPCbgKaon->DrawCopy("same");
729
730 //   TH2F* hdEdxTPCbgProton =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCProton"));
731 //   hdEdxTPCbgProton->SetLineColor(4);
732 //   hdEdxTPCbgProton->SetMarkerColor(4);
733 //   hdEdxTPCbgProton->SetMarkerStyle(7);
734 //   hdEdxTPCbgProton->DrawCopy("same");
735
736 //   TLegend* legend3 = new TLegend(0.66, 0.66, 0.88, 0.88);    
737 //   legend3->SetBorderSize(0);
738 //   legend3->SetFillColor(0);
739 //   legend3->AddEntry(hdEdxTPCbgPion, "Pions","LP");
740 //   legend3->AddEntry(hdEdxTPCbgKaon, "Kaons","LP");
741 //   legend3->AddEntry(hdEdxTPCbgProton, "Protons","LP");
742 //   legend3->Draw();
743
744 //   c3->Update();
745
746 //   //-----------------
747 //   TCanvas* c4 = new TCanvas("c4", "c4", 100, 100, 500, 900);
748 //   c4->Divide(1,3);
749
750 //   c4->cd(1);
751 //   TH1F* hAccepatncePt = dynamic_cast<TH1F*> (fListOfHists->FindObject("fAccPt"));
752 //   TH1F* hAcceptanceP = dynamic_cast<TH1F*> (fListOfHists->FindObject("fAccP"));
753 //   hAccepatncePt->Divide(hAcceptanceP);
754 //   hAccepatncePt->SetTitle("Acceptance correction");
755 //   hAccepatncePt->GetYaxis()->SetTitle("Acceptance correction");
756 //   hAccepatncePt->DrawCopy();
757
758 //   c4->cd(2);
759 //   TH1F* hEfficiencyPID = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffPID"));
760 //   TH1F* hEfficiencyTot = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffTot"));
761 //   hEfficiencyPID->Divide(hEfficiencyTot);
762 //   hEfficiencyPID->SetTitle("Efficiency correction");
763 //   hEfficiencyPID->GetYaxis()->SetTitle("Efficiency correction");
764 //   hEfficiencyPID->DrawCopy();
765
766 //   c4->cd(3);
767 //   TProfile* hKDecayCorr = dynamic_cast<TProfile*> (fListOfHists->FindObject("fKaonDecayCorr"));
768 //   hKDecayCorr->SetTitle("Kaon decay correction");
769 //   hKDecayCorr->GetYaxis()->SetTitle("Kaon decay correction");
770 //   hKDecayCorr->GetXaxis()->SetTitle(" p_{T} [GeV/c]");
771 //   hKDecayCorr->DrawCopy();
772
773 //   c4->Update();
774
775 //   //---------------------
776 //   TCanvas* c5 = new TCanvas("c5", "c5", 100, 100, 600, 900);
777 //   c5->Divide(1,2);
778
779 //   c5->cd(1);
780 //   TH1F* hPtPionNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtPion"));
781 //   hPtPionNoCorr->Scale(1.0/nEvents/normalization/hPtPionNoCorr->GetBinWidth(1));
782 //   hPtPionNoCorr->SetTitle("p_{T} distribution (no corrections)");
783 //   hPtPionNoCorr->SetLineColor(2);
784 //   hPtPionNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
785 //   hPtPionNoCorr->DrawCopy();
786
787 //   TH1F* hPtKaonNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtKaon"));
788 //   hPtKaonNoCorr->Scale(1.0/nEvents/normalization/hPtKaonNoCorr->GetBinWidth(1));
789 //   hPtKaonNoCorr->SetLineColor(3);
790 //   hPtKaonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
791 //   hPtKaonNoCorr->DrawCopy("same");
792
793 //   TH1F* hPtProtonNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtProton"));
794 //   hPtProtonNoCorr->Scale(1.0/nEvents/normalization/hPtProtonNoCorr->GetBinWidth(1));
795 //   hPtProtonNoCorr->SetLineColor(4);
796 //   hPtProtonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
797 //   hPtProtonNoCorr->DrawCopy("same");
798
799 //   TH1F* hPt = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPt"));
800 //   hPt->Scale(1.0/nEvents/normalization/hPt->GetBinWidth(1));
801 //   hPt->GetYaxis()->SetRangeUser(1E-3,1E1);
802 //   hPt->DrawCopy("same");
803
804 //   TLegend* legend4 = new TLegend(0.63, 0.63, 0.88, 0.88);    
805 //   legend4->SetBorderSize(0);
806 //   legend4->SetFillColor(0);
807 //   legend4->AddEntry(hPt, "p_{T} dist (all)","L");
808 //   legend4->AddEntry(hPtPionNoCorr, "p_{T} dist (Pions)","L");
809 //   legend4->AddEntry(hPtKaonNoCorr, "p_{T} dist (Kaons)","L");
810 //   legend4->AddEntry(hPtProtonNoCorr, "p_{T} dist (Protons)","L");
811 //   legend4->Draw();
812
813 //   gPad->SetLogy();
814
815
816 //   c5->cd(2);
817 //   TH1F* hPtPionCorr = static_cast<TH1F*>(hPtPionNoCorr->Clone());
818 //   hPtPionCorr->SetTitle("p_{T} distribution (with corrections)");
819 //   hPtPionCorr->Multiply(hAccepatncePt);
820 //   hPtPionCorr->Divide(hEfficiencyPID);
821 //   hPtPionCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
822 //   hPtPionCorr->DrawCopy();
823
824 //   TH1F* hPtKaonCorr = static_cast<TH1F*>(hPtKaonNoCorr->Clone());
825 //   hPtKaonCorr->Multiply(hAccepatncePt);
826 //   hPtKaonCorr->Divide(hEfficiencyPID);
827 //   hPtKaonCorr->Multiply(hKDecayCorr);
828 //   hPtKaonCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
829 //   hPtKaonCorr->DrawCopy("same");
830
831 //   TH1F* hPtProtonCorr = static_cast<TH1F*>(hPtProtonNoCorr->Clone());
832 //   hPtProtonCorr->Multiply(hAccepatncePt);
833 //   hPtProtonCorr->Divide(hEfficiencyPID);
834 //   hPtProtonCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
835 //   hPtProtonCorr->DrawCopy("same");
836
837 //   hPt->GetYaxis()->SetRangeUser(1E-2,1E1);
838 //   hPt->DrawCopy("same");
839
840 //   TLegend* legend5 = new TLegend(0.63, 0.63, 0.88, 0.88);    
841 //   legend5->SetBorderSize(0);
842 //   legend5->SetFillColor(0);
843 //   legend5->AddEntry(hPt, "p_{T} dist (all)","L");
844 //   legend5->AddEntry(hPtPionCorr, "p_{T} dist (Pions)","L");
845 //   legend5->AddEntry(hPtKaonCorr, "p_{T} dist (Kaons)","L");
846 //   legend5->AddEntry(hPtProtonCorr, "p_{T} dist (Protons)","L");
847 //   legend5->Draw();
848
849 //   gPad->SetLogy();
850
851 //   c5->Update();
852
853 //   //-----------------
854 //   if (fAnalysisMC){
855 //     TCanvas* c6 = new TCanvas("c6", "c6", 100, 100, 1200, 800); 
856 //     c6->Divide(2,2);
857
858 //     c6->cd(1);
859 //     TH1F* hPt_clone = static_cast<TH1F*>(hPt->Clone()); 
860 //     hPt_clone->GetYaxis()->SetRangeUser(1E-5,1E1);
861 //     hPt_clone->SetTitle("p_{T} distribution (all)");
862 //     hPt_clone->SetLineColor(1);
863 //     hPt_clone->DrawCopy();
864  
865 //     TH1F* hPtMC = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMC"));
866 //     hPtMC->Scale(1.0/nEvents/normalization/hPtMC->GetBinWidth(1));
867 //     hPtMC->GetYaxis()->SetRangeUser(1E-5,1E1);
868 //     hPtMC->DrawCopy("same");
869
870 //     TLegend* legend6 = new TLegend(0.57, 0.57, 0.87, 0.87);    
871 //     legend6->SetBorderSize(0);
872 //     legend6->SetFillColor(0);
873 //     legend6->AddEntry(hPt_clone, "p_{T} dist (Rec)", "L");
874 //     legend6->AddEntry(hPtMC, "p_{T} dist (MC)", "L");
875 //     legend6->Draw();
876
877 //     gPad->SetLogy();
878
879 //     c6->cd(2);
880 //     TH1F* hPtPion_clone = static_cast<TH1F*>(hPtPionCorr->Clone()); 
881 //     hPtPion_clone->GetYaxis()->SetRangeUser(1E-2,1E1);
882 //     hPtPion_clone->SetTitle("p_{T} distribution (Pions)");
883 //     hPtPion_clone->SetLineColor(1);
884 //     hPtPion_clone->DrawCopy();
885
886 //     TH1F* hPtMCPion = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCPion"));
887 //     hPtMCPion->Scale(1.0/nEvents/normalization/hPtMCPion->GetBinWidth(1));
888 //     hPtMCPion->GetYaxis()->SetRangeUser(1E-2,1E1);
889 //     hPtMCPion->DrawCopy("same");
890
891 //     TLegend* legend7 = new TLegend(0.57, 0.57, 0.87, 0.87);    
892 //     legend7->SetBorderSize(0);
893 //     legend7->SetFillColor(0);
894 //     legend7->AddEntry(hPtPion_clone, "p_{T} dist (Rec)", "L");
895 //     legend7->AddEntry(hPtMCPion, "p_{T} dist (MC)", "L");
896 //     legend7->Draw();
897
898 //     gPad->SetLogy();
899
900 //     c6->cd(3);
901 //     TH1F* hPtKaon_clone = static_cast<TH1F*>(hPtKaonCorr->Clone()); 
902 //     hPtKaon_clone->GetYaxis()->SetRangeUser(1E-2,1E0);
903 //     hPtKaon_clone->SetLineColor(1);
904 //     hPtKaon_clone->DrawCopy();
905
906 //     TH1F* hPtMCKaon = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCKaon"));
907 //     hPtMCKaon->Scale(1.0/nEvents/normalization/hPtMCKaon->GetBinWidth(1));
908 //     hPtMCKaon->GetYaxis()->SetRangeUser(1E-2,1E0);
909 //     hPtMCKaon->DrawCopy("same");
910
911 //     TLegend* legend8 = new TLegend(0.57, 0.57, 0.87, 0.87);    
912 //     legend8->SetBorderSize(0);
913 //     legend8->SetFillColor(0);
914 //     legend8->AddEntry(hPtKaon_clone, "p_{T} dist (Rec)", "L");
915 //     legend8->AddEntry(hPtMCKaon, "p_{T} dist (MC)", "L");
916 //     legend8->Draw();
917
918 //     gPad->SetLogy(); 
919
920 //     c6->cd(4);
921 //     TH1F* hPtProton_clone = static_cast<TH1F*>(hPtProtonCorr->Clone()); 
922 //     hPtProton_clone->GetYaxis()->SetRangeUser(1E-2,1E-1);
923 //     hPtProton_clone->SetLineColor(1);
924 //     hPtProton_clone->DrawCopy();
925
926 //     TH1F* hPtMCProton = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCProton"));
927 //     hPtMCProton->Scale(1.0/nEvents/normalization/hPtMCProton->GetBinWidth(1));
928 //     hPtMCProton->GetYaxis()->SetRangeUser(1E-2,1E-1);
929 //     hPtMCProton->DrawCopy("same");
930
931 //     TLegend* legend9 = new TLegend(0.2, 0.25, 0.5, 0.55);    
932 //     legend9->SetBorderSize(0);
933 //     legend9->SetFillColor(0);
934 //     legend9->AddEntry(hPtProton_clone, "p_{T} dist (Rec)", "L");
935 //     legend9->AddEntry(hPtMCProton, "p_{T} dist (MC)", "L");
936 //     legend9->Draw();
937
938 //     gPad->SetLogy(); 
939
940 //     c6->Update();
941 //   }
942
943 }