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 = 58;
\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.};
\r
95 // THnSparse track histograms
\r
98 // TPC -> ITS matching efficiency
\r
99 // eta:phi:pt:isPrim:charge:isMatch
\r
100 Int_t binsRecMCTrackHistTPCITS[6]= { 30, 90, ptNbins, 2, 3, 2 };
\r
101 Double_t minRecMCTrackHistTPCITS[6]={-1.5, 0., ptMin, 0., -1., 0 };
\r
102 Double_t maxRecMCTrackHistTPCITS[6]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
\r
104 fRecMCTrackHistTPCITS = new THnSparseF("fRecMCTrackHistTPCITS","eta:phi:pt:isPrim:charge:isMatch",6,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->Sumw2();
\r
114 // ITS -> TPC matching efficiency
\r
115 // eta:phi:pt:isPrim:charge:isMatch
\r
116 Int_t binsRecMCTrackHistITSTPC[6]= { 30, 90, ptNbins, 2, 3, 2 };
\r
117 Double_t minRecMCTrackHistITSTPC[6]={-1.5, 0., ptMin, 0., -1., 0 };
\r
118 Double_t maxRecMCTrackHistITSTPC[6]={ 1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
\r
120 fRecMCTrackHistITSTPC = new THnSparseF("fRecMCTrackHistITSTPC","eta:phi:pt:isPrim:charge:isMatch",6,binsRecMCTrackHistITSTPC,minRecMCTrackHistITSTPC,maxRecMCTrackHistITSTPC);
\r
121 fRecMCTrackHistITSTPC->SetBinEdges(2,binsPt);
\r
122 fRecMCTrackHistITSTPC->GetAxis(0)->SetTitle("#eta");
\r
123 fRecMCTrackHistITSTPC->GetAxis(1)->SetTitle("#phi (rad)");
\r
124 fRecMCTrackHistITSTPC->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
\r
125 fRecMCTrackHistITSTPC->GetAxis(3)->SetTitle("isPrim");
\r
126 fRecMCTrackHistITSTPC->GetAxis(4)->SetTitle("charge");
\r
127 fRecMCTrackHistITSTPC->GetAxis(5)->SetTitle("isMatch");
\r
128 fRecMCTrackHistITSTPC->Sumw2();
\r
130 // init output folder
\r
131 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
156 // trigger selection
\r
157 Bool_t isEventTriggered = kTRUE;
\r
158 AliPhysicsSelection *trigSel = NULL;
\r
159 AliTriggerAnalysis *trigAna = NULL; // needed for andV0
\r
161 if(evtCuts->IsTriggerRequired())
\r
164 trigSel = GetPhysicsTriggerSelection();
\r
166 printf("cannot get trigSel \n");
\r
173 trigSel->SetAnalyzeMC();
\r
176 isEventTriggered = trigSel->IsCollisionCandidate(esdEvent);
\r
178 if(GetTrigger() == AliTriggerAnalysis::kV0AND)
\r
180 trigAna = trigSel->GetTriggerAnalysis();
\r
184 isEventTriggered = trigAna->IsOfflineTriggerFired(esdEvent, GetTrigger());
\r
185 }//if(GetTrigger() == AliTriggerAnalysis::kV0AND)
\r
186 }//if(IsUseMCInfo())
\r
187 }//if(evtCuts->IsTriggerRequired())
\r
189 // use MC information
\r
190 AliHeader* header = 0;
\r
191 AliGenEventHeader* genHeader = 0;
\r
192 AliStack* stack = 0;
\r
195 Int_t multMCTrueTracks = 0;
\r
200 AliDebug(AliLog::kError, "mcEvent not available");
\r
203 // get MC event header
\r
204 header = mcEvent->Header();
\r
206 AliDebug(AliLog::kError, "Header not available");
\r
209 // MC particle stack
\r
210 stack = mcEvent->Stack();
\r
212 AliDebug(AliLog::kError, "Stack not available");
\r
217 genHeader = header->GenEventHeader();
\r
219 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
\r
222 genHeader->PrimaryVertex(vtxMC);
\r
224 // multipliticy of all MC primary tracks
\r
225 // in Zv, pt and eta ranges)
\r
226 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
\r
230 // get reconstructed vertex
\r
231 const AliESDVertex* vtxESD = 0;
\r
232 if(evtCuts->IsRecVertexRequired())
\r
234 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
\r
235 vtxESD = esdEvent->GetPrimaryVertexTPC();
\r
237 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
\r
238 vtxESD = esdEvent->GetPrimaryVertexTracks();
\r
244 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
\r
245 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
\r
246 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
\r
248 TObjArray *allChargedTracks=0;
\r
249 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);
\r
251 // check event cuts
\r
252 if(isEventOK && isEventTriggered)
\r
254 // get all charged tracks
\r
255 allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
\r
256 if(!allChargedTracks) return;
\r
258 Bool_t isMatch = kFALSE;
\r
259 Int_t entries = allChargedTracks->GetEntries();
\r
261 // TPC -> ITS prolongation efficiency
\r
262 for(Int_t iTrack=0; iTrack<entries;++iTrack)
\r
264 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
\r
265 if(!track) continue;
\r
267 if(track->Charge()==0) continue;
\r
268 if (!track->GetTPCInnerParam()) continue;
\r
270 // Get TPC only tracks (must be deleted by user)
\r
271 AliESDtrack* tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
\r
272 if(!tpcTrack) continue;
\r
274 // check loose cuts for TPC tracks
\r
275 if(!esdTrackCuts->AcceptTrack(tpcTrack)) { delete tpcTrack; continue; }
\r
277 if(track->GetNcls(0) >= 2 && (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
\r
283 FillHistograms(tpcTrack, stack, isMatch, 0);
\r
284 if(tpcTrack) delete tpcTrack;
\r
288 // ITS -> TPC prolongation efficiency
\r
290 for(Int_t iTrack=0; iTrack<entries;++iTrack)
\r
292 AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(iTrack);
\r
293 if(!track) continue;
\r
296 // ITS stand alone tracks
\r
298 if(!(track->GetStatus() & AliESDtrack::kITSpureSA)) continue;
\r
299 if(!(track->GetStatus() & AliESDtrack::kITSrefit)) continue;
\r
300 if(track->GetNcls(0)<4) continue;
\r
301 if(!track->HasPointOnITSLayer(0) && !track->HasPointOnITSLayer(1)) continue;
\r
303 // Check matching with TPC only track
\r
305 for(Int_t jTrack=0; jTrack<entries;++jTrack)
\r
307 if(iTrack==jTrack) continue;
\r
309 AliESDtrack *track2 = (AliESDtrack*)allChargedTracks->At(jTrack);
\r
310 if(!track2) continue;
\r
311 if(track2->Charge()==0) continue;
\r
312 if (!track2->GetTPCInnerParam()) continue;
\r
314 // Get TPC only tracks (must be deleted by user)
\r
315 AliESDtrack* tpcTrack2 = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, jTrack);
\r
316 if(!tpcTrack2) continue;
\r
318 // check loose cuts for TPC tracks
\r
319 //if(!esdTrackCuts->AcceptTrack(tpcTrack2)) { delete tpcTrack2; continue; }
\r
322 if (TMath::Abs(track->GetY() - tpcTrack2->GetY()) > 3) { delete tpcTrack2; continue; }
\r
323 if (TMath::Abs(track->GetSnp() - tpcTrack2->GetSnp()) > 0.2) { delete tpcTrack2; continue; }
\r
324 if (TMath::Abs(track->GetTgl() - tpcTrack2->GetTgl()) > 0.2) { delete tpcTrack2; continue; }
\r
334 Bool_t isITSTPC = kTRUE;
\r
335 FillHistograms(track, stack, isMatch, isITSTPC);
\r
339 if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
\r
343 //_____________________________________________________________________________
\r
344 void AlidNdPtEfficiency::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, const Bool_t isMatch, const Bool_t isITSTPC) const
\r
347 // Fill ESD track and MC histograms
\r
349 if(!esdTrack) return;
\r
350 Int_t charge = esdTrack->Charge();
\r
351 if(charge == 0.) return;
\r
353 Float_t pt = esdTrack->Pt();
\r
354 Float_t eta = esdTrack->Eta();
\r
355 Float_t phi = esdTrack->Phi();
\r
358 // Fill rec vs MC information
\r
360 Bool_t isPrim = kTRUE;
\r
362 if(IsUseMCInfo()) {
\r
364 Int_t label = esdTrack->GetLabel();
\r
365 if(label < 0.) return; // fake ITS track
\r
366 TParticle* particle = stack->Particle(label);
\r
367 if(!particle) return;
\r
368 if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;
\r
369 isPrim = stack->IsPhysicalPrimary(label);
\r
373 Double_t vRecMCTrackHist[6] = { eta,phi,pt,isPrim,charge,isMatch };
\r
376 fRecMCTrackHistITSTPC->Fill(vRecMCTrackHist);
\r
379 fRecMCTrackHistTPCITS->Fill(vRecMCTrackHist);
\r
383 //_____________________________________________________________________________
\r
384 Long64_t AlidNdPtEfficiency::Merge(TCollection* const list)
\r
386 // Merge list of objects (needed by PROOF)
\r
391 if (list->IsEmpty())
\r
394 TIterator* iter = list->MakeIterator();
\r
397 // collection of generated histograms
\r
399 while((obj = iter->Next()) != 0) {
\r
400 AlidNdPtEfficiency* entry = dynamic_cast<AlidNdPtEfficiency*>(obj);
\r
401 if (entry == 0) continue;
\r
404 fRecMCTrackHistTPCITS->Add(entry->fRecMCTrackHistTPCITS);
\r
405 fRecMCTrackHistITSTPC->Add(entry->fRecMCTrackHistITSTPC);
\r
413 //_____________________________________________________________________________
\r
414 void AlidNdPtEfficiency::Analyse()
\r
417 // Analyse histograms
\r
419 TH1::AddDirectory(kFALSE);
\r
420 TObjArray *aFolderObj = new TObjArray;
\r
429 AlidNdPtEventCuts *evtCuts = GetEventCuts();
\r
430 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
\r
431 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
\r
433 if(!evtCuts || !accCuts || !esdTrackCuts) {
\r
434 Error("AlidNdPtEfficiency::Analyse()", "cuts not available");
\r
439 // TPC->ITS efficiency
\r
443 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
\r
444 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
445 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(0);
\r
447 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_TPCITS");
\r
448 h1Dc->Divide(h1Dall);
\r
449 aFolderObj->Add(h1Dc);
\r
450 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
453 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
454 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
\r
455 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
456 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(1);
\r
458 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_TPCITS");
\r
459 h1Dc->Divide(h1Dall);
\r
460 aFolderObj->Add(h1Dc);
\r
461 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
464 fRecMCTrackHistTPCITS->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
465 h1Dall = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
\r
466 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(2,2);
\r
467 h1D = (TH1D *)fRecMCTrackHistTPCITS->Projection(2);
\r
469 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_TPCITS");
\r
470 h1Dc->Divide(h1Dall);
\r
471 aFolderObj->Add(h1Dc);
\r
472 fRecMCTrackHistTPCITS->GetAxis(5)->SetRange(1,2);
\r
476 // TPC->ITS efficiency
\r
479 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-1.5, 1.499);
\r
480 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
483 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
\r
484 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
485 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(0);
\r
487 h1Dc = (TH1D *)h1D->Clone("eff_vs_eta_ITSTPC");
\r
488 h1Dc->Divide(h1Dall);
\r
489 aFolderObj->Add(h1Dc);
\r
490 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
493 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
494 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
\r
495 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
496 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(1);
\r
498 h1Dc = (TH1D *)h1D->Clone("eff_vs_phi_ITSTPC");
\r
499 h1Dc->Divide(h1Dall);
\r
500 aFolderObj->Add(h1Dc);
\r
501 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
504 fRecMCTrackHistITSTPC->GetAxis(0)->SetRangeUser(-0.8, 0.799);
\r
505 h1Dall = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
\r
506 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(2,2);
\r
507 h1D = (TH1D *)fRecMCTrackHistITSTPC->Projection(2);
\r
509 h1Dc = (TH1D *)h1D->Clone("eff_vs_pT_ITSTPC");
\r
510 h1Dc->Divide(h1Dall);
\r
511 aFolderObj->Add(h1Dc);
\r
512 fRecMCTrackHistITSTPC->GetAxis(5)->SetRange(1,2);
\r
514 // export objects to analysis folder
\r
515 fAnalysisFolder = ExportToFolder(aFolderObj);
\r
517 // delete only TObjArray
\r
518 if(aFolderObj) delete aFolderObj;
\r
521 //_____________________________________________________________________________
\r
522 TFolder* AlidNdPtEfficiency::ExportToFolder(TObjArray * const array)
\r
524 // recreate folder avery time and export objects to new one
\r
526 AlidNdPtEfficiency * comp=this;
\r
527 TFolder *folder = comp->GetAnalysisFolder();
\r
529 TString name, title;
\r
530 TFolder *newFolder = 0;
\r
532 Int_t size = array->GetSize();
\r
535 // get name and title from old folder
\r
536 name = folder->GetName();
\r
537 title = folder->GetTitle();
\r
543 newFolder = CreateFolder(name.Data(),title.Data());
\r
544 newFolder->SetOwner();
\r
546 // add objects to folder
\r
548 newFolder->Add(array->At(i));
\r
556 //_____________________________________________________________________________
\r
557 TFolder* AlidNdPtEfficiency::CreateFolder(TString name,TString title) {
\r
558 // create folder for analysed histograms
\r
560 TFolder *folder = 0;
\r
561 folder = new TFolder(name.Data(),title.Data());
\r