Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDAnalysisTask.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 //==============================================================================
17 // AliHMPIDAnalysysTask - Class representing a basic analysis tool of HMPID data
18 // A set of histograms is created.
19 //==============================================================================
20 //
21 // By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
22 // to perform the analysis on local data, on data on alien using local machine
23 // and on CAF.
24
25 #ifndef AliHMPIDAnalysisTASK_CXX
26 #define AliHMPIDAnalysisTASK_CXX
27
28
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TFile.h"
32 #include "TCanvas.h"
33 #include "TGraphErrors.h"
34 #include "AliAnalysisManager.h"
35 #include "AliESDInputHandler.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliESDtrack.h"
39 #include "AliPID.h"
40 #include "AliLog.h"
41 #include "AliHMPIDAnalysisTask.h"
42
43 ClassImp(AliHMPIDAnalysisTask)
44
45 //__________________________________________________________________________
46 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
47   fESD(0x0),fMC(0x0),fUseMC(kTRUE),
48   fHmpHistList(0x0),
49   fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
50   fHmpMipTrkDist(0x0),fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),
51   fHmpMipCharge3cm(0x0),fHmpMipCharge1cm(0x0),
52   fHmpNumPhots(0x0),fHmpTrkFlags(0x0),
53   fN1(6),
54   fN2(8),
55   fPionEff(0x0),
56   fKaonEff(0x0),
57   fProtEff(0x0),
58   fPionTot(0x0),
59   fKaonTot(0x0),
60   fProtTot(0x0),
61   fPionNot(0x0),
62   fKaonNot(0x0),
63   fProtNot(0x0),
64   fPionCon(0x0),
65   fKaonCon(0x0),
66   fProtCon(0x0),
67   fThetavsPiFromK(0x0),
68   fThetapivsPesd(0x0),
69   fThetaKvsPesd(0x0),
70   fThetaPvsPesd(0x0),
71   fProtGen(0x0),
72   fPbarGen(0x0),
73   fProtHmp(0x0),
74   fPbarHmp(0x0),
75   fTree(0x0)
76 {
77   //
78   //Default ctor
79   //
80   for (Int_t i=0; i<28; i++) fVar[i]=0;
81 }
82
83 //___________________________________________________________________________
84 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
85   AliAnalysisTaskSE(name),
86   fESD(0x0), fMC(0x0), fUseMC(kTRUE),
87   fHmpHistList(0x0),
88   fHmpPesdPhmp(0x0), fHmpCkovPesd(0x0), fHmpCkovPhmp(0x0),
89   fHmpMipTrkDist(0x0), fHmpMipTrkDistX(0x0), fHmpMipTrkDistY(0x0),
90   fHmpMipCharge3cm(0x0), fHmpMipCharge1cm(0x0),
91   fHmpNumPhots(0x0), fHmpTrkFlags(0x0),
92   fN1(6),
93   fN2(8),
94   fPionEff(0x0),
95   fKaonEff(0x0),
96   fProtEff(0x0),
97   fPionTot(0x0),
98   fKaonTot(0x0),
99   fProtTot(0x0),
100   fPionNot(0x0),
101   fKaonNot(0x0),
102   fProtNot(0x0),
103   fPionCon(0x0),
104   fKaonCon(0x0),
105   fProtCon(0x0),
106   fThetavsPiFromK(0x0),
107   fThetapivsPesd(0x0),
108   fThetaKvsPesd(0x0),
109   fThetaPvsPesd(0x0),
110   fProtGen(0x0),
111   fPbarGen(0x0),
112   fProtHmp(0x0),
113   fPbarHmp(0x0),
114   fTree(0x0)
115 {
116   //
117   // Constructor. Initialization of Inputs and Outputs
118   //
119   for (Int_t i=0; i<28; i++) fVar[i]=0;
120
121   DefineOutput(1,TList::Class());
122   DefineOutput(2,TTree::Class());
123 }
124
125 //___________________________________________________________________________
126 AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
127 {
128   //
129   // Assignment operator
130   //
131   if (this!=&c) {
132     AliAnalysisTaskSE::operator=(c);
133     fESD             = c.fESD;
134     fMC              = c.fMC;
135     fUseMC           = c.fUseMC;
136     fHmpHistList     = c.fHmpHistList;
137     fHmpPesdPhmp     = c.fHmpPesdPhmp;
138     fHmpCkovPesd     = c.fHmpCkovPesd;
139     fHmpCkovPhmp     = c.fHmpCkovPhmp;
140     fHmpMipTrkDist   = c.fHmpMipTrkDist;
141     fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
142     fHmpMipTrkDistY  = c.fHmpMipTrkDistY;
143     fHmpMipCharge3cm = c.fHmpMipCharge3cm;
144     fHmpMipCharge1cm = c.fHmpMipCharge1cm;
145     fHmpNumPhots     = c.fHmpNumPhots;
146     fHmpTrkFlags     = c.fHmpTrkFlags;
147     fN1              = c.fN1;
148     fN2              = c.fN2;
149     fPionEff         = c.fPionEff;
150     fKaonEff         = c.fKaonEff;
151     fProtEff         = c.fProtEff;
152     fPionTot         = c.fPionTot;
153     fKaonTot         = c.fKaonTot;
154     fProtTot         = c.fProtTot;
155     fPionNot         = c.fPionNot;
156     fKaonNot         = c.fKaonNot;
157     fProtNot         = c.fProtNot;
158     fPionCon         = c.fPionCon;
159     fKaonCon         = c.fKaonCon;
160     fProtCon         = c.fProtCon;
161     fThetavsPiFromK  = c.fThetavsPiFromK;
162     fThetapivsPesd   = c.fThetapivsPesd;
163     fThetaKvsPesd    = c.fThetaKvsPesd;
164     fThetaPvsPesd    = c.fThetaPvsPesd;
165     fProtGen         = c.fProtGen;
166     fPbarGen         = c.fPbarGen;
167     fProtHmp         = c.fProtHmp;
168     fPbarHmp         = c.fPbarHmp;
169     fTree            = c.fTree;
170     for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
171   }
172   return *this;
173 }
174
175 //___________________________________________________________________________
176 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
177   AliAnalysisTaskSE(c),
178   fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
179   fHmpHistList(c.fHmpHistList),
180   fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
181   fHmpMipTrkDist(c.fHmpMipTrkDist),fHmpMipTrkDistX(c.fHmpMipTrkDistX),fHmpMipTrkDistY(c.fHmpMipTrkDistY),
182   fHmpMipCharge3cm(c.fHmpMipCharge3cm), fHmpMipCharge1cm(c.fHmpMipCharge1cm),
183   fHmpNumPhots(c.fHmpNumPhots), fHmpTrkFlags(c.fHmpTrkFlags),
184   fN1(c.fN1),
185   fN2(c.fN2),
186   fPionEff(c.fPionEff),
187   fKaonEff(c.fKaonEff),
188   fProtEff(c.fProtEff),
189   fPionTot(c.fPionTot),
190   fKaonTot(c.fKaonTot),
191   fProtTot(c.fProtTot),
192   fPionNot(c.fPionNot),
193   fKaonNot(c.fKaonNot),
194   fProtNot(c.fProtNot),
195   fPionCon(c.fPionCon),
196   fKaonCon(c.fKaonCon),
197   fProtCon(c.fProtCon),
198   fThetavsPiFromK(c.fThetavsPiFromK),
199   fThetapivsPesd(c.fThetapivsPesd),
200   fThetaKvsPesd(c.fThetaKvsPesd),
201   fThetaPvsPesd(c.fThetaPvsPesd),
202   fProtGen(c.fProtGen),
203   fPbarGen(c.fPbarGen),
204   fProtHmp(c.fProtHmp),
205   fPbarHmp(c.fPbarHmp),
206   fTree(c.fTree)
207 {
208   //
209   // Copy Constructor
210   //
211   for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
212 }
213  
214 //___________________________________________________________________________
215 AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
216   //
217   //destructor
218   //
219   Info("~AliHMPIDAnalysisTask","Calling Destructor");
220   if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
221 }
222
223 //___________________________________________________________________________
224 void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
225 {
226   // Connect ESD here
227
228   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
229   if (!esdH) {
230     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
231   } else
232     fESD = esdH->GetEvent();
233
234   if (fUseMC){
235     // Connect MC
236     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
237     if (!mcH) {
238       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
239       fUseMC = kFALSE;
240     } else
241       fMC = mcH->MCEvent();
242       if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
243   }
244 }
245
246 //___________________________________________________________________________
247 void AliHMPIDAnalysisTask::UserExec(Option_t *)
248 {
249   Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
250   Double_t probs[5];
251   AliPID *pPid = new AliPID();
252   pPid->SetPriors(priors);
253   AliESDtrack *track = 0;
254   TParticle *pPart = 0;
255   AliStack* pStack = 0;
256   Int_t label = -1;
257   if (fUseMC){
258     pStack = fMC->Stack();
259   }
260
261   Double_t ktol = 0.001;
262   Double_t rin[3], rout[3];
263   Float_t distCut  = 0.7;
264   Float_t thetaCut = 0.175;
265
266   //
267   // Main loop function, executed on Event basis
268   //
269   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
270
271     track = fESD->GetTrack(iTrack);
272     if(!track) continue;
273     track->GetInnerXYZ(rin);
274     track->GetOuterXYZ(rout);
275
276     if(Equal(track->GetHMPIDsignal(),-20.,ktol))      fHmpTrkFlags->Fill(0);
277     else if(Equal(track->GetHMPIDsignal(),-9.,ktol))  fHmpTrkFlags->Fill(1);
278     else if(Equal(track->GetHMPIDsignal(),-5.,ktol))  fHmpTrkFlags->Fill(2);
279     else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
280     else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
281     else fHmpTrkFlags->Fill(5);
282
283     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
284     if(track->GetHMPIDcluIdx() < 0) continue;
285
286     Int_t q, nph;
287     Float_t x, y;
288     Float_t xpc, ypc, th, ph;
289     track->GetHMPIDmip(x,y,q,nph);
290     track->GetHMPIDtrk(xpc,ypc,th,ph);
291
292     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
293
294     Float_t dist = TMath::Sqrt((xpc-x)*(xpc-x)+(ypc-y)*(ypc-y));
295     fHmpMipTrkDist->Fill(dist);
296     fHmpMipTrkDistX->Fill(xpc - x);
297     fHmpMipTrkDistY->Fill(ypc - y);
298     Double_t pHmp[3] = {0}, pHmp3 = 0;
299     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
300     if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
301
302     Float_t b[2];
303     Float_t bCov[3];
304     track->GetImpactParameters(b,bCov);    
305
306     track->GetHMPIDpid(probs);
307     pPid->SetProbabilities(probs);
308     if (fUseMC){
309       if ((label = track->GetLabel()) < 0) continue;
310       pPart = pStack->Particle(label);
311     }
312
313     if(track->GetHMPIDsignal() > 0 ){
314       Double_t thetaCkov = track->GetHMPIDsignal() - (Int_t)track->GetHMPIDsignal();
315       if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
316       if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
317       fHmpNumPhots->Fill(nph);
318       fHmpCkovPesd->Fill(track->P(),thetaCkov);
319       if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,thetaCkov);
320
321       if (fUseMC && dist<distCut && TMath::Abs(th)<thetaCut){
322         if (!pStack->IsPhysicalPrimary(label)) continue;
323         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
324         if (pdgCode==211){
325           fThetapivsPesd->Fill(track->P(),thetaCkov);
326           Int_t mot=pPart->GetFirstMother();
327           if (mot > -1){
328             TParticle *pMot=pStack->Particle(mot);
329             TString str=pMot->GetName();
330             if (str.Contains("K")) fThetavsPiFromK->Fill(pHmp3,thetaCkov);
331           }
332         }
333         if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),thetaCkov);
334         if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),thetaCkov);
335
336         if (track->P()<1. || track->P()>5.) continue;
337         Int_t pBin=(Int_t) (2*(track->P()-1));
338         if (pdgCode!=2212) fProtCon->Fill(pBin);
339         if (pdgCode==2212){
340           fProtTot->Fill(pBin);
341           fProtEff->Fill(pBin,pPid->GetProbability(AliPID::kProton));
342           fPionNot->Fill(pBin,pPid->GetProbability(AliPID::kPion));
343           fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
344         }
345         if (pdgCode!=211) fPionCon->Fill(pBin);
346         if (pdgCode!=321) fKaonCon->Fill(pBin);
347         if (pdgCode==211){
348           if (pBin < 6){
349             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
350                            pPid->GetProbability(AliPID::kMuon)+
351                            pPid->GetProbability(AliPID::kPion);
352             fPionTot->Fill(pBin);
353             fPionEff->Fill(pBin,weight);
354             fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
355           }
356           fProtNot->Fill(pBin,pPid->GetProbability(AliPID::kProton));
357         }
358         if (pdgCode==321){
359           if (pBin < 6){
360             fKaonTot->Fill(pBin);
361             fKaonEff->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
362             fPionNot->Fill(pBin,(pPid->GetProbability(AliPID::kPion)));
363           }
364           fProtNot->Fill(pBin,(pPid->GetProbability(AliPID::kProton)));
365         }
366       }
367     }//there is signal
368     fVar[0] = track->GetHMPIDcluIdx()/1000000;
369     fVar[1] = pHmp3;
370     fVar[2] = (Float_t)track->P();
371     fVar[3] = xpc;
372     fVar[4] = ypc;
373     fVar[5] = x;
374     fVar[6] = y;
375     fVar[7] = (Float_t)thetaCkov;
376     fVar[8] = q;
377     fVar[9] = th;
378     fVar[10] = ph;
379     fVar[11] = (Float_t)track->GetSign();
380     fVar[12] = (Float_t)nph;
381     fVar[13] = (Float_t)track->GetNcls(1);
382     fVar[14] = (Float_t)probs[0];
383     fVar[15] = (Float_t)probs[1];
384     fVar[16] = (Float_t)probs[2];
385     fVar[17] = (Float_t)probs[3];
386     fVar[18] = (Float_t)probs[4];
387     fVar[19] = (Float_t)track->GetTOFsignal();
388     fVar[20] = (Float_t)track->GetKinkIndex(0);
389     fVar[21] = (Float_t)track->Xv();
390     fVar[22] = (Float_t)track->Yv();
391     fVar[23] = (Float_t)track->Zv();
392     fVar[24] = (Float_t)track->GetTPCchi2();
393     fVar[25] = b[0];
394     fVar[26] = b[1];
395     fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
396     fTree->Fill();
397   }//track loop
398   delete pPid;
399
400   if (fUseMC){
401     Float_t phiMin=(5./180.)*TMath::Pi() , phiMax=(55./180.)*TMath::Pi();
402     for (Int_t iPart=0; iPart<pStack->GetNprimary(); iPart++){
403       TString namepart=pStack->Particle(iPart)->GetName();
404       if (!namepart.Contains("proton")) continue;
405       pPart=pStack->Particle(iPart);
406       if (!pStack->IsPhysicalPrimary(iPart) || pPart->P()<1. || pPart->P()>5.) continue;
407       if (pPart->Phi()< phiMin || pPart->Phi()>phiMax) continue;
408       if (TMath::Abs(pPart->Eta()) > 0.5) continue;
409       Int_t pBin=(Int_t) (2*(pPart->P()-1));
410       if (namepart == "proton") fProtGen->Fill(pBin);
411       if (namepart == "antiproton") fPbarGen->Fill(pBin);
412       for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++){
413         if (fESD->GetTrack(iTrack)->GetLabel() != iPart) continue;
414         if (fESD->GetTrack(iTrack)->GetHMPIDsignal() < 0 ) continue;
415         if (namepart == "proton") fProtHmp->Fill(pBin);
416         if (namepart == "antiproton") fPbarHmp->Fill(pBin);
417       }
418     }
419   }
420
421   /* PostData(0) is taken care of by AliAnalysisTaskSE */
422   PostData(1,fHmpHistList);
423   PostData(2,fTree);
424 }
425
426 //___________________________________________________________________________
427 void AliHMPIDAnalysisTask::Terminate(Option_t*)
428 {
429   // The Terminate() function is the last function to be called during
430   // a query. It always runs on the client, it can be used to present
431   // the results graphically or save the results to file.
432
433   Info("Terminate"," ");
434
435   if (!fUseMC) return;
436
437   fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
438
439   if (!fHmpHistList) {
440     AliError("Histogram List is not available");
441     return;
442   }
443
444   fPionEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionEff"));
445   fPionTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionTot"));
446   fPionNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionNot"));
447   fPionCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionCon"));
448   fKaonEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonEff"));
449   fKaonTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonTot"));
450   fKaonNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonNot"));
451   fKaonCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonCon"));
452   fProtEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtEff"));
453   fProtTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtTot"));
454   fProtNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtNot"));
455   fProtCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtCon"));
456
457   Float_t *pionEff = fPionEff->GetArray();
458   Int_t   *pionTot = fPionTot->GetArray();
459   Float_t *pionNot = fPionNot->GetArray();
460   Int_t   *pionCon = fPionCon->GetArray();
461   Float_t *kaonEff = fKaonEff->GetArray();
462   Int_t   *kaonTot = fKaonTot->GetArray();
463   Float_t *kaonNot = fKaonNot->GetArray();
464   Int_t   *kaonCon = fKaonCon->GetArray();
465   Float_t *protEff = fProtEff->GetArray();
466   Int_t   *protTot = fProtTot->GetArray();
467   Float_t *protNot = fProtNot->GetArray();
468   Int_t   *protCon = fProtCon->GetArray();
469
470   TGraphErrors *effPi = new TGraphErrors(fN1);
471   TGraphErrors *effKa = new TGraphErrors(fN1);
472   TGraphErrors *effPr = new TGraphErrors(fN2);
473   TGraphErrors *conPi = new TGraphErrors(fN1);
474   TGraphErrors *conKa = new TGraphErrors(fN1);
475   TGraphErrors *conPr = new TGraphErrors(fN2);
476
477   Float_t pt[8]={1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75};
478   Float_t eff=0, effErr=0, con=0, conErr=0;
479   for (Int_t i=0; i< fN2; i++){
480     eff = protEff[i+1]/TMath::Max(protTot[i+1],1);
481     effErr = TMath::Sqrt(eff*(1.-eff)/TMath::Max(protTot[i+1],1));
482     con = protNot[i+1]/TMath::Max(protCon[i+1],1);
483     conErr = TMath::Sqrt(con*(1.-con)/protCon[i+1]);
484     effPr->SetPoint(i,pt[i],eff);
485     effPr->SetPointError(i,0,effErr);
486     conPr->SetPoint(i,pt[i],con);
487     conPr->SetPointError(i,0,conErr);
488
489     if (i>=fN1) continue;
490     eff = pionEff[i+1]/pionTot[i+1];
491     effErr = TMath::Sqrt(eff*(1.-eff)/pionTot[i+1]);
492     con=pionNot[i+1]/pionCon[i+1];
493     conErr = TMath::Sqrt(con*(1.-con)/pionCon[i+1]);
494     effPi->SetPoint(i,pt[i],(Float_t)pionEff[i+1]/(Float_t)pionTot[i+1]);
495     effPi->SetPointError(i,0,effErr);
496     conPi->SetPoint(i,pt[i],(Float_t)pionNot[i+1]/(Float_t)pionCon[i+1]);
497     conPi->SetPointError(i,0,conErr);
498
499     eff = kaonEff[i+1]/TMath::Max(kaonTot[i+1],1);
500     effErr = TMath::Sqrt(eff*(1.-eff)/kaonTot[i+1]);
501     con = kaonNot[i+1]/TMath::Max(kaonCon[i+1],1);
502     conErr = TMath::Sqrt(con*(1.-con)/kaonCon[i+1]);
503     effKa->SetPoint(i,pt[i],(Float_t)kaonEff[i+1]/TMath::Max(kaonTot[i+1],1));
504     effKa->SetPointError(i,0,effErr);
505     conKa->SetPoint(i,pt[i],(Float_t)kaonNot[i+1]/TMath::Max(kaonCon[i+1],1));
506     conKa->SetPointError(i,0,conErr);
507   }
508
509   fProtGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtGen"));
510   fPbarGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarGen"));
511   fProtHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtHmp"));
512   fPbarHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarHmp"));
513
514   Int_t   *protGen = fProtGen->GetArray();
515   Int_t   *pbarGen = fPbarGen->GetArray();
516   Int_t   *protHmp = fProtHmp->GetArray();
517   Int_t   *pbarHmp = fPbarHmp->GetArray();
518
519   TGraphErrors *ratPr = new TGraphErrors(fN2);
520   TGraphErrors *ratPb = new TGraphErrors(fN2);
521
522   Float_t rat1=0, rat1Err=0, rat2=0, rat2Err=0;
523   for (Int_t i=0; i< fN2; i++){
524     rat1 = (Float_t)protHmp[i+1]/TMath::Max(protGen[i+1],1);
525     rat1Err = TMath::Sqrt(rat1*(1.-rat1)/TMath::Max(protGen[i+1],1));
526     rat2 = (Float_t)pbarHmp[i+1]/TMath::Max(pbarGen[i+1],1);
527     rat2Err = TMath::Sqrt(rat2*(1.-rat2)/TMath::Max(pbarGen[i+1],1));
528     ratPr->SetPoint(i,pt[i],rat1);
529     ratPr->SetPointError(i,0,rat1Err);
530     ratPb->SetPoint(i,pt[i],rat2);
531     ratPb->SetPointError(i,0,rat2Err);
532   }
533
534   TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
535   pCan->Divide(2,3);
536
537   pCan->cd(1);
538   effPi->SetTitle("Pions");
539   effPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
540   effPi->GetYaxis()->SetRangeUser(0.,1.);
541   effPi->SetMarkerStyle(20);
542   effPi->Draw("ALP");
543   conPi->SetMarkerStyle(21);
544   conPi->Draw("sameLP");
545
546   pCan->cd(3);
547   effKa->SetTitle("Kaons");
548   effKa->GetXaxis()->SetTitle("p_{T} (GeV/c)");
549   effKa->GetYaxis()->SetRangeUser(0.,1.);
550   effKa->SetMarkerStyle(20);
551   effKa->Draw("ALP");
552   conKa->SetMarkerStyle(21);
553   conKa->Draw("sameLP");
554
555   pCan->cd(5);
556   effPr->SetTitle("Protons");
557   effPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
558   effPr->GetYaxis()->SetRangeUser(0.,1.);
559   effPr->SetMarkerStyle(20);
560   effPr->Draw("ALP");
561   conPr->SetMarkerStyle(21);
562   conPr->Draw("sameLP");
563
564   pCan->cd(2);
565   ratPr->SetTitle("Ratios");
566   ratPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
567   ratPr->GetYaxis()->SetRangeUser(0.,1.);
568   ratPr->SetMarkerStyle(20);
569   ratPr->Draw("ALP");
570   ratPb->SetMarkerStyle(21);
571   ratPb->Draw("sameLP");
572
573   TFile *outFile = new TFile("HmpidGraphs.root","recreate");
574   pCan->Write();
575   outFile->Close();
576
577   AliAnalysisTaskSE::Terminate();
578
579 }
580
581 //___________________________________________________________________________
582 void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
583   //
584   //HERE ONE CAN CREATE OUTPUT OBJECTS
585   //TO BE SET BEFORE THE EXECUTION OF THE TASK
586   //
587
588   //slot #1
589 //   OpenFile(1);
590    fHmpHistList = new TList();
591    fHmpHistList->SetOwner();
592
593    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
594    fHmpHistList->Add(fHmpPesdPhmp);
595
596    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
597    fHmpHistList->Add(fHmpCkovPesd);
598
599    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
600    fHmpHistList->Add(fHmpCkovPhmp);
601
602    fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",400,0,20);
603    fHmpHistList->Add(fHmpMipTrkDist);
604
605    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
606    fHmpHistList->Add(fHmpMipTrkDistX);
607
608    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
609    fHmpHistList->Add(fHmpMipTrkDistY);
610
611    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
612    fHmpHistList->Add(fHmpMipCharge3cm);
613
614    fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
615    fHmpHistList->Add(fHmpMipCharge1cm);
616
617    fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
618    fHmpHistList->Add(fHmpNumPhots);
619
620    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
621    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
622    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
623    fHmpHistList->Add(fHmpTrkFlags);
624
625    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
626    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
627    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
628    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
629    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
630    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
631    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
632    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
633    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
634    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
635    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
636    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
637
638    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
639    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
640    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
641    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
642
643    fProtGen = new TH1I("ProtGen","Generated protons",fN2,0,fN2);
644    fPbarGen = new TH1I("PbarGen","Generated antiprotons",fN2,0,fN2);
645    fProtHmp = new TH1I("ProtHmp","HMPID protons",fN2,0,fN2);
646    fPbarHmp = new TH1I("PbarHmp","HMPID antiprotons",fN2,0,fN2);
647
648    fHmpHistList->Add(fProtGen); fHmpHistList->Add(fPbarGen);
649    fHmpHistList->Add(fProtHmp); fHmpHistList->Add(fPbarHmp);
650
651    fThetavsPiFromK= new TH2F("ThetavsPiFromK","Theta vs p of pions from K;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
652    fHmpHistList->Add(fThetavsPiFromK);
653
654    fThetapivsPesd = new TH2F("ThetapivsPesd","Theta of pions vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
655    fHmpHistList->Add(fThetapivsPesd);
656
657    fThetaKvsPesd  = new TH2F("ThetaKvsPesd","Theta of kaons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
658    fHmpHistList->Add(fThetaKvsPesd);
659
660    fThetaPvsPesd  = new TH2F("ThetaPvsPesd","Theta of protons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
661    fHmpHistList->Add(fThetaPvsPesd);
662
663 //   OpenFile(2);
664    fTree = new TTree("Tree","Tree with data");
665    fTree->Branch("Chamber",&fVar[0]);
666    fTree->Branch("pHmp3",&fVar[1]);
667    fTree->Branch("P",&fVar[2]);
668    fTree->Branch("Xpc",&fVar[3]);
669    fTree->Branch("Ypc",&fVar[4]);
670    fTree->Branch("X",&fVar[5]);
671    fTree->Branch("Y",&fVar[6]);
672    fTree->Branch("HMPIDsignal",&fVar[7]);
673    fTree->Branch("Charge",&fVar[8]);
674    fTree->Branch("Theta",&fVar[9]);
675    fTree->Branch("Phi",&fVar[10]);
676    fTree->Branch("Sign",&fVar[11]);
677    fTree->Branch("NumPhotons",&fVar[12]);
678    fTree->Branch("NumTPCclust",&fVar[13]);
679    fTree->Branch("Prob0",&fVar[14]);
680    fTree->Branch("Prob1",&fVar[15]);
681    fTree->Branch("Prob2",&fVar[16]);
682    fTree->Branch("Prob3",&fVar[17]);
683    fTree->Branch("Prob4",&fVar[18]);
684    fTree->Branch("TOFsignal",&fVar[19]);
685    fTree->Branch("KinkIndex",&fVar[20]);
686    fTree->Branch("Xv",&fVar[21]);
687    fTree->Branch("Yv",&fVar[22]);
688    fTree->Branch("Zv",&fVar[23]);
689    fTree->Branch("TPCchi2",&fVar[24]);
690    fTree->Branch("b0",&fVar[25]);
691    fTree->Branch("b1",&fVar[26]);
692    fTree->Branch("ClustSize",&fVar[27]);
693
694    PostData(1,fHmpHistList);
695    PostData(2,fTree);
696 }
697
698 //____________________________________________________________________________________________________________________________________
699 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
700 {
701  return abs(x - y) <= tolerance ;
702 }
703    
704 #endif