1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
15 //------------------------------------------------------------------------------
\r
16 // AlidNdPtEfficiency class.
\r
18 // a. functionality:
\r
19 // - fills generic cut histograms
\r
20 // - generates cuts (selection criteria)
\r
23 // - generic cut histograms
\r
24 // - control histograms
\r
26 // Author: J.Otwinowski 18/11/2010
\r
27 //------------------------------------------------------------------------------
\r
31 #include "AliHeader.h"
\r
32 #include "AliGenEventHeader.h"
\r
33 #include "AliStack.h"
\r
34 #include "AliESDEvent.h"
\r
35 #include "AliMCEvent.h"
\r
36 #include "AliESDtrackCuts.h"
\r
37 #include "AliLog.h"
\r
38 #include "AliTracker.h"
\r
40 #include "AlidNdPtEventCuts.h"
\r
41 #include "AlidNdPtAcceptanceCuts.h"
\r
42 #include "AlidNdPtBackgroundCuts.h"
\r
43 #include "AlidNdPtAnalysis.h"
\r
44 #include "AliPhysicsSelection.h"
\r
46 #include "AliPWG0Helper.h"
\r
47 #include "AlidNdPtHelper.h"
\r
48 #include "AlidNdPtEfficiency.h"
\r
50 using namespace std;
\r
52 ClassImp(AlidNdPtEfficiency)
\r
54 //_____________________________________________________________________________
\r
55 AlidNdPtEfficiency::AlidNdPtEfficiency(): AlidNdPt(),
\r
57 fRecMCTrackHistTPCITS(0),
\r
58 fRecMCTrackHistITSTPC(0)
\r
60 // default constructor
\r
64 //_____________________________________________________________________________
\r
65 AlidNdPtEfficiency::AlidNdPtEfficiency(Char_t* name, Char_t* title): AlidNdPt(name,title),
\r
67 fRecMCTrackHistTPCITS(0),
\r
68 fRecMCTrackHistITSTPC(0)
\r
74 //_____________________________________________________________________________
\r
75 AlidNdPtEfficiency::~AlidNdPtEfficiency() {
\r
77 if(fRecMCTrackHistTPCITS) delete fRecMCTrackHistTPCITS; fRecMCTrackHistTPCITS=0;
\r
78 if(fRecMCTrackHistITSTPC) delete fRecMCTrackHistITSTPC; fRecMCTrackHistITSTPC=0;
\r
80 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
\r
83 //_____________________________________________________________________________
\r
84 void AlidNdPtEfficiency::Init(){
\r
88 const Int_t ptNbins = 63;
\r
89 const Double_t ptMin = 0.;
\r
90 const Double_t ptMax = 20.;
\r
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};
\r
95 // THnSparse track histograms
\r
98 // TPC -> ITS matching efficiency
\r
99 // eta:phi:pt:isPrim:charge:isMatch:isTPC
\r
100 Int_t binsRecMCTrackHistTPCITS[7]= { 30, 90, ptNbins, 2, 3, 2, 2 };
\r
101 Double_t minRecMCTrackHistTPCITS[7]={-1.5, 0., ptMin, 0., -1., 0., 0. };
\r
102 Double_t maxRecMCTrackHistTPCITS[7]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2., 2. };
\r
104 fRecMCTrackHistTPCITS = new THnSparseF("fRecMCTrackHistTPCITS","eta:phi:pt:isPrim:charge:isMatch:isTPC",7,binsRecMCTrackHistTPCITS,minRecMCTrackHistTPCITS,maxRecMCTrackHistTPCITS);
\r
105 fRecMCTrackHistTPCITS->SetBinEdges(2,binsPt);
\r
106 fRecMCTrackHistTPCITS->GetAxis(0)->SetTitle("#eta");
\r
107 fRecMCTrackHistTPCITS->GetAxis(1)->SetTitle("#phi (rad)");
\r
108 fRecMCTrackHistTPCITS->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
\r
109 fRecMCTrackHistTPCITS->GetAxis(3)->SetTitle("isPrim");
\r
110 fRecMCTrackHistTPCITS->GetAxis(4)->SetTitle("charge");
\r
111 fRecMCTrackHistTPCITS->GetAxis(5)->SetTitle("isMatch");
\r
112 fRecMCTrackHistTPCITS->GetAxis(6)->SetTitle("isTPC");
\r
113 fRecMCTrackHistTPCITS->Sumw2();
\r
115 // ITS -> TPC matching efficiency
\r
116 // eta:phi:pt:isPrim:charge:isMatch
\r
117 Int_t binsRecMCTrackHistITSTPC[6]= { 30, 90, ptNbins, 2, 3, 2 };
\r
118 Double_t minRecMCTrackHistITSTPC[6]={-1.5, 0., ptMin, 0., -1., 0 };
\r
119 Double_t maxRecMCTrackHistITSTPC[6]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
\r
121 fRecMCTrackHistITSTPC = new THnSparseF("fRecMCTrackHistITSTPC","eta:phi:pt:isPrim:charge:isMatch",6,binsRecMCTrackHistITSTPC,minRecMCTrackHistITSTPC,maxRecMCTrackHistITSTPC);
\r
122 fRecMCTrackHistITSTPC->SetBinEdges(2,binsPt);
\r
123 fRecMCTrackHistITSTPC->GetAxis(0)->SetTitle("#eta");
\r
124 fRecMCTrackHistITSTPC->GetAxis(1)->SetTitle("#phi (rad)");
\r
125 fRecMCTrackHistITSTPC->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
\r
126 fRecMCTrackHistITSTPC->GetAxis(3)->SetTitle("isPrim");
\r
127 fRecMCTrackHistITSTPC->GetAxis(4)->SetTitle("charge");
\r
128 fRecMCTrackHistITSTPC->GetAxis(5)->SetTitle("isMatch");
\r
129 fRecMCTrackHistITSTPC->Sumw2();
\r
131 // init output folder
\r
132 fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
\r
135 //_____________________________________________________________________________
\r
136 void AlidNdPtEfficiency::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent)
\r
139 // Process real and/or simulated events
\r
142 AliDebug(AliLog::kError, "esdEvent not available");
\r
146 // get selection cuts
\r
147 AlidNdPtEventCuts *evtCuts = GetEventCuts();
\r
148 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
\r
149 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
\r
151 if(!evtCuts || !accCuts || !esdTrackCuts) {
\r
152 AliDebug(AliLog::kError, "cuts not available");
\r
158 // trigger selection
\r
159 Bool_t isEventTriggered = kTRUE;
\r
161 // use MC information
\r
162 AliHeader* header = 0;
\r
163 AliGenEventHeader* genHeader = 0;
\r
164 AliStack* stack = 0;
\r
167 Int_t multMCTrueTracks = 0;
\r
172 AliDebug(AliLog::kError, "mcEvent not available");
\r
175 // get MC event header
\r
176 header = mcEvent->Header();
\r
178 AliDebug(AliLog::kError, "Header not available");
\r
181 // MC particle stack
\r
182 stack = mcEvent->Stack();
\r
184 AliDebug(AliLog::kError, "Stack not available");
\r
189 genHeader = header->GenEventHeader();
\r
191 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
\r
194 genHeader->PrimaryVertex(vtxMC);
\r
196 // multipliticy of all MC primary tracks
\r
197 // in Zv, pt and eta ranges)
\r
198 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
\r
202 // get reconstructed vertex
\r
203 const AliESDVertex* vtxESD = 0;
\r
204 if(evtCuts->IsRecVertexRequired())
\r
206 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
\r
207 vtxESD = esdEvent->GetPrimaryVertexTPC();
\r
209 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
\r
210 vtxESD = esdEvent->GetPrimaryVertexTracks();
\r
216 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
\r
217 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
\r
218 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
\r
220 TObjArray *allChargedTracks=0;
\r
221 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);
\r
223 // check event cuts
\r
224 if(isEventOK && isEventTriggered)
\r
226 // get all charged tracks
\r
227 allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
\r
228 if(!allChargedTracks) return;
\r
230 Int_t entries = allChargedTracks->GetEntries();
\r
231 Bool_t isTPC = kFALSE;
\r
232 Bool_t isMatch = kFALSE;
\r
234 // TPC -> ITS prolongation efficiency
\r
235 for(Int_t iTrack=0; iTrack<entries;++iTrack)
\r
237 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
\r
238 if(!track) continue;
\r
242 if(track->Charge()==0) continue;
\r
243 if(!track->GetTPCInnerParam()) continue;
\r
244 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) continue;
\r
246 // check loose cuts for TPC tracks
\r
247 if(!esdTrackCuts->AcceptTrack(track)) { continue; }
\r
251 if( (track->GetStatus()&AliESDtrack::kITSrefit) &&
\r
252 (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) )
\r
258 FillHistograms(track, stack, isMatch, isTPC, kFALSE);
\r
259 //if(tpcTrack) delete tpcTrack;
\r
263 // ITS -> TPC prolongation efficiency
\r
265 for(Int_t iTrack=0; iTrack<entries;++iTrack)
\r
267 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
\r
268 if(!track) continue;
\r
271 // ITS stand alone tracks
\r
273 if(!(track->GetStatus() & AliESDtrack::kITSpureSA)) continue;
\r
274 if(!(track->GetStatus() & AliESDtrack::kITSrefit)) continue;
\r
275 if(track->GetNcls(0)<4) continue;
\r
276 if(!track->HasPointOnITSLayer(0) && !track->HasPointOnITSLayer(1)) continue;
\r
278 // Check matching with TPC only track
\r
279 for(Int_t jTrack=0; jTrack<entries;++jTrack)
\r
283 if(iTrack==jTrack) continue;
\r
285 AliESDtrack *track2 = (AliESDtrack*)allChargedTracks->At(jTrack);
\r
286 if(!track2) continue;
\r
287 if(track2->Charge()==0) continue;
\r
288 if(!track2->GetTPCInnerParam()) continue;
\r
289 if(!(track2->GetStatus() & AliESDtrack::kTPCrefit)) continue;
\r
291 // Get TPC only tracks (must be deleted by user)
\r
292 AliESDtrack* tpcTrack2 = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, jTrack);
\r
293 if(!tpcTrack2) continue;
\r
294 if(!tpcTrack2->RelateToVertex(vtxESD,esdEvent->GetMagneticField(),100.)) { delete tpcTrack2; continue; }
\r
296 // check loose cuts for TPC tracks
\r
297 if(!esdTrackCuts->AcceptTrack(tpcTrack2)) { delete tpcTrack2; continue; }
\r
300 if (TMath::Abs(track->GetY() - tpcTrack2->GetY()) > 3) { delete tpcTrack2; continue; }
\r
301 if (TMath::Abs(track->GetSnp() - tpcTrack2->GetSnp()) > 0.2) { delete tpcTrack2; continue; }
\r
302 if (TMath::Abs(track->GetTgl() - tpcTrack2->GetTgl()) > 0.2) { delete tpcTrack2; continue; }
\r
312 FillHistograms(track, stack, isMatch, kFALSE, kTRUE);
\r
316 if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
\r
320 //_____________________________________________________________________________
\r
321 void AlidNdPtEfficiency::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, const Bool_t isMatch, const Bool_t isTPC, const Bool_t isITSTPC) const
\r
324 // Fill ESD track and MC histograms
\r
326 if(!esdTrack) return;
\r
327 Int_t charge = esdTrack->Charge();
\r
328 if(charge == 0.) return;
\r
330 Float_t pt = esdTrack->Pt();
\r
331 Float_t eta = esdTrack->Eta();
\r
332 Float_t phi = esdTrack->Phi();
\r
335 // Fill rec vs MC information
\r
337 Bool_t isPrim = kTRUE;
\r
339 if(IsUseMCInfo()) {
\r
341 Int_t label = esdTrack->GetLabel();
\r
342 if(label < 0.) return; // fake ITS track
\r
343 TParticle* particle = stack->Particle(label);
\r
344 if(!particle) return;
\r
345 if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;
\r
346 isPrim = stack->IsPhysicalPrimary(label);
\r
350 Double_t vRecMCTrackHist[6] = { eta,phi,pt,isPrim,charge,isMatch };
\r
351 Double_t vRecMCTrackHistTPCITS[7] = { eta,phi,pt,isPrim,charge,isMatch,isTPC };
\r
354 fRecMCTrackHistITSTPC->Fill(vRecMCTrackHist);
\r
357 fRecMCTrackHistTPCITS->Fill(vRecMCTrackHistTPCITS);
\r
361 //_____________________________________________________________________________
\r
362 Long64_t AlidNdPtEfficiency::Merge(TCollection* const list)
\r
364 // Merge list of objects (needed by PROOF)
\r
369 if (list->IsEmpty())
\r
372 TIterator* iter = list->MakeIterator();
\r
375 // collection of generated histograms
\r
377 while((obj = iter->Next()) != 0) {
\r
378 AlidNdPtEfficiency* entry = dynamic_cast<AlidNdPtEfficiency*>(obj);
\r
379 if (entry == 0) continue;
\r
382 fRecMCTrackHistTPCITS->Add(entry->fRecMCTrackHistTPCITS);
\r
383 fRecMCTrackHistITSTPC->Add(entry->fRecMCTrackHistITSTPC);
\r
391 //_____________________________________________________________________________
\r
392 void AlidNdPtEfficiency::Analyse()
\r
395 // Analyse histograms
\r
397 TH1::AddDirectory(kFALSE);
\r
398 TObjArray *aFolderObj = new TObjArray;
\r
399 if(!aFolderObj) return;
\r
409 AlidNdPtEventCuts *evtCuts = GetEventCuts();
\r
410 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
\r
411 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
\r
413 if(!evtCuts || !accCuts || !esdTrackCuts) {
\r
414 Error("AlidNdPtEfficiency::Analyse()", "cuts not available");
\r
419 // TPC->ITS efficiency
\r
423 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
\r
424 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
\r
425 if(!h1Dall) return;
\r
426 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
427 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
\r
430 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_TPCITS");
\r
431 h1Dc->Divide(h1Dall);
\r
432 aFolderObj->Add(h1Dc);
\r
433 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
434 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
\r
437 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
\r
438 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
439 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
\r
440 if(!h1Dall) return;
\r
441 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
442 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
\r
445 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_TPCITS");
\r
446 h1Dc->Divide(h1Dall);
\r
447 aFolderObj->Add(h1Dc);
\r
448 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
449 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
\r
452 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(2,2);
\r
453 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
454 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
\r
455 if(!h1Dall) return;
\r
456 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
457 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
\r
460 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_TPCITS");
\r
461 h1Dc->Divide(h1Dall);
\r
462 aFolderObj->Add(h1Dc);
\r
463 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
464 fRecMCTrackHistTPCITS->GetAxis(6)->SetRange(1,2);
\r
468 // ITS->TPC efficiency
\r
471 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-1.5, 1.499);
\r
472 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
475 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
\r
476 if(!h1Dall) return;
\r
477 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
478 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
\r
481 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_ITSTPC");
\r
482 h1Dc->Divide(h1Dall);
\r
483 aFolderObj->Add(h1Dc);
\r
484 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
487 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
488 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
\r
489 if(!h1Dall) return;
\r
490 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
491 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
\r
494 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_ITSTPC");
\r
495 h1Dc->Divide(h1Dall);
\r
496 aFolderObj->Add(h1Dc);
\r
497 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
500 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
501 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
\r
502 if(!h1Dall) return;
\r
503 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
504 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
\r
507 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_ITSTPC");
\r
508 h1Dc->Divide(h1Dall);
\r
509 aFolderObj->Add(h1Dc);
\r
510 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
512 // export objects to analysis folder
\r
513 fAnalysisFolder = ExportToFolder(aFolderObj);
\r
514 if(!fAnalysisFolder) {
\r
515 if(aFolderObj) delete aFolderObj;
\r
519 // delete only TObjArray
\r
520 if(aFolderObj) delete aFolderObj;
\r
523 //_____________________________________________________________________________
\r
524 TFolder* AlidNdPtEfficiency::ExportToFolder(TObjArray * const array)
\r
526 // recreate folder avery time and export objects to new one
\r
528 AlidNdPtEfficiency * comp=this;
\r
529 TFolder *folder = comp->GetAnalysisFolder();
\r
531 TString name, title;
\r
532 TFolder *newFolder = 0;
\r
534 Int_t size = array->GetSize();
\r
537 // get name and title from old folder
\r
538 name = folder->GetName();
\r
539 title = folder->GetTitle();
\r
545 newFolder = CreateFolder(name.Data(),title.Data());
\r
546 newFolder->SetOwner();
\r
548 // add objects to folder
\r
550 newFolder->Add(array->At(i));
\r
558 //_____________________________________________________________________________
\r
559 TFolder* AlidNdPtEfficiency::CreateFolder(TString name,TString title) {
\r
560 // create folder for analysed histograms
\r
562 TFolder *folder = 0;
\r
563 folder = new TFolder(name.Data(),title.Data());
\r