Deleted wrong file
[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       if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
315       if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
316       fHmpNumPhots->Fill(nph);
317       fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
318       if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
319
320       if (fUseMC && dist<distCut && TMath::Abs(th)<thetaCut){
321         if (!pStack->IsPhysicalPrimary(label)) continue;
322         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
323         if (pdgCode==211){
324           fThetapivsPesd->Fill(track->P(),track->GetHMPIDsignal());
325           Int_t mot=pPart->GetFirstMother();
326           if (mot > -1){
327             TParticle *pMot=pStack->Particle(mot);
328             TString str=pMot->GetName();
329             if (str.Contains("K")) fThetavsPiFromK->Fill(pHmp3,track->GetHMPIDsignal());
330           }
331         }
332         if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),track->GetHMPIDsignal());
333         if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),track->GetHMPIDsignal());
334
335         if (track->P()<1. || track->P()>5.) continue;
336         Int_t pBin=(Int_t) (2*(track->P()-1));
337         if (pdgCode!=2212) fProtCon->Fill(pBin);
338         if (pdgCode==2212){
339           fProtTot->Fill(pBin);
340           fProtEff->Fill(pBin,pPid->GetProbability(AliPID::kProton));
341           fPionNot->Fill(pBin,pPid->GetProbability(AliPID::kPion));
342           fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
343         }
344         if (pdgCode!=211) fPionCon->Fill(pBin);
345         if (pdgCode!=321) fKaonCon->Fill(pBin);
346         if (pdgCode==211){
347           if (pBin < 6){
348             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
349                            pPid->GetProbability(AliPID::kMuon)+
350                            pPid->GetProbability(AliPID::kPion);
351             fPionTot->Fill(pBin);
352             fPionEff->Fill(pBin,weight);
353             fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
354           }
355           fProtNot->Fill(pBin,pPid->GetProbability(AliPID::kProton));
356         }
357         if (pdgCode==321){
358           if (pBin < 6){
359             fKaonTot->Fill(pBin);
360             fKaonEff->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
361             fPionNot->Fill(pBin,(pPid->GetProbability(AliPID::kPion)));
362           }
363           fProtNot->Fill(pBin,(pPid->GetProbability(AliPID::kProton)));
364         }
365       }
366     }//there is signal
367     fVar[0] = track->GetHMPIDcluIdx()/1000000;
368     fVar[1] = pHmp3;
369     fVar[2] = (Float_t)track->P();
370     fVar[3] = xpc;
371     fVar[4] = ypc;
372     fVar[5] = x;
373     fVar[6] = y;
374     fVar[7] = (Float_t)track->GetHMPIDsignal();
375     fVar[8] = q;
376     fVar[9] = th;
377     fVar[10] = ph;
378     fVar[11] = (Float_t)track->GetSign();
379     fVar[12] = (Float_t)nph;
380     fVar[13] = (Float_t)track->GetNcls(1);
381     fVar[14] = (Float_t)probs[0];
382     fVar[15] = (Float_t)probs[1];
383     fVar[16] = (Float_t)probs[2];
384     fVar[17] = (Float_t)probs[3];
385     fVar[18] = (Float_t)probs[4];
386     fVar[19] = (Float_t)track->GetTOFsignal();
387     fVar[20] = (Float_t)track->GetKinkIndex(0);
388     fVar[21] = (Float_t)track->Xv();
389     fVar[22] = (Float_t)track->Yv();
390     fVar[23] = (Float_t)track->Zv();
391     fVar[24] = (Float_t)track->GetTPCchi2();
392     fVar[25] = b[0];
393     fVar[26] = b[1];
394     fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
395     fTree->Fill();
396   }//track loop
397   delete pPid;
398
399   if (fUseMC){
400     Float_t phiMin=(5./180.)*TMath::Pi() , phiMax=(55./180.)*TMath::Pi();
401     for (Int_t iPart=0; iPart<pStack->GetNprimary(); iPart++){
402       TString namepart=pStack->Particle(iPart)->GetName();
403       if (!namepart.Contains("proton")) continue;
404       pPart=pStack->Particle(iPart);
405       if (!pStack->IsPhysicalPrimary(iPart) || pPart->P()<1. || pPart->P()>5.) continue;
406       if (pPart->Phi()< phiMin || pPart->Phi()>phiMax) continue;
407       if (TMath::Abs(pPart->Eta()) > 0.5) continue;
408       Int_t pBin=(Int_t) (2*(pPart->P()-1));
409       if (namepart == "proton") fProtGen->Fill(pBin);
410       if (namepart == "antiproton") fPbarGen->Fill(pBin);
411       for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++){
412         if (fESD->GetTrack(iTrack)->GetLabel() != iPart) continue;
413         if (fESD->GetTrack(iTrack)->GetHMPIDsignal() < 0 ) continue;
414         if (namepart == "proton") fProtHmp->Fill(pBin);
415         if (namepart == "antiproton") fPbarHmp->Fill(pBin);
416       }
417     }
418   }
419
420   /* PostData(0) is taken care of by AliAnalysisTaskSE */
421   PostData(1,fHmpHistList);
422   PostData(2,fTree);
423 }
424
425 //___________________________________________________________________________
426 void AliHMPIDAnalysisTask::Terminate(Option_t*)
427 {
428   // The Terminate() function is the last function to be called during
429   // a query. It always runs on the client, it can be used to present
430   // the results graphically or save the results to file.
431
432   Info("Terminate"," ");
433
434   if (!fUseMC) return;
435
436   fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
437
438   if (!fHmpHistList) {
439     AliError("Histogram List is not available");
440     return;
441   }
442
443   fPionEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionEff"));
444   fPionTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionTot"));
445   fPionNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionNot"));
446   fPionCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionCon"));
447   fKaonEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonEff"));
448   fKaonTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonTot"));
449   fKaonNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonNot"));
450   fKaonCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonCon"));
451   fProtEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtEff"));
452   fProtTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtTot"));
453   fProtNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtNot"));
454   fProtCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtCon"));
455
456   Float_t *pionEff = fPionEff->GetArray();
457   Int_t   *pionTot = fPionTot->GetArray();
458   Float_t *pionNot = fPionNot->GetArray();
459   Int_t   *pionCon = fPionCon->GetArray();
460   Float_t *kaonEff = fKaonEff->GetArray();
461   Int_t   *kaonTot = fKaonTot->GetArray();
462   Float_t *kaonNot = fKaonNot->GetArray();
463   Int_t   *kaonCon = fKaonCon->GetArray();
464   Float_t *protEff = fProtEff->GetArray();
465   Int_t   *protTot = fProtTot->GetArray();
466   Float_t *protNot = fProtNot->GetArray();
467   Int_t   *protCon = fProtCon->GetArray();
468
469   TGraphErrors *effPi = new TGraphErrors(fN1);
470   TGraphErrors *effKa = new TGraphErrors(fN1);
471   TGraphErrors *effPr = new TGraphErrors(fN2);
472   TGraphErrors *conPi = new TGraphErrors(fN1);
473   TGraphErrors *conKa = new TGraphErrors(fN1);
474   TGraphErrors *conPr = new TGraphErrors(fN2);
475
476   Float_t pt[8]={1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75};
477   Float_t eff=0, effErr=0, con=0, conErr=0;
478   for (Int_t i=0; i< fN2; i++){
479     eff = protEff[i+1]/TMath::Max(protTot[i+1],1);
480     effErr = TMath::Sqrt(eff*(1.-eff)/TMath::Max(protTot[i+1],1));
481     con = protNot[i+1]/TMath::Max(protCon[i+1],1);
482     conErr = TMath::Sqrt(con*(1.-con)/protCon[i+1]);
483     effPr->SetPoint(i,pt[i],eff);
484     effPr->SetPointError(i,0,effErr);
485     conPr->SetPoint(i,pt[i],con);
486     conPr->SetPointError(i,0,conErr);
487
488     if (i>=fN1) continue;
489     eff = pionEff[i+1]/pionTot[i+1];
490     effErr = TMath::Sqrt(eff*(1.-eff)/pionTot[i+1]);
491     con=pionNot[i+1]/pionCon[i+1];
492     conErr = TMath::Sqrt(con*(1.-con)/pionCon[i+1]);
493     effPi->SetPoint(i,pt[i],(Float_t)pionEff[i+1]/(Float_t)pionTot[i+1]);
494     effPi->SetPointError(i,0,effErr);
495     conPi->SetPoint(i,pt[i],(Float_t)pionNot[i+1]/(Float_t)pionCon[i+1]);
496     conPi->SetPointError(i,0,conErr);
497
498     eff = kaonEff[i+1]/TMath::Max(kaonTot[i+1],1);
499     effErr = TMath::Sqrt(eff*(1.-eff)/kaonTot[i+1]);
500     con = kaonNot[i+1]/TMath::Max(kaonCon[i+1],1);
501     conErr = TMath::Sqrt(con*(1.-con)/kaonCon[i+1]);
502     effKa->SetPoint(i,pt[i],(Float_t)kaonEff[i+1]/TMath::Max(kaonTot[i+1],1));
503     effKa->SetPointError(i,0,effErr);
504     conKa->SetPoint(i,pt[i],(Float_t)kaonNot[i+1]/TMath::Max(kaonCon[i+1],1));
505     conKa->SetPointError(i,0,conErr);
506   }
507
508   fProtGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtGen"));
509   fPbarGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarGen"));
510   fProtHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtHmp"));
511   fPbarHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarHmp"));
512
513   Int_t   *protGen = fProtGen->GetArray();
514   Int_t   *pbarGen = fPbarGen->GetArray();
515   Int_t   *protHmp = fProtHmp->GetArray();
516   Int_t   *pbarHmp = fPbarHmp->GetArray();
517
518   TGraphErrors *ratPr = new TGraphErrors(fN2);
519   TGraphErrors *ratPb = new TGraphErrors(fN2);
520
521   Float_t rat1=0, rat1Err=0, rat2=0, rat2Err=0;
522   for (Int_t i=0; i< fN2; i++){
523     rat1 = (Float_t)protHmp[i+1]/TMath::Max(protGen[i+1],1);
524     rat1Err = TMath::Sqrt(rat1*(1.-rat1)/TMath::Max(protGen[i+1],1));
525     rat2 = (Float_t)pbarHmp[i+1]/TMath::Max(pbarGen[i+1],1);
526     rat2Err = TMath::Sqrt(rat2*(1.-rat2)/TMath::Max(pbarGen[i+1],1));
527     ratPr->SetPoint(i,pt[i],rat1);
528     ratPr->SetPointError(i,0,rat1Err);
529     ratPb->SetPoint(i,pt[i],rat2);
530     ratPb->SetPointError(i,0,rat2Err);
531   }
532
533   TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
534   pCan->Divide(2,3);
535
536   pCan->cd(1);
537   effPi->SetTitle("Pions");
538   effPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
539   effPi->GetYaxis()->SetRangeUser(0.,1.);
540   effPi->SetMarkerStyle(20);
541   effPi->Draw("ALP");
542   conPi->SetMarkerStyle(21);
543   conPi->Draw("sameLP");
544
545   pCan->cd(3);
546   effKa->SetTitle("Kaons");
547   effKa->GetXaxis()->SetTitle("p_{T} (GeV/c)");
548   effKa->GetYaxis()->SetRangeUser(0.,1.);
549   effKa->SetMarkerStyle(20);
550   effKa->Draw("ALP");
551   conKa->SetMarkerStyle(21);
552   conKa->Draw("sameLP");
553
554   pCan->cd(5);
555   effPr->SetTitle("Protons");
556   effPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
557   effPr->GetYaxis()->SetRangeUser(0.,1.);
558   effPr->SetMarkerStyle(20);
559   effPr->Draw("ALP");
560   conPr->SetMarkerStyle(21);
561   conPr->Draw("sameLP");
562
563   pCan->cd(2);
564   ratPr->SetTitle("Ratios");
565   ratPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
566   ratPr->GetYaxis()->SetRangeUser(0.,1.);
567   ratPr->SetMarkerStyle(20);
568   ratPr->Draw("ALP");
569   ratPb->SetMarkerStyle(21);
570   ratPb->Draw("sameLP");
571
572   TFile *outFile = new TFile("HmpidGraphs.root","recreate");
573   pCan->Write();
574   outFile->Close();
575
576   AliAnalysisTaskSE::Terminate();
577
578 }
579
580 //___________________________________________________________________________
581 void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
582   //
583   //HERE ONE CAN CREATE OUTPUT OBJECTS
584   //TO BE SET BEFORE THE EXECUTION OF THE TASK
585   //
586
587   //slot #1
588 //   OpenFile(1);
589    fHmpHistList = new TList();
590    fHmpHistList->SetOwner();
591
592    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
593    fHmpHistList->Add(fHmpPesdPhmp);
594
595    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
596    fHmpHistList->Add(fHmpCkovPesd);
597
598    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
599    fHmpHistList->Add(fHmpCkovPhmp);
600
601    fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",400,0,20);
602    fHmpHistList->Add(fHmpMipTrkDist);
603
604    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
605    fHmpHistList->Add(fHmpMipTrkDistX);
606
607    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
608    fHmpHistList->Add(fHmpMipTrkDistY);
609
610    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
611    fHmpHistList->Add(fHmpMipCharge3cm);
612
613    fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
614    fHmpHistList->Add(fHmpMipCharge1cm);
615
616    fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
617    fHmpHistList->Add(fHmpNumPhots);
618
619    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
620    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
621    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
622    fHmpHistList->Add(fHmpTrkFlags);
623
624    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
625    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
626    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
627    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
628    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
629    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
630    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
631    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
632    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
633    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
634    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
635    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
636
637    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
638    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
639    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
640    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
641
642    fProtGen = new TH1I("ProtGen","Generated protons",fN2,0,fN2);
643    fPbarGen = new TH1I("PbarGen","Generated antiprotons",fN2,0,fN2);
644    fProtHmp = new TH1I("ProtHmp","HMPID protons",fN2,0,fN2);
645    fPbarHmp = new TH1I("PbarHmp","HMPID antiprotons",fN2,0,fN2);
646
647    fHmpHistList->Add(fProtGen); fHmpHistList->Add(fPbarGen);
648    fHmpHistList->Add(fProtHmp); fHmpHistList->Add(fPbarHmp);
649
650    fThetavsPiFromK= new TH2F("ThetavsPiFromK","Theta vs p of pions from K;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
651    fHmpHistList->Add(fThetavsPiFromK);
652
653    fThetapivsPesd = new TH2F("ThetapivsPesd","Theta of pions vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
654    fHmpHistList->Add(fThetapivsPesd);
655
656    fThetaKvsPesd  = new TH2F("ThetaKvsPesd","Theta of kaons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
657    fHmpHistList->Add(fThetaKvsPesd);
658
659    fThetaPvsPesd  = new TH2F("ThetaPvsPesd","Theta of protons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
660    fHmpHistList->Add(fThetaPvsPesd);
661
662 //   OpenFile(2);
663    fTree = new TTree("Tree","Tree with data");
664    fTree->Branch("Chamber",&fVar[0]);
665    fTree->Branch("pHmp3",&fVar[1]);
666    fTree->Branch("P",&fVar[2]);
667    fTree->Branch("Xpc",&fVar[3]);
668    fTree->Branch("Ypc",&fVar[4]);
669    fTree->Branch("X",&fVar[5]);
670    fTree->Branch("Y",&fVar[6]);
671    fTree->Branch("HMPIDsignal",&fVar[7]);
672    fTree->Branch("Charge",&fVar[8]);
673    fTree->Branch("Theta",&fVar[9]);
674    fTree->Branch("Phi",&fVar[10]);
675    fTree->Branch("Sign",&fVar[11]);
676    fTree->Branch("NumPhotons",&fVar[12]);
677    fTree->Branch("NumTPCclust",&fVar[13]);
678    fTree->Branch("Prob0",&fVar[14]);
679    fTree->Branch("Prob1",&fVar[15]);
680    fTree->Branch("Prob2",&fVar[16]);
681    fTree->Branch("Prob3",&fVar[17]);
682    fTree->Branch("Prob4",&fVar[18]);
683    fTree->Branch("TOFsignal",&fVar[19]);
684    fTree->Branch("KinkIndex",&fVar[20]);
685    fTree->Branch("Xv",&fVar[21]);
686    fTree->Branch("Yv",&fVar[22]);
687    fTree->Branch("Zv",&fVar[23]);
688    fTree->Branch("TPCchi2",&fVar[24]);
689    fTree->Branch("b0",&fVar[25]);
690    fTree->Branch("b1",&fVar[26]);
691    fTree->Branch("ClustSize",&fVar[27]);
692
693    PostData(1,fHmpHistList);
694    PostData(2,fTree);
695 }
696
697 //____________________________________________________________________________________________________________________________________
698 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
699 {
700  return abs(x - y) <= tolerance ;
701 }
702    
703 #endif