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