]> 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
Deleted wrong file
[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 ){
314 if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
315 if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
316 fHmpNumPhots->Fill(nph);
317 fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
318 if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
319
320 if (fUseMC && dist<distCut && TMath::Abs(th)<thetaCut){
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 }
331 }
332 if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),track->GetHMPIDsignal());
333 if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),track->GetHMPIDsignal());
334
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);
338 if (pdgCode==2212){
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));
343 }
344 if (pdgCode!=211) fPionCon->Fill(pBin);
345 if (pdgCode!=321) fKaonCon->Fill(pBin);
346 if (pdgCode==211){
347 if (pBin < 6){
348 Float_t weight=pPid->GetProbability(AliPID::kElectron)+
349 pPid->GetProbability(AliPID::kMuon)+
350 pPid->GetProbability(AliPID::kPion);
351 fPionTot->Fill(pBin);
352 fPionEff->Fill(pBin,weight);
353 fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
354 }
355 fProtNot->Fill(pBin,pPid->GetProbability(AliPID::kProton));
356 }
357 if (pdgCode==321){
358 if (pBin < 6){
359 fKaonTot->Fill(pBin);
360 fKaonEff->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
361 fPionNot->Fill(pBin,(pPid->GetProbability(AliPID::kPion)));
362 }
363 fProtNot->Fill(pBin,(pPid->GetProbability(AliPID::kProton)));
364 }
365 }
366 }//there is signal
367 fVar[0] = track->GetHMPIDcluIdx()/1000000;
368 fVar[1] = pHmp3;
369 fVar[2] = (Float_t)track->P();
370 fVar[3] = xpc;
371 fVar[4] = ypc;
372 fVar[5] = x;
373 fVar[6] = y;
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;
395 fTree->Fill();
396 }//track loop
397 delete pPid;
398
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
420 /* PostData(0) is taken care of by AliAnalysisTaskSE */
421 PostData(1,fHmpHistList);
422 PostData(2,fTree);
423}
424
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
432 Info("Terminate"," ");
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 }
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
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();
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
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
533 TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
534 pCan->Divide(2,3);
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
545 pCan->cd(3);
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
554 pCan->cd(5);
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
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
572 TFile *outFile = new TFile("HmpidGraphs.root","recreate");
573 pCan->Write();
574 outFile->Close();
575
576 AliAnalysisTaskSE::Terminate();
577
578}
579
580//___________________________________________________________________________
581void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
582 //
583 //HERE ONE CAN CREATE OUTPUT OBJECTS
584 //TO BE SET BEFORE THE EXECUTION OF THE TASK
585 //
586
587 //slot #1
588// OpenFile(1);
589 fHmpHistList = new TList();
590 fHmpHistList->SetOwner();
591
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);
594
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);
597
598 fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
599 fHmpHistList->Add(fHmpCkovPhmp);
600
601 fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",400,0,20);
602 fHmpHistList->Add(fHmpMipTrkDist);
603
604 fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
605 fHmpHistList->Add(fHmpMipTrkDistX);
606
607 fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
608 fHmpHistList->Add(fHmpMipTrkDistY);
609
610 fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
611 fHmpHistList->Add(fHmpMipCharge3cm);
612
613 fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
614 fHmpHistList->Add(fHmpMipCharge1cm);
615
616 fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
617 fHmpHistList->Add(fHmpNumPhots);
618
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);
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
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
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
662// OpenFile(2);
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]);
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);
695}
696
697//____________________________________________________________________________________________________________________________________
698Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
699{
700 return abs(x - y) <= tolerance ;
701}
702
703#endif