Fix Coverity defect
[u/mrichter/AliRoot.git] / PWG1 / HMPID / AliHMPIDTaskQA.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 // AliHMPIDTaskQA - Class representing a quality check tool of HMPID 
18 // A set of histograms is created.
19 //==============================================================================
20
21 #ifndef AliHMPIDTaskQA_CXX
22 #define AliHMPIDTaskQA_CXX
23
24
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TFile.h"
28 #include "TCanvas.h"
29 #include "TGraphErrors.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDInputHandler.h"
32 #include "AliMCEventHandler.h"
33 #include "AliMCEvent.h"
34 #include "AliESDtrack.h"
35 #include "AliPID.h"
36 #include "AliLog.h"
37 #include "AliHMPIDTaskQA.h"
38
39 ClassImp(AliHMPIDTaskQA)
40
41 //__________________________________________________________________________
42 AliHMPIDTaskQA::AliHMPIDTaskQA() :
43   fESD(0x0),fMC(0x0),fUseMC(kTRUE),
44   fHmpHistList(0x0),
45   fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),fHmpMipCharge3cm(0x0),
46   fHmpTrkFlags(0x0),
47   fN1(6),
48   fN2(8),
49   fPionEff(0x0),
50   fKaonEff(0x0),
51   fProtEff(0x0),
52   fPionTot(0x0),
53   fKaonTot(0x0),
54   fProtTot(0x0),
55   fPionNot(0x0),
56   fKaonNot(0x0),
57   fProtNot(0x0),
58   fPionCon(0x0),
59   fKaonCon(0x0),
60   fProtCon(0x0),
61   fTree(0x0)
62 {
63   //
64   //Default ctor
65   //
66   for (Int_t i=0; i<28; i++) fVar[i]=0;
67 }
68
69 //___________________________________________________________________________
70 AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
71   AliAnalysisTaskSE(name),
72   fESD(0x0), fMC(0x0), fUseMC(kTRUE),
73   fHmpHistList(0x0),
74   fHmpMipTrkDistX(0x0), fHmpMipTrkDistY(0x0),
75   fHmpMipCharge3cm(0x0),
76   fHmpTrkFlags(0x0),
77   fN1(6),
78   fN2(8),
79   fPionEff(0x0),
80   fKaonEff(0x0),
81   fProtEff(0x0),
82   fPionTot(0x0),
83   fKaonTot(0x0),
84   fProtTot(0x0),
85   fPionNot(0x0),
86   fKaonNot(0x0),
87   fProtNot(0x0),
88   fPionCon(0x0),
89   fKaonCon(0x0),
90   fProtCon(0x0),
91   fTree(0x0)
92 {
93   //
94   // Constructor. Initialization of Inputs and Outputs
95   //
96   for (Int_t i=0; i<28; i++) fVar[i]=0;
97
98   DefineOutput(1,TList::Class());
99   DefineOutput(2,TTree::Class());
100 }
101
102 //___________________________________________________________________________
103 AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
104 {
105   //
106   // Assignment operator
107   //
108   if (this!=&c) {
109     AliAnalysisTaskSE::operator=(c);
110     fESD             = c.fESD;
111     fMC              = c.fMC;
112     fUseMC           = c.fUseMC;
113     fHmpHistList     = c.fHmpHistList;
114     fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
115     fHmpMipTrkDistY  = c.fHmpMipTrkDistY;
116     fHmpMipCharge3cm = c.fHmpMipCharge3cm;
117     fHmpTrkFlags     = c.fHmpTrkFlags;
118     fN1              = c.fN1;
119     fN2              = c.fN2;
120     fPionEff         = c.fPionEff;
121     fKaonEff         = c.fKaonEff;
122     fProtEff         = c.fProtEff;
123     fPionTot         = c.fPionTot;
124     fKaonTot         = c.fKaonTot;
125     fProtTot         = c.fProtTot;
126     fPionNot         = c.fPionNot;
127     fKaonNot         = c.fKaonNot;
128     fProtNot         = c.fProtNot;
129     fPionCon         = c.fPionCon;
130     fKaonCon         = c.fKaonCon;
131     fProtCon         = c.fProtCon;
132     fTree            = c.fTree;
133     for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
134   }
135   return *this;
136 }
137
138 //___________________________________________________________________________
139 AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
140   AliAnalysisTaskSE(c),
141   fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
142   fHmpHistList(c.fHmpHistList),
143   fHmpMipTrkDistX(c.fHmpMipTrkDistX),fHmpMipTrkDistY(c.fHmpMipTrkDistY),
144   fHmpMipCharge3cm(c.fHmpMipCharge3cm),
145   fHmpTrkFlags(c.fHmpTrkFlags),
146   fN1(c.fN1),
147   fN2(c.fN2),
148   fPionEff(c.fPionEff),
149   fKaonEff(c.fKaonEff),
150   fProtEff(c.fProtEff),
151   fPionTot(c.fPionTot),
152   fKaonTot(c.fKaonTot),
153   fProtTot(c.fProtTot),
154   fPionNot(c.fPionNot),
155   fKaonNot(c.fKaonNot),
156   fProtNot(c.fProtNot),
157   fPionCon(c.fPionCon),
158   fKaonCon(c.fKaonCon),
159   fProtCon(c.fProtCon),
160   fTree(c.fTree)
161 {
162   //
163   // Copy Constructor
164   //
165   for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
166 }
167  
168 //___________________________________________________________________________
169 AliHMPIDTaskQA::~AliHMPIDTaskQA() {
170   //
171   //destructor
172   //
173   Info("~AliHMPIDTaskQA","Calling Destructor");
174   if (fHmpHistList /*&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()*/) delete fHmpHistList;
175 }
176
177 //___________________________________________________________________________
178 void AliHMPIDTaskQA::ConnectInputData(Option_t *)
179 {
180   // Connect ESD here
181
182   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
183   if (!esdH) {
184     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
185   } else
186     fESD = esdH->GetEvent();
187
188   if (fUseMC){
189     // Connect MC
190     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
191     if (!mcH) {
192       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
193       fUseMC = kFALSE;
194     } else
195       fMC = mcH->MCEvent();
196       if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
197   }
198 }
199
200 //___________________________________________________________________________
201 void AliHMPIDTaskQA::UserExec(Option_t *)
202 {
203   Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
204   Double_t probs[5];
205   AliPID *pPid = new AliPID();
206   pPid->SetPriors(priors);
207   Double_t n = 1.293;
208   Double_t dGeVMass[] = {0.000511,0.105658,0.13957018,0.493677,0.938272};
209   AliESDtrack *track=0;
210   TParticle *pPart=0;
211   AliStack* pStack = 0;
212   Int_t label;
213   if (fUseMC){
214     pStack = fMC->Stack();
215   }
216
217   //
218   // Main loop function, executed on Event basis
219   //
220   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
221
222     track = fESD->GetTrack(iTrack);
223     if(!track) continue;
224     Double_t dPtr = track->P();
225     Double_t rin[3], rout[3];
226     track->GetInnerXYZ(rin);
227     track->GetOuterXYZ(rout);
228     Double_t ktol = 0.001;
229
230     if(Equal(track->GetHMPIDsignal(),-20.,ktol))      fHmpTrkFlags->Fill(0);
231     else if(Equal(track->GetHMPIDsignal(),-9.,ktol))  fHmpTrkFlags->Fill(1);
232     else if(Equal(track->GetHMPIDsignal(),-5.,ktol))  fHmpTrkFlags->Fill(2);
233     else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
234     else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
235     else fHmpTrkFlags->Fill(4);
236
237     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
238     if(track->GetHMPIDcluIdx() < 0) continue;
239
240     Int_t q, nph;
241     Float_t x, y;
242     Float_t xpc, ypc, th, ph;
243     track->GetHMPIDmip(x,y,q,nph);
244     track->GetHMPIDtrk(xpc,ypc,th,ph);
245
246     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
247
248     Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
249     if(q>100.) fHmpMipTrkDistX->Fill(xpc - x);
250     if(q>100.) fHmpMipTrkDistY->Fill(ypc - y);
251     Double_t pHmp[3] = {0}, pHmp3 = 0;
252     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
253     if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
254
255     Float_t b[2];
256     Float_t bCov[3];
257     track->GetImpactParameters(b,bCov);    
258
259     track->GetHMPIDpid(probs);
260     pPid->SetProbabilities(probs);
261     if (fUseMC){
262       if ((label = track->GetLabel()) < 0) continue;
263       pPart = pStack->Particle(label);
264     }
265
266     if(track->GetHMPIDsignal() > 0 ){
267
268       if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
269         if (!pStack->IsPhysicalPrimary(label)) continue;
270         if (track->Pt()<1. || track->Pt()>5.) continue;
271         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
272         Int_t ptBin=(Int_t) (2*(track->Pt()-1));
273         if (pdgCode!=2212) fProtCon->Fill(ptBin);
274         if (pdgCode==2212){
275           fProtTot->Fill(ptBin);
276           fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
277           fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
278           fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
279         }
280         if (pdgCode!=211) fPionCon->Fill(ptBin);
281         if (pdgCode!=321) fKaonCon->Fill(ptBin);
282         if (pdgCode==211){
283           if (ptBin < 6){
284             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
285                            pPid->GetProbability(AliPID::kMuon)+
286                            pPid->GetProbability(AliPID::kPion);
287             fPionTot->Fill(ptBin);
288             fPionEff->Fill(ptBin,weight);
289             fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
290           }
291           fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
292         }
293         if (pdgCode==321){
294           if (ptBin < 6){
295             fKaonTot->Fill(ptBin);
296             fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
297             fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
298           }
299           fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
300         }
301       }
302     }//there is signal
303
304     fVar[0] = track->GetHMPIDcluIdx()/1000000;
305     fVar[1] = pHmp3;
306     fVar[2] = dPtr;
307     fVar[3] = xpc;
308     fVar[4] = ypc;
309     fVar[5] = x;
310     fVar[6] = y;
311     fVar[7] = (Float_t)track->GetHMPIDsignal();
312     fVar[8] = q;
313     fVar[9] = th;
314     fVar[10] = ph;
315     fVar[11] = (Float_t)track->GetSign();
316     fVar[12] = (Float_t)nph;
317     fVar[13] = (Float_t)track->GetNcls(1);
318     fVar[14] = (Float_t)probs[0];
319     fVar[15] = (Float_t)probs[1];
320     fVar[16] = (Float_t)probs[2];
321     fVar[17] = (Float_t)probs[3];
322     fVar[18] = (Float_t)probs[4];
323     fVar[19] = (Float_t)track->GetTOFsignal();
324     fVar[20] = (Float_t)track->GetKinkIndex(0);
325     fVar[21] = (Float_t)track->Xv();
326     fVar[22] = (Float_t)track->Yv();
327     fVar[23] = (Float_t)track->Zv();
328     fVar[24] = (Float_t)track->GetTPCchi2();
329     fVar[25] = b[0];
330     fVar[26] = b[1];
331     fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
332     fTree->Fill();
333   }//track loop
334   delete pPid;
335
336   /* PostData(0) is taken care of by AliAnalysisTaskSE */
337   PostData(1,fHmpHistList);
338   PostData(2,fTree);
339 }
340
341 //___________________________________________________________________________
342 void AliHMPIDTaskQA::Terminate(Option_t*)
343 {
344   // The Terminate() function is the last function to be called during
345   // a query. It always runs on the client, it can be used to present
346   // the results graphically or save the results to file.
347
348   Info("Terminate"," ");
349   AliAnalysisTaskSE::Terminate();
350
351 }
352
353 //___________________________________________________________________________
354 void AliHMPIDTaskQA::UserCreateOutputObjects() {
355   //
356   //HERE ONE CAN CREATE OUTPUT OBJECTS
357   //TO BE SET BEFORE THE EXECUTION OF THE TASK
358   //
359
360   //slot #1
361    OpenFile(1);
362    fHmpHistList = new TList();
363
364    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
365    fHmpHistList->Add(fHmpMipTrkDistX);
366
367    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
368    fHmpHistList->Add(fHmpMipTrkDistY);
369
370    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
371    fHmpHistList->Add(fHmpMipCharge3cm);
372
373    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
374    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
375    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
376    fHmpHistList->Add(fHmpTrkFlags);
377
378    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
379    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
380    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
381    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
382    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
383    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
384    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
385    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
386    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
387    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
388    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
389    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
390
391    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
392    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
393    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
394    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
395
396    OpenFile(2);
397    fTree = new TTree("Tree","Tree with data");
398    fTree->Branch("Chamber",&fVar[0]);
399    fTree->Branch("pHmp3",&fVar[1]);
400    fTree->Branch("P",&fVar[2]);
401    fTree->Branch("Xpc",&fVar[3]);
402    fTree->Branch("Ypc",&fVar[4]);
403    fTree->Branch("X",&fVar[5]);
404    fTree->Branch("Y",&fVar[6]);
405    fTree->Branch("HMPIDsignal",&fVar[7]);
406    fTree->Branch("Charge",&fVar[8]);
407    fTree->Branch("Theta",&fVar[9]);
408    fTree->Branch("Phi",&fVar[10]);
409    fTree->Branch("Sign",&fVar[11]);
410    fTree->Branch("NumPhotons",&fVar[12]);
411    fTree->Branch("NumTPCclust",&fVar[13]);
412    fTree->Branch("Prob0",&fVar[14]);
413    fTree->Branch("Prob1",&fVar[15]);
414    fTree->Branch("Prob2",&fVar[16]);
415    fTree->Branch("Prob3",&fVar[17]);
416    fTree->Branch("Prob4",&fVar[18]);
417    fTree->Branch("TOFsignal",&fVar[19]);
418    fTree->Branch("KinkIndex",&fVar[20]);
419    fTree->Branch("Xv",&fVar[21]);
420    fTree->Branch("Yv",&fVar[22]);
421    fTree->Branch("Zv",&fVar[23]);
422    fTree->Branch("TPCchi2",&fVar[24]);
423    fTree->Branch("b0",&fVar[25]);
424    fTree->Branch("b1",&fVar[26]);
425    fTree->Branch("ClustSize",&fVar[27]);
426 }
427
428 //____________________________________________________________________________________________________________________________________
429 Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
430 {
431  return abs(x - y) <= tolerance ;
432 }
433    
434 #endif