]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HMPID/AliHMPIDAnalysisTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
addapt array of histograms for different SM combinations to new EMCAL SM
[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),
dcb23f7e 49 fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
e7c33863 50 fHmpMipTrkDist(0x0),fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),
51 fHmpMipCharge3cm(0x0),fHmpMipCharge1cm(0x0),
52 fHmpNumPhots(0x0),fHmpTrkFlags(0x0),
740b6bde 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),
f5dd9ecb 70 fThetaPvsPesd(0x0),
e7c33863 71 fProtGen(0x0),
72 fPbarGen(0x0),
73 fProtHmp(0x0),
74 fPbarHmp(0x0),
f5dd9ecb 75 fTree(0x0)
9a8aafa5 76{
77 //
78 //Default ctor
79 //
e7c33863 80 for (Int_t i=0; i<28; i++) fVar[i]=0;
9a8aafa5 81}
740b6bde 82
9a8aafa5 83//___________________________________________________________________________
84AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
85 AliAnalysisTaskSE(name),
66cdae53 86 fESD(0x0), fMC(0x0), fUseMC(kTRUE),
e7c33863 87 fHmpHistList(0x0),
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),
e7c33863 110 fProtGen(0x0),
111 fPbarGen(0x0),
112 fProtHmp(0x0),
113 fPbarHmp(0x0),
f5dd9ecb 114 fTree(0x0)
9a8aafa5 115{
116 //
117 // Constructor. Initialization of Inputs and Outputs
118 //
e7c33863 119 for (Int_t i=0; i<28; i++) fVar[i]=0;
740b6bde 120
66cdae53 121 DefineOutput(1,TList::Class());
f5dd9ecb 122 DefineOutput(2,TTree::Class());
740b6bde 123}
124
dcb23f7e 125//___________________________________________________________________________
740b6bde 126AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
dcb23f7e 127{
128 //
129 // Assignment operator
130 //
131 if (this!=&c) {
740b6bde 132 AliAnalysisTaskSE::operator=(c);
133 fESD = c.fESD;
66cdae53 134 fMC = c.fMC;
135 fUseMC = c.fUseMC;
740b6bde 136 fHmpHistList = c.fHmpHistList;
740b6bde 137 fHmpPesdPhmp = c.fHmpPesdPhmp;
138 fHmpCkovPesd = c.fHmpCkovPesd;
139 fHmpCkovPhmp = c.fHmpCkovPhmp;
dcb23f7e 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;
740b6bde 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;
e7c33863 165 fProtGen = c.fProtGen;
166 fPbarGen = c.fPbarGen;
167 fProtHmp = c.fProtHmp;
168 fPbarHmp = c.fPbarHmp;
f5dd9ecb 169 fTree = c.fTree;
e7c33863 170 for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
dcb23f7e 171 }
172 return *this;
740b6bde 173}
dcb23f7e 174
175//___________________________________________________________________________
176AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
177 AliAnalysisTaskSE(c),
66cdae53 178 fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
e7c33863 179 fHmpHistList(c.fHmpHistList),
dcb23f7e 180 fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
740b6bde 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),
f5dd9ecb 201 fThetaPvsPesd(c.fThetaPvsPesd),
e7c33863 202 fProtGen(c.fProtGen),
203 fPbarGen(c.fPbarGen),
204 fProtHmp(c.fProtHmp),
205 fPbarHmp(c.fPbarHmp),
f5dd9ecb 206 fTree(c.fTree)
dcb23f7e 207{
208 //
209 // Copy Constructor
210 //
e7c33863 211 for (Int_t i=0; i<28; i++) fVar[i]=c.fVar[i];
dcb23f7e 212}
213
9a8aafa5 214//___________________________________________________________________________
215AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
216 //
217 //destructor
218 //
740b6bde 219 Info("~AliHMPIDAnalysisTask","Calling Destructor");
e7c33863 220 if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
9a8aafa5 221}
740b6bde 222
38b8c336 223//___________________________________________________________________________
f5dd9ecb 224void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
38b8c336 225{
226 // Connect ESD here
740b6bde 227
38b8c336 228 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
229 if (!esdH) {
66cdae53 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) {
740b6bde 238 AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
66cdae53 239 fUseMC = kFALSE;
740b6bde 240 } else
241 fMC = mcH->MCEvent();
f5dd9ecb 242 if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
66cdae53 243 }
740b6bde 244}
245
66cdae53 246//___________________________________________________________________________
247void AliHMPIDAnalysisTask::UserExec(Option_t *)
9a8aafa5 248{
740b6bde 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);
e7c33863 253 AliESDtrack *track = 0;
254 TParticle *pPart = 0;
66cdae53 255 AliStack* pStack = 0;
e7c33863 256 Int_t label = -1;
66cdae53 257 if (fUseMC){
258 pStack = fMC->Stack();
259 }
260
e7c33863 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
9a8aafa5 266 //
267 // Main loop function, executed on Event basis
268 //
9a8aafa5 269 for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
270
740b6bde 271 track = fESD->GetTrack(iTrack);
9a8aafa5 272 if(!track) continue;
9a8aafa5 273 track->GetInnerXYZ(rin);
274 track->GetOuterXYZ(rout);
740b6bde 275
276 if(Equal(track->GetHMPIDsignal(),-20.,ktol)) fHmpTrkFlags->Fill(0);
277 else if(Equal(track->GetHMPIDsignal(),-9.,ktol)) fHmpTrkFlags->Fill(1);
66cdae53 278 else if(Equal(track->GetHMPIDsignal(),-5.,ktol)) fHmpTrkFlags->Fill(2);
279 else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
740b6bde 280 else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
e7c33863 281 else fHmpTrkFlags->Fill(5);
740b6bde 282
dcb23f7e 283 if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
38b8c336 284 if(track->GetHMPIDcluIdx() < 0) continue;
740b6bde 285
740b6bde 286 Int_t q, nph;
287 Float_t x, y;
9a8aafa5 288 Float_t xpc, ypc, th, ph;
38b8c336 289 track->GetHMPIDmip(x,y,q,nph);
290 track->GetHMPIDtrk(xpc,ypc,th,ph);
740b6bde 291
dcb23f7e 292 if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
740b6bde 293
e7c33863 294 Float_t dist = TMath::Sqrt((xpc-x)*(xpc-x)+(ypc-y)*(ypc-y));
38b8c336 295 fHmpMipTrkDist->Fill(dist);
740b6bde 296 fHmpMipTrkDistX->Fill(xpc - x);
38b8c336 297 fHmpMipTrkDistY->Fill(ypc - y);
740b6bde 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
f5dd9ecb 302 Float_t b[2];
303 Float_t bCov[3];
304 track->GetImpactParameters(b,bCov);
305
740b6bde 306 track->GetHMPIDpid(probs);
307 pPid->SetProbabilities(probs);
66cdae53 308 if (fUseMC){
309 if ((label = track->GetLabel()) < 0) continue;
310 pPart = pStack->Particle(label);
311 }
740b6bde 312
66cdae53 313 if(track->GetHMPIDsignal() > 0 ){
f5dd9ecb 314 if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
740b6bde 315 if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
38b8c336 316 fHmpNumPhots->Fill(nph);
317 fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
f5dd9ecb 318 if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
740b6bde 319
e7c33863 320 if (fUseMC && dist<distCut && TMath::Abs(th)<thetaCut){
66cdae53 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 }
740b6bde 331 }
66cdae53 332 if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),track->GetHMPIDsignal());
333 if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),track->GetHMPIDsignal());
334
e7c33863 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);
66cdae53 338 if (pdgCode==2212){
e7c33863 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));
740b6bde 343 }
e7c33863 344 if (pdgCode!=211) fPionCon->Fill(pBin);
345 if (pdgCode!=321) fKaonCon->Fill(pBin);
66cdae53 346 if (pdgCode==211){
e7c33863 347 if (pBin < 6){
66cdae53 348 Float_t weight=pPid->GetProbability(AliPID::kElectron)+
349 pPid->GetProbability(AliPID::kMuon)+
350 pPid->GetProbability(AliPID::kPion);
e7c33863 351 fPionTot->Fill(pBin);
352 fPionEff->Fill(pBin,weight);
353 fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
66cdae53 354 }
e7c33863 355 fProtNot->Fill(pBin,pPid->GetProbability(AliPID::kProton));
66cdae53 356 }
357 if (pdgCode==321){
e7c33863 358 if (pBin < 6){
359 fKaonTot->Fill(pBin);
360 fKaonEff->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
361 fPionNot->Fill(pBin,(pPid->GetProbability(AliPID::kPion)));
66cdae53 362 }
e7c33863 363 fProtNot->Fill(pBin,(pPid->GetProbability(AliPID::kProton)));
740b6bde 364 }
740b6bde 365 }
38b8c336 366 }//there is signal
f5dd9ecb 367 fVar[0] = track->GetHMPIDcluIdx()/1000000;
368 fVar[1] = pHmp3;
e7c33863 369 fVar[2] = (Float_t)track->P();
f5dd9ecb 370 fVar[3] = xpc;
371 fVar[4] = ypc;
372 fVar[5] = x;
373 fVar[6] = y;
e7c33863 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;
f5dd9ecb 395 fTree->Fill();
740b6bde 396 }//track loop
740b6bde 397 delete pPid;
398
e7c33863 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
9a8aafa5 420 /* PostData(0) is taken care of by AliAnalysisTaskSE */
66cdae53 421 PostData(1,fHmpHistList);
f5dd9ecb 422 PostData(2,fTree);
9a8aafa5 423}
424
9a8aafa5 425//___________________________________________________________________________
426void 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
e7c33863 432 Info("Terminate"," ");
66cdae53 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 }
740b6bde 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
e7c33863 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();
740b6bde 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
e7c33863 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
740b6bde 533 TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
e7c33863 534 pCan->Divide(2,3);
740b6bde 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
e7c33863 545 pCan->cd(3);
740b6bde 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
e7c33863 554 pCan->cd(5);
740b6bde 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
e7c33863 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
740b6bde 572 TFile *outFile = new TFile("HmpidGraphs.root","recreate");
573 pCan->Write();
574 outFile->Close();
575
9a8aafa5 576 AliAnalysisTaskSE::Terminate();
577
578}
579
9a8aafa5 580//___________________________________________________________________________
66cdae53 581void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
9a8aafa5 582 //
583 //HERE ONE CAN CREATE OUTPUT OBJECTS
584 //TO BE SET BEFORE THE EXECUTION OF THE TASK
585 //
9a8aafa5 586
587 //slot #1
e7c33863 588// OpenFile(1);
38b8c336 589 fHmpHistList = new TList();
e7c33863 590 fHmpHistList->SetOwner();
740b6bde 591
38b8c336 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);
740b6bde 594
38b8c336 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);
740b6bde 597
38b8c336 598 fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
740b6bde 599 fHmpHistList->Add(fHmpCkovPhmp);
600
e7c33863 601 fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",400,0,20);
38b8c336 602 fHmpHistList->Add(fHmpMipTrkDist);
740b6bde 603
38b8c336 604 fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
605 fHmpHistList->Add(fHmpMipTrkDistX);
740b6bde 606
38b8c336 607 fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
608 fHmpHistList->Add(fHmpMipTrkDistY);
740b6bde 609
38b8c336 610 fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
611 fHmpHistList->Add(fHmpMipCharge3cm);
740b6bde 612
38b8c336 613 fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
614 fHmpHistList->Add(fHmpMipCharge1cm);
740b6bde 615
38b8c336 616 fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
617 fHmpHistList->Add(fHmpNumPhots);
740b6bde 618
38b8c336 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);
740b6bde 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
e7c33863 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
740b6bde 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
e7c33863 662// OpenFile(2);
f5dd9ecb 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]);
e7c33863 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);
9a8aafa5 695}
f5dd9ecb 696
dcb23f7e 697//____________________________________________________________________________________________________________________________________
698Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
699{
740b6bde 700 return abs(x - y) <= tolerance ;
701}
dcb23f7e 702
9a8aafa5 703#endif