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