]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/AliPerformanceEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
trigger class added (M.Ivanov)
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliPerformanceEff.cxx
CommitLineData
e95d4942 1//------------------------------------------------------------------------------
2// Implementation of AliPerformanceEff class. It keeps information from
3// comparison of reconstructed and MC particle tracks. In addtion,
4// it keeps selection cuts used during comparison. The comparison
5// information is stored in the ROOT histograms. Analysis of these
6// histograms can be done by using Analyse() class function. The result of
7// the analysis (histograms/graphs) are stored in the folder which is
8// a data member of AliPerformanceEff.
9//
10// Author: J.Otwinowski 04/02/2008
11//------------------------------------------------------------------------------
12
13/*
14
15 // after running comparison task, read the file, and get component
16 gROOT->LoadMacro("$ALICE_ROOT/PWGPP/Macros/LoadMyLibs.C");
17 LoadMyLibs();
18 TFile f("Output.root");
19 AliPerformanceEff * compObj = (AliPerformanceEff*)coutput->FindObject("AliPerformanceEff");
20
21 // Analyse comparison data
22 compObj->Analyse();
23
24 // the output histograms/graphs will be stored in the folder "folderEff"
25 compObj->GetAnalysisFolder()->ls("*");
26
27 // user can save whole comparison object (or only folder with anlysed histograms)
28 // in the seperate output file (e.g.)
29 TFile fout("Analysed_Eff.root","recreate");
30 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
31 fout.Close();
32
33*/
34
35#include <TAxis.h>
36#include <TH1D.h>
37#include "THnSparse.h"
38
39//
40#include "AliESDtrack.h"
41#include "AliRecInfoCuts.h"
42#include "AliMCInfoCuts.h"
43#include "AliLog.h"
44#include "AliESDVertex.h"
45#include "AliExternalTrackParam.h"
46#include "AliTracker.h"
47#include "AliESDEvent.h"
48#include "AliMCEvent.h"
49#include "AliMCParticle.h"
50#include "AliHeader.h"
51#include "AliGenEventHeader.h"
52#include "AliStack.h"
53#include "AliPerformanceEff.h"
54
55using namespace std;
56
57ClassImp(AliPerformanceEff)
58
59//_____________________________________________________________________________
60AliPerformanceEff::AliPerformanceEff():
61 AliPerformanceObject("AliPerformanceEff"),
62
63 // histograms
64 fEffHisto(0),
65 fEffSecHisto(0),
66
67 // Cuts
68 fCutsRC(0),
69 fCutsMC(0),
70
71 // histogram folder
72 fAnalysisFolder(0)
73{
74 // default consttructor
75 Init();
76}
77
78//_____________________________________________________________________________
79AliPerformanceEff::AliPerformanceEff(Char_t* name="AliPerformanceEff",Char_t*title="AliPerformanceEff",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
80 AliPerformanceObject(name,title),
81
82 // histograms
83 fEffHisto(0),
84 fEffSecHisto(0),
85
86 // Cuts
87 fCutsRC(0),
88 fCutsMC(0),
89
90 // histogram folder
91 fAnalysisFolder(0)
92{
93 // named constructor
94 //
95 SetAnalysisMode(analysisMode);
96 SetHptGenerator(hptGenerator);
97
98 Init();
99}
100
101
102//_____________________________________________________________________________
103AliPerformanceEff::~AliPerformanceEff()
104{
105// destructor
106
107 if(fEffHisto) delete fEffHisto; fEffHisto=0;
108 if(fEffSecHisto) delete fEffSecHisto; fEffSecHisto=0;
109 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
110}
111
112//_____________________________________________________________________________
113void AliPerformanceEff::Init()
114{
115 // Init histograms
116 //
117
118 // set pt bins
119 Int_t nPtBins = 50;
120 Double_t ptMin = 1.e-2, ptMax = 20.;
121
122 Double_t *binsPt = 0;
123
124 if (IsHptGenerator()) {
125 ptMax = 100.;
126 }
127 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
128
129 /*
130 Int_t nPtBins = 31;
131 Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
132 Double_t ptMin = 0., ptMax = 10.;
133
134 if(IsHptGenerator() == kTRUE) {
135 nPtBins = 100;
136 ptMin = 0.; ptMax = 100.;
137 }
138 */
139
140 //mceta:mcphi:mcpt:pid:recStatus:findable:charge
141 Int_t binsEffHisto[9]={30,144,nPtBins,5,2,2,3,fgkMaxClones+1,fgkMaxFakes+1};
142 Double_t minEffHisto[9]={-1.5,0.,ptMin,0.,0.,0.,-1.5,0,0};
143 Double_t maxEffHisto[9]={ 1.5,2.*TMath::Pi(), ptMax,5.,2.,2.,1.5,fgkMaxClones,fgkMaxFakes};
144
145 fEffHisto = new THnSparseF("fEffHisto","mceta:mcphi:mcpt:pid:recStatus:findable:charge:nclones:nfakes",9,binsEffHisto,minEffHisto,maxEffHisto);
146 fEffHisto->SetBinEdges(2,binsPt);
147
148 fEffHisto->GetAxis(0)->SetTitle("#eta_{mc}");
149 fEffHisto->GetAxis(1)->SetTitle("#phi_{mc} (rad)");
150 fEffHisto->GetAxis(2)->SetTitle("p_{Tmc} (GeV/c)");
151 fEffHisto->GetAxis(3)->SetTitle("pid");
152 fEffHisto->GetAxis(4)->SetTitle("recStatus");
153 fEffHisto->GetAxis(5)->SetTitle("findable");
154 fEffHisto->GetAxis(6)->SetTitle("charge");
155 fEffHisto->GetAxis(7)->SetTitle("nClones");
156 fEffHisto->GetAxis(8)->SetTitle("nFakes");
157 fEffHisto->Sumw2();
158
159 //mceta:mcphi:mcpt:pid:recStatus:findable:mcR:mother_phi:mother_eta:charge
160 Int_t binsEffSecHisto[12]={30,60,nPtBins,5,2,2,100,60,30,3,fgkMaxClones+1,fgkMaxFakes+1};
161 Double_t minEffSecHisto[12]={-1.5,0.,ptMin,0.,0.,0.,0.,0.,-1.5,-1.5,0,0};
162 Double_t maxEffSecHisto[12]={ 1.5,2.*TMath::Pi(), ptMax,5.,2.,2.,200,2.*TMath::Pi(),1.5,1.5,fgkMaxClones,fgkMaxFakes};
163
164 fEffSecHisto = new THnSparseF("fEffSecHisto","mceta:mcphi:mcpt:pid:recStatus:findable:mcR:mother_phi:mother_eta:charge:nclones:nfakes",12,binsEffSecHisto,minEffSecHisto,maxEffSecHisto);
165 fEffSecHisto->SetBinEdges(2,binsPt);
166
167 fEffSecHisto->GetAxis(0)->SetTitle("#eta_{mc}");
168 fEffSecHisto->GetAxis(1)->SetTitle("#phi_{mc} (rad)");
169 fEffSecHisto->GetAxis(2)->SetTitle("p_{Tmc} (GeV/c)");
170 fEffSecHisto->GetAxis(3)->SetTitle("pid");
171 fEffSecHisto->GetAxis(4)->SetTitle("recStatus");
172 fEffSecHisto->GetAxis(5)->SetTitle("findable");
173 fEffSecHisto->GetAxis(6)->SetTitle("mcR (cm)");
174 fEffSecHisto->GetAxis(7)->SetTitle("mother_phi (rad)");
175 fEffSecHisto->GetAxis(8)->SetTitle("mother_eta");
176 fEffSecHisto->GetAxis(9)->SetTitle("charge");
177 fEffSecHisto->GetAxis(10)->SetTitle("nClones");
178 fEffSecHisto->GetAxis(11)->SetTitle("nFakes");
179 fEffSecHisto->Sumw2();
180
181 // init cuts
182 if(!fCutsMC)
183 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
184 if(!fCutsRC)
185 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
186
187 // init folder
188 fAnalysisFolder = CreateFolder("folderEff","Analysis Efficiency Folder");
189}
190
191//_____________________________________________________________________________
192void AliPerformanceEff::ProcessTPC(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
193{
194 // Fill TPC only efficiency comparison information
195 if(!esdEvent) return;
196 if(!mcEvent) return;
197
198 AliStack *stack = mcEvent->Stack();
199 if (!stack) {
200 AliDebug(AliLog::kError, "Stack not available");
201 return;
202 }
203
204 Int_t *labelsRec = NULL;
205 labelsRec = new Int_t[esdEvent->GetNumberOfTracks()];
206 if(!labelsRec)
207 {
208 Printf("Cannot create labelsRec");
209 return;
210 }
211 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRec[i] = 0; }
212
213 Int_t *labelsAllRec = NULL;
214 labelsAllRec = new Int_t[esdEvent->GetNumberOfTracks()];
215 if(!labelsAllRec) {
216 delete [] labelsRec;
217 Printf("Cannot create labelsAllRec");
218 return;
219 }
220 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRec[i] = 0; }
221
222 // loop over rec. tracks
223 AliESDtrack *track=0;
224 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
225 {
226 track = esdEvent->GetTrack(iTrack);
227 if(!track) continue;
228 if(track->Charge()==0) continue;
229
230 // if not fUseKinkDaughters don't use tracks with kink index > 0
231 if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
232
233 //Int_t label = TMath::Abs(track->GetLabel());
234 Int_t label = track->GetLabel();
235 labelsAllRec[iTrack]=label;
236
237 // TPC only
238 if(IsRecTPC(track) != 0) {
239 labelsRec[iTrack]=label;
240 }
241
242 }
243
244 //
245 // MC histograms for efficiency studies
246 //
247 Int_t nPart = stack->GetNtrack();
248 //Int_t nPart = stack->GetNprimary();
249 for (Int_t iMc = 0; iMc < nPart; ++iMc)
250 {
251 if (iMc == 0) continue; //Cannot distinguish between track or fake track
252 TParticle* particle = stack->Particle(iMc);
253 if (!particle) continue;
254 if (!particle->GetPDG()) continue;
255 if (particle->GetPDG()->Charge() == 0.0) continue;
256
257 // physical primary
258 Bool_t prim = stack->IsPhysicalPrimary(iMc);
259 if(!prim) continue;
260
261 // --- check for double filling in stack
262 // use only particles with no daughters in the list of primaries
263 Int_t nDaughters = 0;// particle->GetNDaughters();
264
265 for( Int_t iDaught=0; iDaught < particle->GetNDaughters(); iDaught++ ) {
266 if( particle->GetDaughter(iDaught) < stack->GetNprimary() )
267 nDaughters++;
268 }
269
270 if( nDaughters > 0 )
271 continue;
272 // --- check for double filling in stack
273
274 /*Bool_t findable = kFALSE;
275 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
276 {
277 // check findable
278 if(iMc == labelsAllRec[iRec])
279 {
280 findable = IsFindable(mcEvent,iMc);
281 break;
282 }
283 }*/
284 Bool_t findable = IsFindable(mcEvent,iMc);
285
286 Bool_t recStatus = kFALSE;
287 Int_t nClones = 0, nFakes = 0;
288 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
289 {
290 // check reconstructed
291 if(iMc == labelsRec[iRec])
292 {
293 if (recStatus && nClones < fgkMaxClones) nClones++;
294 recStatus = kTRUE;
295 }
296 //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
297 if (labelsRec[iRec] < 0 && -labelsRec[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
298 }
299
300 // Only 5 charged particle species (e,mu,pi,K,p)
301 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
302
303 // transform Pdg to Pid
304 Int_t pid = TransformToPID(particle);
305
306 Float_t mceta = particle->Eta();
307 Float_t mcphi = particle->Phi();
308 if(mcphi<0) mcphi += 2.*TMath::Pi();
309 Float_t mcpt = particle->Pt();
310 Float_t charge = 0.;
311 if (particle->GetPDG()->Charge() < 0) charge = -1.;
312 else if (particle->GetPDG()->Charge() > 0) charge = 1.;
313
314 // Fill histograms
315 Double_t vEffHisto[9] = {mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes};
316 fEffHisto->Fill(vEffHisto);
317 }
318 if(labelsRec) delete [] labelsRec; labelsRec = 0;
319 if(labelsAllRec) delete [] labelsAllRec; labelsAllRec = 0;
320}
321
322//_____________________________________________________________________________
323void AliPerformanceEff::ProcessTPCSec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
324{
325 // Fill TPC only efficiency comparison information for secondaries
326
327 if(!esdEvent) return;
328
329 Int_t *labelsRecSec = NULL;
330 labelsRecSec = new Int_t[esdEvent->GetNumberOfTracks()];
331 if(!labelsRecSec)
332 {
333 Printf("Cannot create labelsRecSec");
334 return;
335 }
336 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecSec[i] = 0; }
337
338 Int_t *labelsAllRecSec = NULL;
339 labelsAllRecSec = new Int_t[esdEvent->GetNumberOfTracks()];
340 if(!labelsAllRecSec) {
341 delete [] labelsRecSec;
342 Printf("Cannot create labelsAllRecSec");
343 return;
344 }
345 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecSec[i] = 0; }
346
347 // loop over rec. tracks
348 AliESDtrack *track=0;
349 Int_t multAll=0, multRec=0;
350 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
351 {
352 track = esdEvent->GetTrack(iTrack);
353 if(!track) continue;
354 if(track->Charge()==0) continue;
355
356 // if not fUseKinkDaughters don't use tracks with kink index > 0
357 if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
358
359 //Int_t label = TMath::Abs(track->GetLabel());
360 Int_t label = track->GetLabel();
361 labelsAllRecSec[multAll]=label;
362 multAll++;
363
364 // TPC only
365 if(IsRecTPC(track) != 0) {
366 labelsRecSec[multRec]=label;
367 multRec++;
368 }
369 }
370
371 //
372 // MC histograms for efficiency studies
373 //
374 if(mcEvent) {
375
376 AliStack *stack = mcEvent->Stack();
377 if (stack) {
378
379 Int_t nPart = stack->GetNtrack();
380 //Int_t nPart = stack->GetNprimary();
381 for (Int_t iMc = 0; iMc < nPart; ++iMc)
382 {
383 if (iMc == 0) continue; //Cannot distinguish between track or fake track
384 TParticle* particle = stack->Particle(iMc);
385 if (!particle) continue;
386 if (!particle->GetPDG()) continue;
387 if (particle->GetPDG()->Charge() == 0.0) continue;
388
389 // physical primary
390 Bool_t prim = stack->IsPhysicalPrimary(iMc);
391
392 // only secondaries which can be reconstructed at TPC
393 if(prim) continue;
394
395 //Float_t radius = TMath::Sqrt(particle->Vx()*particle->Vx()+particle->Vy()*particle->Vy()+particle->Vz()*particle->Vz());
396 //if(radius > fCutsMC->GetMaxR()) continue;
397
398 // only secondary electrons from gamma conversion
399 //if( TMath::Abs(particle->GetPdgCode())!=fCutsMC->GetEM() || particle->GetUniqueID() != 5) continue;
400
401 /*Bool_t findable = kFALSE;
402 for(Int_t iRec=0; iRec<multAll; ++iRec)
403 {
404 // check findable
405 if(iMc == labelsAllRecSec[iRec])
406 {
407 findable = IsFindable(mcEvent,iMc);
408 break;
409 }
410 }*/
411 Bool_t findable = IsFindable(mcEvent,iMc);
412
413 Bool_t recStatus = kFALSE;
414 Int_t nClones = 0, nFakes = 0;
415 for(Int_t iRec=0; iRec<multRec; ++iRec)
416 {
417 // check reconstructed
418 if(iMc == labelsRecSec[iRec])
419 {
420 if (recStatus && nClones < fgkMaxClones) nClones++;
421 recStatus = kTRUE;
422 }
423 //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
424 if (labelsRecSec[iRec] < 0 && -labelsRecSec[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
425 }
426
427 // Only 5 charged particle species (e,mu,pi,K,p)
428 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
429
430 // transform Pdg to Pid
431 Int_t pid = TransformToPID(particle);
432
433 Float_t mceta = particle->Eta();
434 Float_t mcphi = particle->Phi();
435 if(mcphi<0) mcphi += 2.*TMath::Pi();
436 Float_t mcpt = particle->Pt();
437 Float_t mcR = particle->R();
438
439 // get info about mother
440 Int_t motherLabel = particle->GetMother(0);
441 if(motherLabel < 0) continue;
442 TParticle *mother = stack->Particle(motherLabel);
443 if(!mother) continue;
444
445 Float_t mother_eta = mother->Eta();
446 Float_t mother_phi = mother->Phi();
447 if(mother_phi<0) mother_phi += 2.*TMath::Pi();
448
449 Float_t charge = 0.;
450 if (particle->GetPDG()->Charge() < 0) charge = -1.;
451 else if (particle->GetPDG()->Charge() > 0) charge = 1.;
452
453 // Fill histograms
454 Double_t vEffSecHisto[12] = { mceta, mcphi, mcpt, pid, recStatus, findable, mcR, mother_phi, mother_eta, charge, nClones, nFakes };
455 fEffSecHisto->Fill(vEffSecHisto);
456 }
457 }
458 }
459
460 if(labelsRecSec) delete [] labelsRecSec; labelsRecSec = 0;
461 if(labelsAllRecSec) delete [] labelsAllRecSec; labelsAllRecSec = 0;
462}
463
464
465
466
467//_____________________________________________________________________________
468void AliPerformanceEff::ProcessTPCITS(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
469{
470 // Fill efficiency comparison information
471
472 if(!esdEvent) return;
473 if(!mcEvent) return;
474
475 AliStack *stack = mcEvent->Stack();
476 if (!stack) {
477 AliDebug(AliLog::kError, "Stack not available");
478 return;
479 }
480
481 Int_t *labelsRecTPCITS = NULL;
482 labelsRecTPCITS = new Int_t[esdEvent->GetNumberOfTracks()];
483 if(!labelsRecTPCITS)
484 {
485 Printf("Cannot create labelsRecTPCITS");
486 return;
487 }
488 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecTPCITS[i] = 0; }
489
490 Int_t *labelsAllRecTPCITS = NULL;
491 labelsAllRecTPCITS = new Int_t[esdEvent->GetNumberOfTracks()];
492 if(!labelsAllRecTPCITS) {
493 delete [] labelsRecTPCITS;
494 Printf("Cannot create labelsAllRecTPCITS");
495 return;
496 }
497 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecTPCITS[i] = 0; }
498
499 // loop over rec. tracks
500 AliESDtrack *track=0;
501 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
502 {
503 track = esdEvent->GetTrack(iTrack);
504 if(!track) continue;
505 if(track->Charge()==0) continue;
506
507 // if not fUseKinkDaughters don't use tracks with kink index > 0
508 if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
509
510 //Int_t label = TMath::Abs(track->GetLabel());
511 Int_t label = track->GetLabel();
512 labelsAllRecTPCITS[iTrack]=label;
513
514 // iTPC+ITS
515 if(IsRecTPCITS(track) != 0)
516 labelsRecTPCITS[iTrack]=label;
517 }
518
519 //
520 // MC histograms for efficiency studies
521 //
522 //Int_t nPart = stack->GetNtrack();
523 Int_t nPart = stack->GetNprimary();
524 for (Int_t iMc = 0; iMc < nPart; ++iMc)
525 {
526 if (iMc == 0) continue; //Cannot distinguish between track or fake track
527 TParticle* particle = stack->Particle(iMc);
528 if (!particle) continue;
529 if (!particle->GetPDG()) continue;
530 if (particle->GetPDG()->Charge() == 0.0) continue;
531
532 // physical primary
533 Bool_t prim = stack->IsPhysicalPrimary(iMc);
534 if(!prim) continue;
535
536 /*Bool_t findable = kFALSE;
537 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
538 {
539 // check findable
540 if(iMc == labelsAllRecTPCITS[iRec])
541 {
542 findable = IsFindable(mcEvent,iMc);
543 break;
544 }
545 }*/
546 Bool_t findable = IsFindable(mcEvent,iMc);
547
548 Bool_t recStatus = kFALSE;
549 Int_t nClones = 0, nFakes = 0;
550 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
551 {
552 // check reconstructed
553 if(iMc == labelsRecTPCITS[iRec])
554 {
555 if (recStatus && nClones < fgkMaxClones) nClones++;
556 recStatus = kTRUE;
557 }
558 //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
559 if (labelsRecTPCITS[iRec] < 0 && -labelsRecTPCITS[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
560 }
561
562 // Only 5 charged particle species (e,mu,pi,K,p)
563 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
564
565 // transform Pdg to Pid
566 Int_t pid = TransformToPID(particle);
567
568 Float_t mceta = particle->Eta();
569 Float_t mcphi = particle->Phi();
570 if(mcphi<0) mcphi += 2.*TMath::Pi();
571 Float_t mcpt = particle->Pt();
572
573 Float_t charge = 0.;
574 if (particle->GetPDG()->Charge() < 0) charge = -1.;
575 else if (particle->GetPDG()->Charge() > 0) charge = 1.;
576
577 // Fill histograms
578 Double_t vEffHisto[9] = { mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes};
579 fEffHisto->Fill(vEffHisto);
580 }
581
582 if(labelsRecTPCITS) delete [] labelsRecTPCITS; labelsRecTPCITS = 0;
583 if(labelsAllRecTPCITS) delete [] labelsAllRecTPCITS; labelsAllRecTPCITS = 0;
584}
585
586//_____________________________________________________________________________
587void AliPerformanceEff::ProcessConstrained(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
588{
589 // Process comparison information
590 if(!esdEvent) return;
591 if(!mcEvent) return;
592
593 AliStack *stack = mcEvent->Stack();
594 if (!stack) {
595 AliDebug(AliLog::kError, "Stack not available");
596 return;
597 }
598
599 Int_t *labelsRecConstrained = NULL;
600 labelsRecConstrained = new Int_t[esdEvent->GetNumberOfTracks()];
601 if(!labelsRecConstrained)
602 {
603 Printf("Cannot create labelsRecConstrained");
604 return;
605 }
606 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecConstrained[i] = 0; }
607
608 Int_t *labelsAllRecConstrained = NULL;
609 labelsAllRecConstrained = new Int_t[esdEvent->GetNumberOfTracks()];
610 if(!labelsAllRecConstrained) {
611 delete [] labelsRecConstrained;
612 Printf("Cannot create labelsAllRecConstrained");
613 return;
614 }
615 for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecConstrained[i] = 0; }
616
617 // loop over rec. tracks
618 AliESDtrack *track=0;
619 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
620 {
621 track = esdEvent->GetTrack(iTrack);
622 if(!track) continue;
623 if(track->Charge()==0) continue;
624
625 // if not fUseKinkDaughters don't use tracks with kink index > 0
626 if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
627
628 //Int_t label = TMath::Abs(track->GetLabel());
629 Int_t label = track->GetLabel();
630 labelsAllRecConstrained[iTrack]=label;
631
632 // Constrained
633 if(IsRecConstrained(track) != 0)
634 labelsRecConstrained[iTrack]=label;
635
636 }
637
638 //
639 // MC histograms for efficiency studies
640 //
641
642
643 //Int_t nPart = stack->GetNtrack();
644 Int_t nPart = stack->GetNprimary();
645 for (Int_t iMc = 0; iMc < nPart; ++iMc)
646 {
647 if (iMc == 0) continue; //Cannot distinguish between track or fake track
648 TParticle* particle = stack->Particle(iMc);
649 if (!particle) continue;
650 if (!particle->GetPDG()) continue;
651 if (particle->GetPDG()->Charge() == 0.0) continue;
652
653 // physical primary
654 Bool_t prim = stack->IsPhysicalPrimary(iMc);
655 if(!prim) continue;
656
657 /*Bool_t findable = kFALSE;
658 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
659 {
660 // check findable
661 if(iMc == labelsAllRecConstrained[iRec])
662 {
663 findable = IsFindable(mcEvent,iMc);
664 break;
665 }
666 }*/
667 Bool_t findable = IsFindable(mcEvent,iMc);
668
669 Bool_t recStatus = kFALSE;
670 Int_t nClones = 0, nFakes = 0;
671 for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
672 {
673 // check reconstructed
674 if(iMc == labelsRecConstrained[iRec])
675 {
676 if (recStatus && nClones < fgkMaxClones) nClones++;
677 recStatus = kTRUE;
678 }
679 //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
680 if (labelsRecConstrained[iRec] < 0 && -labelsRecConstrained[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
681 }
682
683 // Only 5 charged particle species (e,mu,pi,K,p)
684 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
685
686 // transform Pdg to Pid
687 Int_t pid = TransformToPID(particle);
688
689 Float_t mceta = particle->Eta();
690 Float_t mcphi = particle->Phi();
691 if(mcphi<0) mcphi += 2.*TMath::Pi();
692 Float_t mcpt = particle->Pt();
693
694 Float_t charge = 0.;
695 if (particle->GetPDG()->Charge() < 0) charge = -1.;
696 else if (particle->GetPDG()->Charge() > 0) charge = 1.;
697
698 // Fill histograms
699 Double_t vEffHisto[9] = { mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes };
700 fEffHisto->Fill(vEffHisto);
701 }
702
703 if(labelsRecConstrained) delete [] labelsRecConstrained; labelsRecConstrained = 0;
704 if(labelsAllRecConstrained) delete [] labelsAllRecConstrained; labelsAllRecConstrained = 0;
705}
706
707//_____________________________________________________________________________
708void AliPerformanceEff::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
709{
710 // Process comparison information
711 //
712 if(!esdEvent)
713 {
714 Error("Exec","esdEvent not available");
715 return;
716 }
717 AliHeader* header = 0;
718 AliGenEventHeader* genHeader = 0;
719 AliStack* stack = 0;
720 TArrayF vtxMC(3);
721
722 if(bUseMC)
723 {
724 if(!mcEvent) {
725 Error("Exec","mcEvent not available");
726 return;
727 }
728 // get MC event header
729 header = mcEvent->Header();
730 if (!header) {
731 Error("Exec","Header not available");
732 return;
733 }
734 // MC particle stack
735 stack = mcEvent->Stack();
736 if (!stack) {
737 Error("Exec","Stack not available");
738 return;
739 }
740 // get MC vertex
741 genHeader = header->GenEventHeader();
742 if (!genHeader) {
743 Error("Exec","Could not retrieve genHeader from Header");
744 return;
745 }
746 genHeader->PrimaryVertex(vtxMC);
747 }
748 else {
749 Error("Exec","MC information required!");
750 return;
751 }
752
753 // use ESD friends
754 if(bUseESDfriend) {
755 if(!esdFriend) {
756 Error("Exec","esdFriend not available");
757 return;
758 }
759 }
760
761 //
762 // Process events
763 //
764 if(GetAnalysisMode() == 0) ProcessTPC(mcEvent,esdEvent);
765 else if(GetAnalysisMode() == 1) ProcessTPCITS(mcEvent,esdEvent);
766 else if(GetAnalysisMode() == 2) ProcessConstrained(mcEvent,esdEvent);
767 else if(GetAnalysisMode() == 5) ProcessTPCSec(mcEvent,esdEvent);
768 else {
769 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
770 return;
771 }
772}
773
774//_____________________________________________________________________________
775Int_t AliPerformanceEff::TransformToPID(TParticle *particle)
776{
777// transform Pdg to Pid
778// Pdg convension is different for hadrons and leptons
779// (e.g. K+/K- = 321/-321; e+/e- = -11/11 )
780
781 Int_t pid = -1;
782 if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetEM() ) pid = 0;
783 if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetMuM() ) pid = 1;
784 if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
785 if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
786 if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
787
788return pid;
789}
790
791//_____________________________________________________________________________
792Bool_t AliPerformanceEff::IsFindable(const AliMCEvent *mcEvent, Int_t label)
793{
794//
795// Findfindable tracks
796//
797if(!mcEvent) return kFALSE;
798
799 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
800 if(!mcParticle) return kFALSE;
801
802 Int_t counter;
803 Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(),0.05,counter,3.0);
804 //printf("tpcTrackLength %f \n", tpcTrackLength);
805
806return (tpcTrackLength>fCutsMC->GetMinTrackLength());
807}
808
809//_____________________________________________________________________________
810Bool_t AliPerformanceEff::IsRecTPC(AliESDtrack *esdTrack)
811{
812//
813// Check whether track is reconstructed in TPC
814//
815if(!esdTrack) return kFALSE;
816
817 const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
818 if(!track) return kFALSE;
819
820 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
821 esdTrack->GetImpactParametersTPC(dca,cov);
822
823 Bool_t recStatus = kFALSE;
824 if(esdTrack->GetTPCNcls()>fCutsRC->GetMinNClustersTPC()) recStatus = kTRUE;
825
826 /*
827 if( TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
828 TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
829 {
830 recStatus = kTRUE;
831 }
832 */
833
834return recStatus;
835}
836
837//_____________________________________________________________________________
838Bool_t AliPerformanceEff::IsRecTPCITS(AliESDtrack *esdTrack)
839{
840//
841// Check whether track is reconstructed in TPCITS
842//
843if(!esdTrack) return kFALSE;
844
845 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
846 esdTrack->GetImpactParameters(dca,cov);
847
848 Bool_t recStatus = kFALSE;
849
850 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return kFALSE; // TPC refit
851 if ((esdTrack->GetStatus()&AliESDtrack::kITSrefit)==0) return kFALSE; // ITS refit
852 if (esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return kFALSE; // min. nb. ITS clusters
853 //if ((esdTrack->GetStatus()&AliESDtrack::kITSrefit)==0) return kFALSE; // ITS refit
854 //Int_t clusterITS[200];
855 //if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return kFALSE; // min. nb. ITS clusters
856
857 recStatus = kTRUE;
858 /*
859 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
860 TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
861 {
862 recStatus = kTRUE;
863 }
864 */
865
866return recStatus;
867}
868
869//_____________________________________________________________________________
870Bool_t AliPerformanceEff::IsRecConstrained(AliESDtrack *esdTrack)
871{
872//
873// Check whether track is reconstructed in IsRecConstrained
874//
875 if(!esdTrack) return kFALSE;
876
877 const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
878 if(!track) return kFALSE;
879
880 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
881 esdTrack->GetImpactParameters(dca,cov);
882 //Int_t label = TMath::Abs(esdTrack->GetLabel());
883
884 Bool_t recStatus = kFALSE;
885
886 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return kFALSE; // TPC refit
887 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return kFALSE; // min. nb. TPC clusters
888 Int_t clusterITS[200];
889 if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return kFALSE; // min. nb. ITS clusters
890
891 recStatus = kTRUE;
892
893 /*
894 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
895 TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
896 {
897 recStatus = kTRUE;
898 }
899 */
900
901return recStatus;
902}
903
904//_____________________________________________________________________________
905Long64_t AliPerformanceEff::Merge(TCollection* const list)
906{
907 // Merge list of objects (needed by PROOF)
908
909 if (!list)
910 return 0;
911
912 if (list->IsEmpty())
913 return 1;
914
915 TIterator* iter = list->MakeIterator();
916 TObject* obj = 0;
917
918 // collection of generated histograms
919
920 Int_t count=0;
921 while((obj = iter->Next()) != 0)
922 {
923 AliPerformanceEff* entry = dynamic_cast<AliPerformanceEff*>(obj);
924 if (entry == 0) continue;
925
926 fEffHisto->Add(entry->fEffHisto);
927 fEffSecHisto->Add(entry->fEffSecHisto);
928 count++;
929 }
930
931return count;
932}
933
934//_____________________________________________________________________________
935void AliPerformanceEff::Analyse()
936{
937 // Analyse comparison information and store output histograms
938 // in the folder "folderEff"
939 //
940 TH1::AddDirectory(kFALSE);
941 TObjArray *aFolderObj = new TObjArray;
942 if(!aFolderObj) return;
943 char title[256];
944
945 //
946 // efficiency vs pt
947 //
948
949 if(GetAnalysisMode() != 5) {
950
951 fEffHisto->GetAxis(0)->SetRangeUser(-0.9,0.89); // eta range
952 fEffHisto->GetAxis(2)->SetRangeUser(0.1,10.); // pt range
953
954 // rec efficiency vs pt
955 fEffHisto->GetAxis(3)->SetRangeUser(0.,3.99); // reconstructed
956
957 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
958 aFolderObj->Add(AddHistoEff(2, "ptRecEff", "rec. efficiency", 0));
959 aFolderObj->Add(AddHistoEff(2, "ptClone", "clone rate", 1));
960 aFolderObj->Add(AddHistoEff(2, "ptFake", "fake rate", 2));
961
962 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
963 aFolderObj->Add(AddHistoEff(2, "ptRecEffNeg", "rec. efficiency neg.", 0));
964
965 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
966 aFolderObj->Add(AddHistoEff(2, "ptRecEffPos", "rec. efficiency pos.", 0));
967
968 // rec efficiency vs pid vs pt
969 fEffHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
970
971 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
972 aFolderObj->Add(AddHistoEff(2, "ptRecEffPi", "rec. efficiency (pions)", 0));
973
974 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
975 aFolderObj->Add(AddHistoEff(2, "ptRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
976
977 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
978 aFolderObj->Add(AddHistoEff(2, "ptRecEffPiPos", "rec. efficiency (pions) pos.", 0));
979
980
981 fEffHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
982
983 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
984 aFolderObj->Add(AddHistoEff(2, "ptRecEffK", "rec. efficiency (kaons)", 0));
985
986 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
987 aFolderObj->Add(AddHistoEff(2, "ptRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
988
989 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
990 aFolderObj->Add(AddHistoEff(2, "ptRecEffKPos", "rec. efficiency (kaons) pos.", 0));
991
992
993 fEffHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
994
995 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
996 aFolderObj->Add(AddHistoEff(2, "ptRecEffP", "rec. efficiency (protons)", 0));
997
998 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
999 aFolderObj->Add(AddHistoEff(2, "ptRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1000
1001 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1002 aFolderObj->Add(AddHistoEff(2, "ptRecEffPPos", "rec. efficiency (protons) pos.", 0));
1003
1004 // findable efficiency vs pt
1005 fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1006 fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1007
1008 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1009 aFolderObj->Add(AddHistoEff(2, "ptRecEffF", "rec. efficiency (findable)", 0));
1010 aFolderObj->Add(AddHistoEff(2, "ptCloneF", "clone rate (findable)", 1));
1011 aFolderObj->Add(AddHistoEff(2, "ptFakeF", "fake rate (findable)", 2));
1012
1013
1014 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1015 aFolderObj->Add(AddHistoEff(2, "ptRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1016
1017 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1018 aFolderObj->Add(AddHistoEff(2, "ptRecEffFPos", "rec. efficiency (findable) pos.", 0));
1019
1020 //
1021 // efficiency vs eta
1022 //
1023
1024 fEffHisto->GetAxis(0)->SetRangeUser(-1.5,1.5); // eta range
1025 fEffHisto->GetAxis(2)->SetRangeUser(0.2,10.); // pt range
1026 fEffHisto->GetAxis(5)->SetRangeUser(0.,1.); // all
1027
1028 // rec efficiency vs eta
1029 fEffHisto->GetAxis(3)->SetRangeUser(0.,4.); // reconstructed
1030
1031 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1032 aFolderObj->Add(AddHistoEff(0, "etaRecEff", "rec. efficiency", 0));
1033 aFolderObj->Add(AddHistoEff(0, "etaClone", "clone rate", 1));
1034 aFolderObj->Add(AddHistoEff(0, "etaFake", "fake rate", 2));
1035
1036
1037 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1038 aFolderObj->Add(AddHistoEff(0, "etaRecEffNeg", "rec. efficiency neg.", 0));
1039
1040 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1041 aFolderObj->Add(AddHistoEff(0, "etaRecEffPos", "rec. efficiency pos.", 0));
1042
1043 // rec efficiency vs pid vs eta
1044 fEffHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1045
1046 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1047 aFolderObj->Add(AddHistoEff(0, "etaRecEffPi", "rec. efficiency (pions)", 0));
1048
1049 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1050 aFolderObj->Add(AddHistoEff(0, "etaRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
1051
1052 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1053 aFolderObj->Add(AddHistoEff(0, "etaRecEffPiPos", "rec. efficiency (pions) pos.", 0));
1054
1055
1056 fEffHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1057
1058 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1059 aFolderObj->Add(AddHistoEff(0, "etaRecEffK", "rec. efficiency (kaons)", 0));
1060
1061 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1062 aFolderObj->Add(AddHistoEff(0, "etaRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
1063
1064 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1065 aFolderObj->Add(AddHistoEff(0, "etaRecEffKPos", "rec. efficiency (kaons) pos.", 0));
1066
1067
1068 fEffHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1069
1070 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1071 aFolderObj->Add(AddHistoEff(0, "etaRecEffP", "rec. efficiency (protons)", 0));
1072
1073 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1074 aFolderObj->Add(AddHistoEff(0, "etaRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1075
1076 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1077 aFolderObj->Add(AddHistoEff(0, "etaRecEffPPos", "rec. efficiency (protons) pos.", 0));
1078
1079 // findable efficiency vs eta
1080 fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1081 fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1082
1083 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1084 aFolderObj->Add(AddHistoEff(0, "etaRecEffF", "rec. efficiency (findable)", 0));
1085 aFolderObj->Add(AddHistoEff(0, "etaCloneF", "clone rate (findable)", 1));
1086 aFolderObj->Add(AddHistoEff(0, "etaFakeF", "fake rate (findable)", 2));
1087
1088
1089 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1090 aFolderObj->Add(AddHistoEff(0, "etaRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1091
1092 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1093 aFolderObj->Add(AddHistoEff(0, "etaRecEffFPos", "rec. efficiency (findable) pos.", 0));
1094
1095 //
1096 // efficiency vs phi
1097 //
1098
1099 fEffHisto->GetAxis(0)->SetRangeUser(-0.9,0.9); // eta range
1100 fEffHisto->GetAxis(2)->SetRangeUser(0.2,10.); // pt range
1101 fEffHisto->GetAxis(5)->SetRangeUser(0.,1.); // all
1102
1103 // rec efficiency vs phi
1104 fEffHisto->GetAxis(3)->SetRangeUser(0.,4.); // reconstructed
1105
1106 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1107 aFolderObj->Add(AddHistoEff(1, "phiRecEff", "rec. efficiency", 0));
1108 aFolderObj->Add(AddHistoEff(1, "phiClone", "clone rate", 1));
1109 aFolderObj->Add(AddHistoEff(1, "phiFake", "fake rate", 2));
1110
1111 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1112 aFolderObj->Add(AddHistoEff(1, "phiRecEffNeg", "rec. efficiency neg.", 0));
1113
1114 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1115 aFolderObj->Add(AddHistoEff(1, "phiRecEffPos", "rec. efficiency pos.", 0));
1116
1117 // rec efficiency vs pid vs phi
1118 fEffHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1119
1120 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1121 aFolderObj->Add(AddHistoEff(1, "phiRecEffPi", "rec. efficiency (pions)", 0));
1122
1123 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1124 aFolderObj->Add(AddHistoEff(1, "phiRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
1125
1126 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1127 aFolderObj->Add(AddHistoEff(1, "phiRecEffPiPos", "rec. efficiency (pions) pos.", 0));
1128
1129
1130 fEffHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1131
1132 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1133 aFolderObj->Add(AddHistoEff(1, "phiRecEffK", "rec. efficiency (kaons)", 0));
1134
1135 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1136 aFolderObj->Add(AddHistoEff(1, "phiRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
1137
1138 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1139 aFolderObj->Add(AddHistoEff(1, "phiRecEffKPos", "rec. efficiency (kaons) pos.", 0));
1140
1141
1142 fEffHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1143
1144 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1145 aFolderObj->Add(AddHistoEff(1, "phiRecEffP", "rec. efficiency (protons)", 0));
1146
1147 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1148 aFolderObj->Add(AddHistoEff(1, "phiRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1149
1150 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1151 aFolderObj->Add(AddHistoEff(1, "phiRecEffPPos", "rec. efficiency (protons) pos.", 0));
1152
1153 // findable efficiency vs phi
1154 fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1155 fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1156
1157 fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.); // charge all
1158 aFolderObj->Add(AddHistoEff(1, "phiRecEffF", "rec. efficiency (findable)", 0));
1159 aFolderObj->Add(AddHistoEff(1, "phiCloneF", "clone rate (findable)", 1));
1160 aFolderObj->Add(AddHistoEff(1, "phiFakeF", "fake rate (findable)", 2));
1161
1162
1163 fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.); // charge negativ
1164 aFolderObj->Add(AddHistoEff(1, "phiRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1165
1166 fEffHisto->GetAxis(6)->SetRangeUser(0.,2.); // charge positiv
1167 aFolderObj->Add(AddHistoEff(1, "phiRecEffFPos", "rec. efficiency (findable) pos.", 0));
1168 }
1169 else {
1170 //
1171 Float_t minEta=-1.5, maxEta=1.5;
1172 Float_t minR=0.0, maxR=150.0;
1173 Float_t minPt=0.15, maxPt=100.0;
1174
1175 // mother eta range
1176 fEffSecHisto->GetAxis(8)->SetRangeUser(minEta,maxEta);
1177
1178 // particle creation radius range
1179 fEffSecHisto->GetAxis(6)->SetRangeUser(minR,maxR);
1180
1181 //
1182 fEffSecHisto->GetAxis(0)->SetRangeUser(minEta,maxEta);
1183 fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1184
1185 // rec efficiency vs pt
1186
1187 aFolderObj->Add(AddHistoEff(2, "ptRecEff", "rec. efficiency", 0, 1));
1188 aFolderObj->Add(AddHistoEff(2, "ptClone", "clone rate", 1, 1));
1189 aFolderObj->Add(AddHistoEff(2, "ptFake", "fake rate", 2, 1));
1190
1191 // rec efficiency vs pid vs pt
1192
1193 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.); // electrons
1194 aFolderObj->Add(AddHistoEff(2, "ptRecEffEle", "rec. efficiency (electrons)", 0, 1));
1195
1196 fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1197 aFolderObj->Add(AddHistoEff(2, "ptRecEffPi", "rec. efficiency (pions)", 0, 1));
1198
1199 fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1200 aFolderObj->Add(AddHistoEff(2, "ptRecEffK", "rec. efficiency (kaons)", 0, 1));
1201
1202 fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1203 aFolderObj->Add(AddHistoEff(2, "ptRecEffP", "rec. efficiency (protons)", 0, 1));
1204
1205 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1206
1207 // findable efficiency vs pt
1208
1209 fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1210 aFolderObj->Add(AddHistoEff(2, "ptRecEffF", "rec. efficiency (findable)", 0, 1));
1211 aFolderObj->Add(AddHistoEff(2, "ptCloneF", "clone rate (findable)", 1, 1));
1212 aFolderObj->Add(AddHistoEff(2, "ptFakeF", "fake rate (findable)", 2, 1));
1213
1214 //
1215 // efficiency vs eta
1216 //
1217 fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1218 fEffSecHisto->GetAxis(4)->SetRangeUser(0.,1.); // all
1219 fEffSecHisto->GetAxis(5)->SetRangeUser(0.,1.); // all
1220
1221 aFolderObj->Add(AddHistoEff(0, "etaRecEff", "rec. efficiency", 0, 1));
1222 aFolderObj->Add(AddHistoEff(0, "etaClone", "clone rate", 1, 1));
1223 aFolderObj->Add(AddHistoEff(0, "etaFake", "fake rate", 2, 1));
1224
1225 // rec efficiency vs pid vs eta
1226 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.); // electrons
1227 aFolderObj->Add(AddHistoEff(0, "etaRecEffEle", "rec. efficiency (electrons)", 0, 1));
1228
1229 fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1230 aFolderObj->Add(AddHistoEff(0, "etaRecEffPi", "rec. efficiency (pions)", 0, 1));
1231
1232 fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1233 aFolderObj->Add(AddHistoEff(0, "etaRecEffK", "rec. efficiency (kaons)", 0, 1));
1234
1235 fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1236 aFolderObj->Add(AddHistoEff(0, "etaRecEffP", "rec. efficiency (protons)", 0, 1));
1237
1238 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1239
1240 // findable efficiency vs eta
1241
1242 fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1243 aFolderObj->Add(AddHistoEff(0, "etaRecEffF", "rec. efficiency (findable)", 0, 1));
1244 aFolderObj->Add(AddHistoEff(0, "etaCloneF", "clone rate (findable)", 1, 1));
1245 aFolderObj->Add(AddHistoEff(0, "etaFakeF", "fake rate (findable)", 2, 1));
1246
1247 //
1248 // efficiency vs phi
1249 //
1250
1251 fEffSecHisto->GetAxis(0)->SetRangeUser(minEta,maxEta);
1252 fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1253
1254 fEffSecHisto->GetAxis(4)->SetRangeUser(0.,1.); // all
1255 fEffSecHisto->GetAxis(5)->SetRangeUser(0.,1.); // all
1256
1257 aFolderObj->Add(AddHistoEff(1, "phiRecEff", "rec. efficiency", 0, 1));
1258 aFolderObj->Add(AddHistoEff(1, "phiClone", "clone rate", 1, 1));
1259 aFolderObj->Add(AddHistoEff(1, "phiFake", "fake rate", 2, 1));
1260
1261 // rec efficiency vs pid vs phi
1262 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.);
1263 aFolderObj->Add(AddHistoEff(1, "phiRecEffEle", "rec. efficiency (electrons)", 0, 1));
1264
1265 fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1266 aFolderObj->Add(AddHistoEff(1, "phiRecEffPi", "rec. efficiency (pions)", 0, 1));
1267
1268 fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1269 aFolderObj->Add(AddHistoEff(1, "phiRecEffK", "rec. efficiency (kaons)", 0, 1));
1270
1271 fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1272 aFolderObj->Add(AddHistoEff(1, "phiRecEffP", "rec. efficiency (protons)", 0, 1));
1273
1274 fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1275
1276 // findable efficiency vs phi
1277
1278 fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1279 aFolderObj->Add(AddHistoEff(1, "phiRecEffF", "rec. efficiency (findable)", 0, 1));
1280 aFolderObj->Add(AddHistoEff(1, "phiCloneF", "clone rate (findable)", 1, 1));
1281 aFolderObj->Add(AddHistoEff(1, "phiFakeF", "fake rate (findable)", 2, 1));
1282 }
1283
1284 for (Int_t i = 0;i < fEffHisto->GetNdimensions();i++)
1285 {
1286 fEffHisto->GetAxis(i)->SetRange(1,0); //Reset Range
1287 }
1288
1289 // export objects to analysis folder
1290 fAnalysisFolder = ExportToFolder(aFolderObj);
1291
1292 // delete only TObjArray
1293 if(aFolderObj) delete aFolderObj;
1294}
1295
1296//_____________________________________________________________________________
1297TFolder* AliPerformanceEff::ExportToFolder(TObjArray * array)
1298{
1299 // recreate folder avery time and export objects to new one
1300 //
1301 AliPerformanceEff * comp=this;
1302 TFolder *folder = comp->GetAnalysisFolder();
1303
1304 TString name, title;
1305 TFolder *newFolder = 0;
1306 Int_t i = 0;
1307 Int_t size = array->GetSize();
1308
1309 if(folder) {
1310 // get name and title from old folder
1311 name = folder->GetName();
1312 title = folder->GetTitle();
1313
1314 // delete old one
1315 delete folder;
1316
1317 // create new one
1318 newFolder = CreateFolder(name.Data(),title.Data());
1319 newFolder->SetOwner();
1320
1321 // add objects to folder
1322 while(i < size) {
1323 newFolder->Add(array->At(i));
1324 i++;
1325 }
1326 }
1327
1328return newFolder;
1329}
1330
1331
1332//_____________________________________________________________________________
1333TFolder* AliPerformanceEff::CreateFolder(TString name,TString title) {
1334// create folder for analysed histograms
1335//
1336TFolder *folder = 0;
1337 folder = new TFolder(name.Data(),title.Data());
1338
1339 return folder;
1340}
1341
1342TH1D* WeightedProjection(THnSparseF* src, Int_t axis, Int_t nWeights, Int_t* weightCoords)
1343{
1344 THnSparseF* tmp = (THnSparseF*) src->Clone();
1345 Int_t i;
1346 for (i = 0;i < tmp->GetNbins();i++)
1347 {
1348 Int_t coords[12];
1349 tmp->GetBinContent(i, coords);
1350 Int_t weight = 0, j;
1351 for (j = 0;j < nWeights;j++)
1352 {
1353 //The coordinates range from 1 to maxClones / maxFakes + 1, so we have to subtract one
1354 weight += coords[weightCoords[j]] - 1;
1355 }
1356 tmp->SetBinContent(i, weight);
1357 }
1358
1359 TH1D* ret = tmp->Projection(axis);
1360 delete tmp;
1361 return(ret);
1362}
1363
1364//_____________________________________________________________________________
1365TH1D* AliPerformanceEff::AddHistoEff(Int_t axis, const Char_t *name, const Char_t* vsTitle,
1366 const Int_t type, const Int_t secondary) {
1367 // Create and add rec efficiency vs pt, eta, phi
1368
1369 char title[256];
1370
1371 TH1D *recc = NULL;
1372
1373 THnSparseF* EffHisto = secondary ? fEffSecHisto : fEffHisto;
1374
1375 Int_t axis_clone = secondary ? 10 : 7;
1376 Int_t axis_fake = secondary ? 11 : 8;
1377 Int_t axis_all[3] = {4, axis_clone, axis_fake};
1378
1379
1380 if (type == 0) // Efficiency
1381 {
1382 EffHisto->GetAxis(4)->SetRangeUser(0.,1.); // all
1383 TH1D *all = EffHisto->Projection(axis);
1384
1385 EffHisto->GetAxis(4)->SetRangeUser(1.,1.); // reconstructed
1386 TH1D *rec = EffHisto->Projection(axis);
1387 recc = (TH1D*)rec->Clone();
1388
1389 recc->Divide(rec,all,1,1,"B");
1390 recc->GetYaxis()->SetTitle("efficiency");
1391 }
1392 else if (type == 1) // Clone Rate
1393 {
1394 EffHisto->GetAxis(4)->SetRangeUser(1.,1.); // reconstructed
1395 TH1D *all = WeightedProjection(EffHisto, axis, 3, axis_all);
1396
1397 EffHisto->GetAxis(4)->SetRangeUser(1.,1.); // reconstructed
1398 TH1D *clone = WeightedProjection(EffHisto, axis, 1, &axis_clone);
1399 recc = (TH1D*) clone->Clone();
1400
1401 recc->Divide(clone,all,1,1,"B");
1402 recc->GetYaxis()->SetTitle("clone rate");
1403 }
1404 else if (type == 2) // Fake Rate
1405 {
1406 EffHisto->GetAxis(4)->SetRangeUser(1.,1.); // reconstructed
1407 TH1D *all = WeightedProjection(EffHisto, axis, 3, axis_all);
1408
1409 EffHisto->GetAxis(4)->SetRangeUser(1.,1.); // reconstructed
1410 TH1D *fake = WeightedProjection(EffHisto, axis, 1, &axis_fake);
1411 recc = (TH1D*) fake->Clone();
1412
1413 recc->Divide(fake,all,1,1,"B");
1414 recc->GetYaxis()->SetTitle("fake rate");
1415 }
1416
1417 EffHisto->GetAxis(4)->SetRange(1,0); //Reset Range
1418
1419 recc->SetName(name);
1420
1421 recc->GetXaxis()->SetTitle(fEffHisto->GetAxis(axis)->GetTitle());
1422
1423 snprintf(title,256,"%s vs %s",vsTitle, fEffHisto->GetAxis(axis)->GetTitle());
1424 recc->SetTitle(title);
1425
1426 if (axis == 2 ) recc->SetBit(TH1::kLogX);
1427
1428 return recc;
1429}