Move HMPID QA code to PWG1.
[u/mrichter/AliRoot.git] / PWG1 / HMPID / AliHMPIDTaskQA.cxx
CommitLineData
3882540a 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
39ClassImp(AliHMPIDTaskQA)
40
41//__________________________________________________________________________
42AliHMPIDTaskQA::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//___________________________________________________________________________
70AliHMPIDTaskQA::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//___________________________________________________________________________
103AliHMPIDTaskQA& 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//___________________________________________________________________________
139AliHMPIDTaskQA::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//___________________________________________________________________________
169AliHMPIDTaskQA::~AliHMPIDTaskQA() {
170 //
171 //destructor
172 //
173 Info("~AliHMPIDTaskQA","Calling Destructor");
174 if (fHmpHistList /*&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()*/) delete fHmpHistList;
175}
176
177//___________________________________________________________________________
178void 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//___________________________________________________________________________
201void 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//___________________________________________________________________________
342void 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//___________________________________________________________________________
354void 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//____________________________________________________________________________________________________________________________________
429Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
430{
431 return abs(x - y) <= tolerance ;
432}
433
434#endif