Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtEfficiency.cxx
CommitLineData
a65a7e70 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//------------------------------------------------------------------------------
16// AlidNdPtEfficiency class.
17//
18// a. functionality:
19// - fills generic cut histograms
20// - generates cuts (selection criteria)
21//
22// b. data members:
23// - generic cut histograms
24// - control histograms
25//
26// Author: J.Otwinowski 18/11/2010
27//------------------------------------------------------------------------------
28#include "TH1.h"
29#include "TH2.h"
30
31#include "AliHeader.h"
32#include "AliGenEventHeader.h"
33#include "AliStack.h"
34#include "AliESDEvent.h"
35#include "AliMCEvent.h"
36#include "AliESDtrackCuts.h"
37#include "AliLog.h"
38#include "AliTracker.h"
39
40#include "AlidNdPtEventCuts.h"
41#include "AlidNdPtAcceptanceCuts.h"
42#include "AlidNdPtBackgroundCuts.h"
43#include "AlidNdPtAnalysis.h"
44#include "AliPhysicsSelection.h"
45
46#include "AliPWG0Helper.h"
47#include "AlidNdPtHelper.h"
48#include "AlidNdPtEfficiency.h"
49
50using namespace std;
51
52ClassImp(AlidNdPtEfficiency)
53
54//_____________________________________________________________________________
55 AlidNdPtEfficiency::AlidNdPtEfficiency(): AlidNdPt(),
56 fAnalysisFolder(0),
57 fRecMCTrackHistTPCITS(0),
58 fRecMCTrackHistITSTPC(0)
59{
60 // default constructor
61 Init();
62}
63
64//_____________________________________________________________________________
65AlidNdPtEfficiency::AlidNdPtEfficiency(Char_t* name, Char_t* title): AlidNdPt(name,title),
66 fAnalysisFolder(0),
67 fRecMCTrackHistTPCITS(0),
68 fRecMCTrackHistITSTPC(0)
69{
70 // constructor
71 Init();
72}
73
74//_____________________________________________________________________________
75AlidNdPtEfficiency::~AlidNdPtEfficiency() {
76 //
77 if(fRecMCTrackHistTPCITS) delete fRecMCTrackHistTPCITS; fRecMCTrackHistTPCITS=0;
78 if(fRecMCTrackHistITSTPC) delete fRecMCTrackHistITSTPC; fRecMCTrackHistITSTPC=0;
79
80 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
81}
82
83//_____________________________________________________________________________
84void AlidNdPtEfficiency::Init(){
85 //
86 // Init histograms
87 //
88 const Int_t ptNbins = 63;
89 const Double_t ptMin = 0.;
90 const Double_t ptMax = 20.;
91
92 Double_t binsPt[ptNbins+1] = {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.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0, 20.,25.,30.,35.,40.,50};
93
94 //
95 // THnSparse track histograms
96 //
97
98 // TPC -> ITS matching efficiency
99 // eta:phi:pt:isPrim:charge:isMatch:isTPC
100 Int_t binsRecMCTrackHistTPCITS[7]= { 30, 90, ptNbins, 2, 3, 2, 2 };
101 Double_t minRecMCTrackHistTPCITS[7]={-1.5, 0., ptMin, 0., -1., 0., 0. };
102 Double_t maxRecMCTrackHistTPCITS[7]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2., 2. };
103
104 fRecMCTrackHistTPCITS = new THnSparseF("fRecMCTrackHistTPCITS","eta:phi:pt:isPrim:charge:isMatch:isTPC",7,binsRecMCTrackHistTPCITS,minRecMCTrackHistTPCITS,maxRecMCTrackHistTPCITS);
105 fRecMCTrackHistTPCITS->SetBinEdges(2,binsPt);
106 fRecMCTrackHistTPCITS->GetAxis(0)->SetTitle("#eta");
107 fRecMCTrackHistTPCITS->GetAxis(1)->SetTitle("#phi (rad)");
108 fRecMCTrackHistTPCITS->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
109 fRecMCTrackHistTPCITS->GetAxis(3)->SetTitle("isPrim");
110 fRecMCTrackHistTPCITS->GetAxis(4)->SetTitle("charge");
111 fRecMCTrackHistTPCITS->GetAxis(5)->SetTitle("isMatch");
112 fRecMCTrackHistTPCITS->GetAxis(6)->SetTitle("isTPC");
113 fRecMCTrackHistTPCITS->Sumw2();
114
115 // ITS -> TPC matching efficiency
116 // eta:phi:pt:isPrim:charge:isMatch
117 Int_t binsRecMCTrackHistITSTPC[6]= { 30, 90, ptNbins, 2, 3, 2 };
118 Double_t minRecMCTrackHistITSTPC[6]={-1.5, 0., ptMin, 0., -1., 0 };
119 Double_t maxRecMCTrackHistITSTPC[6]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
120
121 fRecMCTrackHistITSTPC = new THnSparseF("fRecMCTrackHistITSTPC","eta:phi:pt:isPrim:charge:isMatch",6,binsRecMCTrackHistITSTPC,minRecMCTrackHistITSTPC,maxRecMCTrackHistITSTPC);
122 fRecMCTrackHistITSTPC->SetBinEdges(2,binsPt);
123 fRecMCTrackHistITSTPC->GetAxis(0)->SetTitle("#eta");
124 fRecMCTrackHistITSTPC->GetAxis(1)->SetTitle("#phi (rad)");
125 fRecMCTrackHistITSTPC->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
126 fRecMCTrackHistITSTPC->GetAxis(3)->SetTitle("isPrim");
127 fRecMCTrackHistITSTPC->GetAxis(4)->SetTitle("charge");
128 fRecMCTrackHistITSTPC->GetAxis(5)->SetTitle("isMatch");
129 fRecMCTrackHistITSTPC->Sumw2();
130
131 // init output folder
132 fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
133}
134
135//_____________________________________________________________________________
136void AlidNdPtEfficiency::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent)
137{
138 //
139 // Process real and/or simulated events
140 //
141 if(!esdEvent) {
142 AliDebug(AliLog::kError, "esdEvent not available");
143 return;
144 }
145
146 // get selection cuts
147 AlidNdPtEventCuts *evtCuts = GetEventCuts();
148 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
149 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
150
151 if(!evtCuts || !accCuts || !esdTrackCuts) {
152 AliDebug(AliLog::kError, "cuts not available");
153 return;
154 }
155
156
157
158 // trigger selection
159 Bool_t isEventTriggered = kTRUE;
160
161 // use MC information
162 AliHeader* header = 0;
163 AliGenEventHeader* genHeader = 0;
164 AliStack* stack = 0;
165 TArrayF vtxMC(3);
166
167 Int_t multMCTrueTracks = 0;
168 if(IsUseMCInfo())
169 {
170 //
171 if(!mcEvent) {
172 AliDebug(AliLog::kError, "mcEvent not available");
173 return;
174 }
175 // get MC event header
176 header = mcEvent->Header();
177 if (!header) {
178 AliDebug(AliLog::kError, "Header not available");
179 return;
180 }
181 // MC particle stack
182 stack = mcEvent->Stack();
183 if (!stack) {
184 AliDebug(AliLog::kError, "Stack not available");
185 return;
186 }
187
188 // get MC vertex
189 genHeader = header->GenEventHeader();
190 if (!genHeader) {
191 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
192 return;
193 }
194 genHeader->PrimaryVertex(vtxMC);
195
196 // multipliticy of all MC primary tracks
197 // in Zv, pt and eta ranges)
198 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
199
200 } // end bUseMC
201
202 // get reconstructed vertex
203 const AliESDVertex* vtxESD = 0;
204 if(evtCuts->IsRecVertexRequired())
205 {
206 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
207 vtxESD = esdEvent->GetPrimaryVertexTPC();
208 }
209 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
210 vtxESD = esdEvent->GetPrimaryVertexTracks();
211 }
212 else {
213 return;
214 }
215 }
216 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
217 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
218 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
219
220 TObjArray *allChargedTracks=0;
221 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);
222
223 // check event cuts
224 if(isEventOK && isEventTriggered)
225 {
226 // get all charged tracks
227 allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
228 if(!allChargedTracks) return;
229
230 Int_t entries = allChargedTracks->GetEntries();
231 Bool_t isTPC = kFALSE;
232 Bool_t isMatch = kFALSE;
233
234 // TPC -> ITS prolongation efficiency
235 for(Int_t iTrack=0; iTrack<entries;++iTrack)
236 {
237 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
238 if(!track) continue;
239
240 isTPC = kFALSE;
241
242 if(track->Charge()==0) continue;
243 if(!track->GetTPCInnerParam()) continue;
244 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) continue;
245
246 // check loose cuts for TPC tracks
247 if(!esdTrackCuts->AcceptTrack(track)) { continue; }
248
249 isTPC = kTRUE;
250 isMatch = kFALSE;
251 if( (track->GetStatus()&AliESDtrack::kITSrefit) &&
252 (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) )
253 {
254 isMatch = kTRUE;
255 }
256
257 //
258 FillHistograms(track, stack, isMatch, isTPC, kFALSE);
259 //if(tpcTrack) delete tpcTrack;
260 }
261
262 //
263 // ITS -> TPC prolongation efficiency
264 //
265 for(Int_t iTrack=0; iTrack<entries;++iTrack)
266 {
267 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
268 if(!track) continue;
269
270 //
271 // ITS stand alone tracks
272 //
273 if(!(track->GetStatus() & AliESDtrack::kITSpureSA)) continue;
274 if(!(track->GetStatus() & AliESDtrack::kITSrefit)) continue;
275 if(track->GetNcls(0)<4) continue;
276 if(!track->HasPointOnITSLayer(0) && !track->HasPointOnITSLayer(1)) continue;
277
278 // Check matching with TPC only track
279 for(Int_t jTrack=0; jTrack<entries;++jTrack)
280 {
281 isMatch = kFALSE;
282
283 if(iTrack==jTrack) continue;
284
285 AliESDtrack *track2 = (AliESDtrack*)allChargedTracks->At(jTrack);
286 if(!track2) continue;
287 if(track2->Charge()==0) continue;
288 if(!track2->GetTPCInnerParam()) continue;
289 if(!(track2->GetStatus() & AliESDtrack::kTPCrefit)) continue;
290
291 // Get TPC only tracks (must be deleted by user)
292 AliESDtrack* tpcTrack2 = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, jTrack);
293 if(!tpcTrack2) continue;
294 if(!tpcTrack2->RelateToVertex(vtxESD,esdEvent->GetMagneticField(),100.)) { delete tpcTrack2; continue; }
295
296 // check loose cuts for TPC tracks
297 if(!esdTrackCuts->AcceptTrack(tpcTrack2)) { delete tpcTrack2; continue; }
298
299 // check matching
300 if (TMath::Abs(track->GetY() - tpcTrack2->GetY()) > 3) { delete tpcTrack2; continue; }
301 if (TMath::Abs(track->GetSnp() - tpcTrack2->GetSnp()) > 0.2) { delete tpcTrack2; continue; }
302 if (TMath::Abs(track->GetTgl() - tpcTrack2->GetTgl()) > 0.2) { delete tpcTrack2; continue; }
303
304 isMatch = kTRUE;
305 if(tpcTrack2) {
306 delete tpcTrack2;
307 }
308 break;
309 }
310
311 //
312 FillHistograms(track, stack, isMatch, kFALSE, kTRUE);
313 }
314 }
315
316 if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
317
318}
319
320//_____________________________________________________________________________
321void AlidNdPtEfficiency::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, const Bool_t isMatch, const Bool_t isTPC, const Bool_t isITSTPC) const
322{
323 //
324 // Fill ESD track and MC histograms
325 //
326 if(!esdTrack) return;
327 Int_t charge = esdTrack->Charge();
328 if(charge == 0.) return;
329
330 Float_t pt = esdTrack->Pt();
331 Float_t eta = esdTrack->Eta();
332 Float_t phi = esdTrack->Phi();
333
334 //
335 // Fill rec vs MC information
336 //
337 Bool_t isPrim = kTRUE;
338
339 if(IsUseMCInfo()) {
340 if(!stack) return;
341 Int_t label = esdTrack->GetLabel();
342 if(label < 0.) return; // fake ITS track
343 TParticle* particle = stack->Particle(label);
344 if(!particle) return;
345 if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;
346 isPrim = stack->IsPhysicalPrimary(label);
347 }
348
349 // fill histo
2942f542 350 Double_t vRecMCTrackHist[6] = { eta,phi,pt,static_cast<Double_t>(isPrim),static_cast<Double_t>(charge),static_cast<Double_t>(isMatch) };
351 Double_t vRecMCTrackHistTPCITS[7] = { eta,phi,pt,static_cast<Double_t>(isPrim),static_cast<Double_t>(charge),static_cast<Double_t>(isMatch),static_cast<Double_t>(isTPC) };
a65a7e70 352
353 if(isITSTPC) {
354 fRecMCTrackHistITSTPC->Fill(vRecMCTrackHist);
355 }
356 else {
357 fRecMCTrackHistTPCITS->Fill(vRecMCTrackHistTPCITS);
358 }
359}
360
361//_____________________________________________________________________________
362Long64_t AlidNdPtEfficiency::Merge(TCollection* const list)
363{
364 // Merge list of objects (needed by PROOF)
365
366 if (!list)
367 return 0;
368
369 if (list->IsEmpty())
370 return 1;
371
372 TIterator* iter = list->MakeIterator();
373 TObject* obj = 0;
374
375 // collection of generated histograms
376 Int_t count=0;
377 while((obj = iter->Next()) != 0) {
378 AlidNdPtEfficiency* entry = dynamic_cast<AlidNdPtEfficiency*>(obj);
379 if (entry == 0) continue;
380
381 // track histo
382 fRecMCTrackHistTPCITS->Add(entry->fRecMCTrackHistTPCITS);
383 fRecMCTrackHistITSTPC->Add(entry->fRecMCTrackHistITSTPC);
384
385 count++;
386 }
387
388return count;
389}
390
391//_____________________________________________________________________________
392void AlidNdPtEfficiency::Analyse()
393{
394 //
395 // Analyse histograms
396 //
397 TH1::AddDirectory(kFALSE);
398 TObjArray *aFolderObj = new TObjArray;
399 if(!aFolderObj) return;
400
401 TH1D *h1Dall = 0;
402 TH1D *h1D = 0;
403 TH1D *h1Dc = 0;
404
405
406 //
407 // get cuts
408 //
409 AlidNdPtEventCuts *evtCuts = GetEventCuts();
410 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
411 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
412
413 if(!evtCuts || !accCuts || !esdTrackCuts) {
414 Error("AlidNdPtEfficiency::Analyse()", "cuts not available");
415 return;
416 }
417
418 //
419 // TPC->ITS efficiency
420 //
421
422 //eff vs eta
423 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
424 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
425 if(!h1Dall) return;
426 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
427 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
428 if(!h1D) return;
429
430 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_TPCITS");
431 h1Dc->Divide(h1Dall);
432 aFolderObj->Add(h1Dc);
433 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
434 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
435
436 //eff vs phi
437 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
438 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
439 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
440 if(!h1Dall) return;
441 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
442 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
443 if(!h1D) return;
444
445 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_TPCITS");
446 h1Dc->Divide(h1Dall);
447 aFolderObj->Add(h1Dc);
448 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
449 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
450
451 //eff vs pT
452 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
453 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
454 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
455 if(!h1Dall) return;
456 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
457 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
458 if(!h1D) return;
459
460 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_TPCITS");
461 h1Dc->Divide(h1Dall);
462 aFolderObj->Add(h1Dc);
463 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
464 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
465
466
467 //
468 // ITS->TPC efficiency
469 //
470
471 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-1.5, 1.499);
472 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
473
474 //eff vs eta
475 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
476 if(!h1Dall) return;
477 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
478 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
479 if(!h1D) return;
480
481 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_ITSTPC");
482 h1Dc->Divide(h1Dall);
483 aFolderObj->Add(h1Dc);
484 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
485
486 //eff vs phi
487 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
488 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
489 if(!h1Dall) return;
490 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
491 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
492 if(!h1D) return;
493
494 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_ITSTPC");
495 h1Dc->Divide(h1Dall);
496 aFolderObj->Add(h1Dc);
497 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
498
499 //eff vs pT
500 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
501 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
502 if(!h1Dall) return;
503 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
504 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
505 if(!h1D) return;
506
507 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_ITSTPC");
508 h1Dc->Divide(h1Dall);
509 aFolderObj->Add(h1Dc);
510 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
511
512 // export objects to analysis folder
513 fAnalysisFolder = ExportToFolder(aFolderObj);
514 if(!fAnalysisFolder) {
515 if(aFolderObj) delete aFolderObj;
516 return;
517 }
518
519 // delete only TObjArray
520 if(aFolderObj) delete aFolderObj;
521}
522
523//_____________________________________________________________________________
524TFolder* AlidNdPtEfficiency::ExportToFolder(TObjArray * const array)
525{
526 // recreate folder avery time and export objects to new one
527 //
528 AlidNdPtEfficiency * comp=this;
529 TFolder *folder = comp->GetAnalysisFolder();
530
531 TString name, title;
532 TFolder *newFolder = 0;
533 Int_t i = 0;
534 Int_t size = array->GetSize();
535
536 if(folder) {
537 // get name and title from old folder
538 name = folder->GetName();
539 title = folder->GetTitle();
540
541 // delete old one
542 delete folder;
543
544 // create new one
545 newFolder = CreateFolder(name.Data(),title.Data());
546 newFolder->SetOwner();
547
548 // add objects to folder
549 while(i < size) {
550 newFolder->Add(array->At(i));
551 i++;
552 }
553 }
554
555return newFolder;
556}
557
558//_____________________________________________________________________________
559TFolder* AlidNdPtEfficiency::CreateFolder(TString name,TString title) {
560// create folder for analysed histograms
561//
562TFolder *folder = 0;
563 folder = new TFolder(name.Data(),title.Data());
564
565 return folder;
566}