1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //==============================================================================
17 // AliHMPIDTaskQA - Class representing a quality check tool of HMPID
18 // A set of histograms is created.
19 //==============================================================================
21 #ifndef AliHMPIDTaskQA_CXX
22 #define AliHMPIDTaskQA_CXX
29 #include "TGraphErrors.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDInputHandler.h"
32 #include "AliMCEventHandler.h"
33 #include "AliMCEvent.h"
34 #include "AliESDtrack.h"
37 #include "AliHMPIDTaskQA.h"
39 ClassImp(AliHMPIDTaskQA)
41 //__________________________________________________________________________
42 AliHMPIDTaskQA::AliHMPIDTaskQA() :
43 fESD(0x0),fMC(0x0),fUseMC(kTRUE),
45 fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
46 fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
66 for (Int_t i=0; i<7; i++){
69 fHmpPhotSin2th[i] = 0x0;
70 fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
71 fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
72 fHmpMipCharge[i] = 0x0;
76 //___________________________________________________________________________
77 AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
78 AliAnalysisTaskSE(name),
79 fESD(0x0), fMC(0x0), fUseMC(kTRUE),
81 fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
82 fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
100 // Constructor. Initialization of Inputs and Outputs
102 for (Int_t i=0; i<7; i++){
103 fHmpPhotons[i] = 0x0;
105 fHmpPhotSin2th[i] = 0x0;
106 fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
107 fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
108 fHmpMipCharge[i] = 0x0;
111 DefineOutput(1,TList::Class());
114 //___________________________________________________________________________
115 AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
118 // Assignment operator
121 AliAnalysisTaskSE::operator=(c);
125 fHmpHistList = c.fHmpHistList;
126 fHmpPesdPhmp = c.fHmpPesdPhmp;
127 fHmpCkovPesd = c.fHmpCkovPesd;
128 fHmpCkovPhmp = c.fHmpCkovPhmp;
129 fHmpMipCharge3cm = c.fHmpMipCharge3cm;
130 fHmpMipTrkDist = c.fHmpMipTrkDist;
131 fHmpTrkFlags = c.fHmpTrkFlags;
134 fPionEff = c.fPionEff;
135 fKaonEff = c.fKaonEff;
136 fProtEff = c.fProtEff;
137 fPionTot = c.fPionTot;
138 fKaonTot = c.fKaonTot;
139 fProtTot = c.fProtTot;
140 fPionNot = c.fPionNot;
141 fKaonNot = c.fKaonNot;
142 fProtNot = c.fProtNot;
143 fPionCon = c.fPionCon;
144 fKaonCon = c.fKaonCon;
145 fProtCon = c.fProtCon;
146 for (Int_t i=0; i<7; i++){
147 fHmpPhotons[i] = c.fHmpPhotons[i];
148 fHmpPhotP[i] = c.fHmpPhotP[i];
149 fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
150 fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
151 fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
152 fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
153 fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
154 fHmpMipCharge[i] = c.fHmpMipCharge[i];
160 //___________________________________________________________________________
161 AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
162 AliAnalysisTaskSE(c),
163 fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
164 fHmpHistList(c.fHmpHistList),
165 fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
166 fHmpMipCharge3cm(c.fHmpMipCharge3cm),fHmpMipTrkDist(c.fHmpMipTrkDist),
167 fHmpTrkFlags(c.fHmpTrkFlags),
170 fPionEff(c.fPionEff),
171 fKaonEff(c.fKaonEff),
172 fProtEff(c.fProtEff),
173 fPionTot(c.fPionTot),
174 fKaonTot(c.fKaonTot),
175 fProtTot(c.fProtTot),
176 fPionNot(c.fPionNot),
177 fKaonNot(c.fKaonNot),
178 fProtNot(c.fProtNot),
179 fPionCon(c.fPionCon),
180 fKaonCon(c.fKaonCon),
186 for (Int_t i=0; i<7; i++){
187 fHmpPhotons[i] = c.fHmpPhotons[i];
188 fHmpPhotP[i] = c.fHmpPhotP[i];
189 fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
190 fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
191 fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
192 fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
193 fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
194 fHmpMipCharge[i] = c.fHmpMipCharge[i];
198 //___________________________________________________________________________
199 AliHMPIDTaskQA::~AliHMPIDTaskQA() {
203 Info("~AliHMPIDTaskQA","Calling Destructor");
204 if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
207 //___________________________________________________________________________
208 void AliHMPIDTaskQA::ConnectInputData(Option_t *option)
210 AliAnalysisTaskSE::ConnectInputData(option);
213 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
215 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
217 fESD = esdH->GetEvent();
221 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
223 AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
229 //___________________________________________________________________________
230 void AliHMPIDTaskQA::UserExec(Option_t *)
232 Double_t priors[10]={1.,1.,1.,1.,1.,0.,0.,0.,0.,0.}; //{0.01,0.01,0.83,0.10,0.5};
234 AliPID *pPid = new AliPID();
235 pPid->SetPriors(priors);
236 AliESDtrack *track=0;
238 AliStack* pStack = 0;
241 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
243 fMC = mcH->MCEvent();
244 pStack = fMC->Stack();
246 AliFatal("fMC flag set but MC handler not available");
252 // Main loop function, executed on Event basis
255 if( !fESD->GetPrimaryVertex()) return;
256 if( fESD->GetPrimaryVertex()->GetNContributors() < 1 ) return;
257 if( TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > 10.0 /* cm */) return;
260 for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) {
262 track = fESD->GetTrack(iTrack);
264 if(! track->IsOn(AliESDtrack::kITSrefit) || !track->IsOn(AliESDtrack::kITSrefit) ) continue;
266 Double_t rin[3], rout[3];
267 track->GetInnerXYZ(rin);
268 track->GetOuterXYZ(rout);
269 Double_t ktol = 0.001;
270 Double_t thetaCh = track->GetHMPIDsignal();
272 if(Equal(thetaCh,-20.,ktol)) fHmpTrkFlags->Fill(0);
273 else if(Equal(thetaCh,-9.,ktol)) fHmpTrkFlags->Fill(1);
274 else if(Equal(thetaCh,-5.,ktol)) fHmpTrkFlags->Fill(2);
275 else if(Equal(thetaCh,-11.,ktol)) fHmpTrkFlags->Fill(3);
276 else if(Equal(thetaCh,-22.,ktol)) fHmpTrkFlags->Fill(4);
277 else fHmpTrkFlags->Fill(4);
279 if(Equal(thetaCh,-20.,ktol)) continue;
280 if(track->GetHMPIDcluIdx() < 0) continue;
284 Float_t xpc, ypc, th, ph;
285 track->GetHMPIDmip(x,y,q,nph);
286 track->GetHMPIDtrk(xpc,ypc,th,ph);
288 if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
290 Float_t sign = track->GetSign();
291 Int_t ch = track->GetHMPIDcluIdx()/1000000;
292 if( ch < 0 || ch > 6 ) continue; //Chamber numbering is [0,6]
294 if (sign > 0.) fHmpMipTrkDistPosX[ch]->Fill(xpc - x), fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
295 if (sign < 0.) fHmpMipTrkDistNegX[ch]->Fill(xpc - x), fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
296 fHmpMipCharge[ch]->Fill(q);
297 Double_t dist = TMath::Sqrt((xpc - x)*(xpc - x) + (ypc - y)*(ypc - y));
298 if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
299 fHmpMipTrkDist->Fill(dist);
301 track->GetHMPIDpid(probs);
302 pPid->SetProbabilities(probs);
304 if ((label = track->GetLabel()) < 0) continue;
305 pPart = pStack->Particle(label);
309 Double_t pHmp[3] = {0}, pHmp3 = 0;
310 if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
311 if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
312 fHmpCkovPesd->Fill(track->P(),thetaCh);
313 if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,thetaCh);
314 fHmpPhotons[ch]->Fill(nph);
315 fHmpPhotP[ch]->Fill(pHmp3,nph);
316 Double_t sin2 = TMath::Power(TMath::Sin(th),2);
317 fHmpPhotSin2th[ch]->Fill(sin2,nph);
319 if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
320 if (!pStack->IsPhysicalPrimary(label)) continue;
321 if (track->Pt()<1. || track->Pt()>5.) continue;
322 Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
323 Int_t ptBin=(Int_t) (2*(track->Pt()-1));
324 if (pdgCode!=2212) fProtCon->Fill(ptBin);
326 fProtTot->Fill(ptBin);
327 fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
328 fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
329 fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
331 if (pdgCode!=211) fPionCon->Fill(ptBin);
332 if (pdgCode!=321) fKaonCon->Fill(ptBin);
335 Float_t weight=pPid->GetProbability(AliPID::kElectron)+
336 pPid->GetProbability(AliPID::kMuon)+
337 pPid->GetProbability(AliPID::kPion);
338 fPionTot->Fill(ptBin);
339 fPionEff->Fill(ptBin,weight);
340 fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
342 fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
346 fKaonTot->Fill(ptBin);
347 fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
348 fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
350 fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
357 /* PostData(0) is taken care of by AliAnalysisTaskSE */
358 PostData(1,fHmpHistList);
361 //___________________________________________________________________________
362 void AliHMPIDTaskQA::Terminate(Option_t*)
364 // The Terminate() function is the last function to be called during
365 // a query. It always runs on the client, it can be used to present
366 // the results graphically or save the results to file.
368 Info("Terminate"," ");
369 AliAnalysisTaskSE::Terminate();
373 //___________________________________________________________________________
374 void AliHMPIDTaskQA::UserCreateOutputObjects() {
376 //HERE ONE CAN CREATE OUTPUT OBJECTS
377 //TO BE SET BEFORE THE EXECUTION OF THE TASK
382 fHmpHistList = new TList();
383 fHmpHistList->SetOwner();
385 fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
386 fHmpHistList->Add(fHmpPesdPhmp);
388 fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",500,0,10,500,0,1.0);
389 fHmpHistList->Add(fHmpCkovPesd);
391 fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",500,0,10,500,0,1.0);
392 fHmpHistList->Add(fHmpCkovPhmp);
394 for (Int_t i=0; i<7; i++){
395 TString title=Form("MIP-Track distance in local X, Chamber %d;distance (cm);Entries",i);
396 fHmpMipTrkDistPosX[i] = new TH1F(Form("fHmpMipTrkDistPosX%d",i),title.Data(),800,-20,20);
397 fHmpHistList->Add(fHmpMipTrkDistPosX[i]);
398 fHmpMipTrkDistNegX[i] = new TH1F(Form("fHmpMipTrkDistNegX%d",i),title.Data(),800,-20,20);
399 fHmpHistList->Add(fHmpMipTrkDistNegX[i]);
401 title=Form("MIP-Track distance in local Y, Chamber %d;distance (cm);Entries",i);
402 fHmpMipTrkDistPosY[i] = new TH1F(Form("fHmpMipTrkDistPosY%d",i),title.Data(),800,-20,20);
403 fHmpHistList->Add(fHmpMipTrkDistPosY[i]);
404 fHmpMipTrkDistNegY[i] = new TH1F(Form("fHmpMipTrkDistNegY%d",i),title.Data(),800,-20,20);
405 fHmpHistList->Add(fHmpMipTrkDistNegY[i]);
407 title=Form("Mip charge distribution, Chamber %d;Charge;Entries",i);
408 fHmpMipCharge[i] = new TH1F(Form("fHmpMipCharge%d",i),title.Data(),5001,-0.5,5000.5);
409 fHmpHistList->Add(fHmpMipCharge[i]);
411 title=Form("Photons, Chamber %d;N of Photons;Entries",i);
412 fHmpPhotons[i] = new TH1F(Form("fHmpPhotons%d",i),title.Data(),100,0,100);
413 fHmpHistList->Add(fHmpPhotons[i]);
415 title=Form("Photons versus HMP momentum, Chamber %d;P (GeV);N of Photons",i);
416 fHmpPhotP[i] = new TH2F(Form("fHmpPhotP%d",i),title.Data(),200,0.,20.,100,0,100);
417 fHmpHistList->Add(fHmpPhotP[i]);
419 title=Form("Photons versus sin(th)^2 (uncorrected), Chamber %d;P (GeV);N of Photons",i);
420 fHmpPhotSin2th[i] = new TH2F(Form("fHmpPhotSin2th%d",i),title.Data(),100,0.,1.,100,0,100);
421 fHmpHistList->Add(fHmpPhotSin2th[i]);
424 fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
425 fHmpHistList->Add(fHmpMipCharge3cm);
427 fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","Mip-track distance in all the chambers",100,0.,10.);
428 fHmpHistList->Add(fHmpMipTrkDist);
430 fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
431 TString summary[6] = {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
432 for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i %s",ibin+1,summary[ibin].Data()));
433 fHmpHistList->Add(fHmpTrkFlags);
435 fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
436 fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
437 fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
438 fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
439 fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
440 fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
441 fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
442 fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
443 fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
444 fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
445 fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
446 fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
448 fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
449 fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
450 fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
451 fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
453 PostData(1,fHmpHistList);
456 //____________________________________________________________________________________________________________________________________
457 Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
459 return abs(x - y) <= tolerance ;