1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*****************************************************
17 // Class AliEPSelectionTask
18 // Class to determine event plane
19 // author: Alberica Toia, Johanna Gramling
20 //*****************************************************
22 #include "AliEPSelectionTask.h"
28 #include <THnSparse.h>
31 #include <TObjString.h>
35 #include <TDirectory.h>
41 #include "AliAnalysisManager.h"
42 #include "AliVEvent.h"
44 #include "AliESDEvent.h"
45 #include "AliMCEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliESDHeader.h"
49 #include "AliESDInputHandler.h"
50 #include "AliAODHandler.h"
51 #include "AliAODEvent.h"
52 #include "AliHeader.h"
53 #include "AliGenHijingEventHeader.h"
54 #include "AliAnalysisTaskSE.h"
55 #include "AliPhysicsSelectionTask.h"
56 #include "AliPhysicsSelection.h"
57 #include "AliBackgroundSelection.h"
58 #include "AliESDUtils.h"
59 #include "AliOADBContainer.h"
60 #include "AliAODMCHeader.h"
61 #include "AliAODTrack.h"
62 #include "AliVTrack.h"
63 #include "AliEventplane.h"
67 ClassImp(AliEPSelectionTask)
69 //________________________________________________________________________
70 AliEPSelectionTask::AliEPSelectionTask():
72 fAnalysisInput("ESD"),
76 fUsePhiWeight(kFALSE),
78 fSaveTrackContribution(kFALSE),
106 // Default constructor
107 AliInfo("Event Plane Selection enabled.");
108 for(Int_t i = 0; i < 4; ++i) {
113 //________________________________________________________________________
114 AliEPSelectionTask::AliEPSelectionTask(const char *name):
115 AliAnalysisTaskSE(name),
116 fAnalysisInput("ESD"),
120 fUsePhiWeight(kFALSE),
121 fUsePtWeight(kFALSE),
122 fSaveTrackContribution(kFALSE),
123 fUserphidist(kFALSE),
150 // Default constructor
151 AliInfo("Event Plane Selection enabled.");
152 DefineOutput(1, TList::Class());
153 for(Int_t i = 0; i < 4; i++) {
158 //________________________________________________________________________
159 AliEPSelectionTask::~AliEPSelectionTask()
162 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
167 delete fESDtrackCuts;
180 if (fPeriod.CompareTo("LHC11h")==0){
181 for(Int_t i = 0; i < 4; i++) {
187 if(fHruns) delete fHruns;
191 //________________________________________________________________________
192 void AliEPSelectionTask::UserCreateOutputObjects()
194 // Create the output containers
195 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
196 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
198 fOutputList = new TList();
199 fOutputList->SetOwner();
200 fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
201 fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
202 fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
203 fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
204 fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
205 fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
206 fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
207 fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
209 fOutputList->Add(fHOutEventplaneQ);
210 fOutputList->Add(fHOutPhi);
211 fOutputList->Add(fHOutPhiCorr);
212 fOutputList->Add(fHOutsub1sub2);
213 fOutputList->Add(fHOutNTEPRes);
214 fOutputList->Add(fHOutPTPsi);
215 fOutputList->Add(fHOutleadPTPsi);
216 fOutputList->Add(fHOutDiff);
218 PostData(1, fOutputList);
223 //________________________________________________________________________
224 void AliEPSelectionTask::UserExec(Option_t */*option*/)
226 // Execute analysis for current event:
227 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
231 AliEventplane *esdEP;
234 Double_t fRP = 0.; // monte carlo reaction plane angle
236 if (fAnalysisInput.CompareTo("ESD")==0){
238 AliVEvent* event = InputEvent();
239 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
241 if (!(fRunNumber == esd->GetRunNumber())) {
242 fRunNumber = esd->GetRunNumber();
243 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
250 AliMCEvent* mcEvent = MCEvent();
251 if (mcEvent && mcEvent->GenEventHeader()) {
252 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
253 if (headerH) fRP = headerH->ReactionPlaneAngle();
257 esdEP = esd->GetEventplane();
258 if (fSaveTrackContribution) {
259 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
260 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
261 esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
262 esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
263 esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
264 esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
267 TObjArray* tracklist = new TObjArray;
268 if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
269 if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
270 else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd);
271 const int nt = tracklist->GetEntries();
274 fQVector = new TVector2(GetQ(esdEP,tracklist));
275 fEventplaneQ = fQVector->Phi()/2;
276 GetQsub(qq1, qq2, tracklist, esdEP);
277 fQsub1 = new TVector2(qq1);
278 fQsub2 = new TVector2(qq2);
279 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
281 esdEP->SetQVector(fQVector);
282 esdEP->SetEventplaneQ(fEventplaneQ);
283 esdEP->SetQsub(fQsub1,fQsub2);
284 esdEP->SetQsubRes(fQsubRes);
286 fHOutEventplaneQ->Fill(fEventplaneQ);
287 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
288 fHOutNTEPRes->Fill(nt,fQsubRes);
290 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
292 for (int iter = 0; iter<nt;iter++){
293 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
295 float delta = track->Phi()-fEventplaneQ;
296 while (delta < 0) delta += TMath::Pi();
297 while (delta > TMath::Pi()) delta -= TMath::Pi();
298 fHOutPTPsi->Fill(track->Pt(),delta);
299 fHOutPhi->Fill(track->Phi());
300 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
304 AliESDtrack* trmax = esd->GetTrack(0);
305 for (int iter = 1; iter<nt;iter++){
306 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
307 if (track && (track->Pt() > trmax->Pt())) trmax = track;
309 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
317 else if (fAnalysisInput.CompareTo("AOD")==0){
318 AliVEvent* event = InputEvent();
319 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
322 if (!(fRunNumber == aod->GetRunNumber())) {
323 fRunNumber = aod->GetRunNumber();
324 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
330 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
331 if (headerH) fRP = headerH->GetReactionPlaneAngle();
334 esdEP = aod->GetHeader()->GetEventplaneP();
335 if(!esdEP) return; // protection against missing EP branch (nanoAODs)
339 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
341 if (fSaveTrackContribution) {
342 esdEP->GetQContributionXArray()->Set(maxID+1);
343 esdEP->GetQContributionYArray()->Set(maxID+1);
344 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
345 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
346 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
347 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
350 const int NT = tracklist->GetEntries();
353 fQVector = new TVector2(GetQ(esdEP,tracklist));
354 fEventplaneQ = fQVector->Phi()/2.;
355 GetQsub(qq1, qq2, tracklist, esdEP);
356 fQsub1 = new TVector2(qq1);
357 fQsub2 = new TVector2(qq2);
358 fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
360 esdEP->SetQVector(fQVector);
361 esdEP->SetEventplaneQ(fEventplaneQ);
362 esdEP->SetQsub(fQsub1,fQsub2);
363 esdEP->SetQsubRes(fQsubRes);
365 fHOutEventplaneQ->Fill(fEventplaneQ);
366 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
367 fHOutNTEPRes->Fill(NT,fQsubRes);
369 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
371 for (int iter = 0; iter<NT;iter++){
372 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
374 float delta = track->Phi()-fEventplaneQ;
375 while (delta < 0) delta += TMath::Pi();
376 while (delta > TMath::Pi()) delta -= TMath::Pi();
377 fHOutPTPsi->Fill(track->Pt(),delta);
378 fHOutPhi->Fill(track->Phi());
379 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
383 AliAODTrack* trmax = aod->GetTrack(0);
384 for (int iter = 1; iter<NT;iter++){
385 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
386 if (track && (track->Pt() > trmax->Pt())) trmax = track;
388 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
399 printf(" Analysis Input not known!\n\n ");
402 PostData(1, fOutputList);
405 //________________________________________________________________________
406 void AliEPSelectionTask::Terminate(Option_t */*option*/)
408 // Terminate analysis
411 //__________________________________________________________________________
412 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
421 int nt = tracklist->GetEntries();
423 for (int i=0; i<nt; i++){
425 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
427 weight = GetWeight(track);
428 if (fSaveTrackContribution){
429 idtemp = track->GetID();
430 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
431 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
432 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
434 mQx += (weight*cos(2*track->Phi()));
435 mQy += (weight*sin(2*track->Phi()));
442 //________________________________________________________________________
443 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
447 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
453 int nt = tracklist->GetEntries();
454 int trackcounter1=0, trackcounter2=0;
457 if (fSplitMethod == AliEPSelectionTask::kRandom){
459 for (Int_t i = 0; i < nt; i++) {
461 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
462 if (!track) continue;
463 weight = GetWeight(track);
464 idtemp = track->GetID();
465 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
467 // This loop splits the track set into 2 random subsets
468 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
469 float random = rn.Rndm();
471 mQx1 += (weight*cos(2*track->Phi()));
472 mQy1 += (weight*sin(2*track->Phi()));
473 if (fSaveTrackContribution){
474 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
475 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
480 mQx2 += (weight*cos(2*track->Phi()));
481 mQy2 += (weight*sin(2*track->Phi()));
482 if (fSaveTrackContribution){
483 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
484 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
489 else if( trackcounter1 >= int(nt/2.)){
490 mQx2 += (weight*cos(2*track->Phi()));
491 mQy2 += (weight*sin(2*track->Phi()));
492 if (fSaveTrackContribution){
493 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
494 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
499 mQx1 += (weight*cos(2*track->Phi()));
500 mQy1 += (weight*sin(2*track->Phi()));
501 if (fSaveTrackContribution){
502 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
503 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
508 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
510 for (Int_t i = 0; i < nt; i++) {
512 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
513 if (!track) continue;
514 weight = GetWeight(track);
515 Double_t eta = track->Eta();
516 idtemp = track->GetID();
517 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
519 if (eta > fEtaGap/2.) {
520 mQx1 += (weight*cos(2*track->Phi()));
521 mQy1 += (weight*sin(2*track->Phi()));
522 if (fSaveTrackContribution){
523 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
524 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
526 } else if (eta < -1.*fEtaGap/2.) {
527 mQx2 += (weight*cos(2*track->Phi()));
528 mQy2 += (weight*sin(2*track->Phi()));
529 if (fSaveTrackContribution){
530 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
531 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
535 } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
537 for (Int_t i = 0; i < nt; i++) {
539 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
540 if (!track) continue;
541 weight = GetWeight(track);
542 Short_t cha = track->Charge();
543 idtemp = track->GetID();
544 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
547 mQx1 += (weight*cos(2*track->Phi()));
548 mQy1 += (weight*sin(2*track->Phi()));
549 if (fSaveTrackContribution){
550 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
551 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
553 } else if (cha < 0) {
554 mQx2 += (weight*cos(2*track->Phi()));
555 mQy2 += (weight*sin(2*track->Phi()));
556 if (fSaveTrackContribution){
557 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
558 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
563 printf("plane resolution determination method not available!\n\n ");
567 mQ[0].Set(mQx1,mQy1);
568 mQ[1].Set(mQx2,mQy2);
573 //________________________________________________________________________
574 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
577 delete fESDtrackCuts;
580 if (fAnalysisInput.CompareTo("AOD")==0){
581 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
587 fESDtrackCuts = trackcuts;
590 //________________________________________________________________________
591 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
594 delete fESDtrackCuts;
597 if (fAnalysisInput.CompareTo("ESD")==0){
598 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
604 fESDtrackCuts = new AliESDtrackCuts();
605 fESDtrackCuts->SetPtRange(ptlow,ptup);
606 fESDtrackCuts->SetMinNClustersTPC(ntpc);
607 fESDtrackCuts->SetEtaRange(etalow,etaup);
608 fAODfilterbit = filterbit;
611 //_____________________________________________________________________________
613 void AliEPSelectionTask::SetTrackType(TString tracktype){
614 fTrackType = tracktype;
616 if (fTrackType.CompareTo("GLOBAL")==0){
617 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
620 if (fTrackType.CompareTo("TPC")==0){
621 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
624 fESDtrackCuts->SetPtRange(0.15,20.);
625 fESDtrackCuts->SetEtaRange(-0.8,0.8);
629 //________________________________________________________________________
630 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
633 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
634 if (fUsePtWeight && track) {
635 if (track->Pt()<2) ptweight=track->Pt();
638 return ptweight*GetPhiWeight(track);
641 //________________________________________________________________________
642 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
644 Double_t phiweight=1;
645 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
648 if(track) phiDist = SelectPhiDist(track);
650 if (fUsePhiWeight && phiDist && track) {
651 Double_t nParticles = phiDist->Integral();
652 Double_t nPhibins = phiDist->GetNbinsX();
654 Double_t Phi = track->Phi();
656 while (Phi<0) Phi += TMath::TwoPi();
657 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
659 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
661 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
666 //__________________________________________________________________________
667 void AliEPSelectionTask::SetPhiDist()
669 if(!fUserphidist && (fPeriod.CompareTo("LHC10h") == 0 || fPeriod.CompareTo("LHC11h") == 0)) { // if it's already set and custom class is required, we use the one provided by the user
671 if (fPeriod.CompareTo("LHC10h")==0)
673 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
674 else if(fPeriod.CompareTo("LHC11h")==0){
675 Int_t runbin=fHruns->FindBin(fRunNumber);
676 if (fHruns->GetBinContent(runbin) > 1){
677 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
679 else if(fHruns->GetBinContent(runbin) < 2){
680 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
681 AliInfo("Using integrated Phi-weights for this run");
683 for (Int_t i = 0; i<4 ;i++)
690 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
691 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
693 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
694 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
696 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
697 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
699 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
700 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
701 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
702 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
703 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
704 fSparseDist->GetAxis(2)->SetRange(1,2);
706 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
709 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
714 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
721 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
722 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
727 cout << "empty bins - rebinning!" << endl;
728 fPhiDist[0]->Rebin();
735 AliError("After Maximum of rebinning still empty Phi-bins!!!");
738 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
739 AliInfo("No Phi-weights available. All Phi weights set to 1");
740 SetUsePhiWeight(kFALSE);
744 //__________________________________________________________________________
745 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
748 fUserphidist = kTRUE;
751 TObject* list = f.Get(listname);
752 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
753 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
759 //_________________________________________________________________________
760 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
762 TObjArray *acctracks = new TObjArray();
766 Int_t maxidtemp = -1;
771 fESDtrackCuts->GetPtRange(ptlow,ptup);
772 fESDtrackCuts->GetEtaRange(etalow,etaup);
773 Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
775 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
776 tr = aod->GetTrack(i);
777 maxidtemp = tr->GetID();
778 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
779 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
780 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
781 if (maxidtemp > maxid1) maxid1 = maxidtemp;
782 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
793 //_________________________________________________________________________
794 void AliEPSelectionTask::SetOADBandPeriod()
796 TString oadbfilename;
798 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
800 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
802 if (fAnalysisInput.CompareTo("AOD")==0){
803 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
804 } else if (fAnalysisInput.CompareTo("ESD")==0){
805 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
808 TFile foadb(oadbfilename);
809 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
811 AliInfo("Using Standard OADB");
812 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
813 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
818 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
821 // if it's already set and custom class is required, we use the one provided by the user
823 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
824 TFile *foadb = TFile::Open(oadbfilename);
825 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
827 AliInfo("Using Standard OADB");
828 fSparseDist = (THnSparse*) foadb->Get("Default");
829 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
832 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
833 fHruns->SetName("runsHisto");
839 //_________________________________________________________________________
840 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
842 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
843 else if(fPeriod.CompareTo("LHC11h")==0)
845 if (track->Charge() < 0)
847 if(track->Eta() < 0.) return fPhiDist[0];
848 else if (track->Eta() > 0.) return fPhiDist[2];
850 else if (track->Charge() > 0)
852 if(track->Eta() < 0.) return fPhiDist[1];
853 else if (track->Eta() > 0.) return fPhiDist[3];
860 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
862 // Need to do this hack beacuse only this type of TPC only tracks in AOD is available and one type of Phi-weights is used
863 TObjArray *acctracks = new TObjArray();
864 acctracks->SetOwner(kTRUE);
866 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
872 fESDtrackCuts->GetPtRange(ptlow,ptup);
873 fESDtrackCuts->GetEtaRange(etalow,etaup);
875 Double_t pDCA[3] = { 0. }; // momentum at DCA
876 Double_t rDCA[3] = { 0. }; // position at DCA
877 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
878 Float_t cDCA[3] = {0.}; // covariance of impact parameters
882 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
884 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
887 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
889 // create a tpc only tracl
890 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
895 // only constrain tracks above threshold
896 AliExternalTrackParam exParam;
897 // take the B-field from the ESD, no 3D fieldMap available at this point
898 Bool_t relate = false;
899 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
904 // fetch the track parameters at the DCA (unconstraint)
905 if(track->GetTPCInnerParam()){
906 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
907 track->GetTPCInnerParam()->GetXYZ(rDCA);
909 // get the DCA to the vertex:
910 track->GetImpactParametersTPC(dDCA,cDCA);
911 // set the constrained parameters to the track
912 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
916 Float_t eta = track->Eta();
917 Float_t pT = track->Pt();
919 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
924 acctracks->Add(track);