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