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"),
77 fUsePhiWeight(kFALSE),
79 fUseRecentering(kFALSE),
80 fSaveTrackContribution(kFALSE),
110 // Default constructor
111 AliInfo("Event Plane Selection enabled.");
112 for(Int_t i = 0; i < 4; ++i) {
115 for(Int_t i = 0; i < 2; ++i) {
120 //________________________________________________________________________
121 AliEPSelectionTask::AliEPSelectionTask(const char *name):
122 AliAnalysisTaskSE(name),
123 fAnalysisInput("ESD"),
128 fUsePhiWeight(kFALSE),
129 fUsePtWeight(kFALSE),
130 fUseRecentering(kFALSE),
131 fSaveTrackContribution(kFALSE),
132 fUserphidist(kFALSE),
161 // Default constructor
162 AliInfo("Event Plane Selection enabled.");
163 DefineOutput(1, TList::Class());
164 for(Int_t i = 0; i < 4; i++) {
167 for(Int_t i = 0; i < 2; ++i) {
172 //________________________________________________________________________
173 AliEPSelectionTask::~AliEPSelectionTask()
176 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
181 delete fESDtrackCuts;
194 if (fPeriod.CompareTo("LHC11h")==0){
195 for(Int_t i = 0; i < 4; i++) {
201 if(fHruns) delete fHruns;
203 if(fQDist[0] && fQDist[1]) {
204 for(Int_t i = 0; i < 2; i++) {
213 //________________________________________________________________________
214 void AliEPSelectionTask::UserCreateOutputObjects()
216 // Create the output containers
217 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
218 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
220 fOutputList = new TList();
221 fOutputList->SetOwner();
222 fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
223 fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
224 fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
225 fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
226 fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
227 fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
228 fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
229 fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
231 fOutputList->Add(fHOutEventplaneQ);
232 fOutputList->Add(fHOutPhi);
233 fOutputList->Add(fHOutPhiCorr);
234 fOutputList->Add(fHOutsub1sub2);
235 fOutputList->Add(fHOutNTEPRes);
236 fOutputList->Add(fHOutPTPsi);
237 fOutputList->Add(fHOutleadPTPsi);
238 fOutputList->Add(fHOutDiff);
240 PostData(1, fOutputList);
245 //________________________________________________________________________
246 void AliEPSelectionTask::UserExec(Option_t */*option*/)
248 // Execute analysis for current event:
249 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
253 AliEventplane *esdEP;
256 Double_t fRP = 0.; // monte carlo reaction plane angle
258 if (fAnalysisInput.CompareTo("ESD")==0){
260 AliVEvent* event = InputEvent();
261 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
263 if (!(fRunNumber == esd->GetRunNumber())) {
264 fRunNumber = esd->GetRunNumber();
265 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
268 SetQvectorDist(); // recentring objects
273 AliMCEvent* mcEvent = MCEvent();
274 if (mcEvent && mcEvent->GenEventHeader()) {
275 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
276 if (headerH) fRP = headerH->ReactionPlaneAngle();
280 esdEP = esd->GetEventplane();
281 if (fSaveTrackContribution) {
282 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
283 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
284 esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
285 esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
286 esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
287 esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
290 TObjArray* tracklist = new TObjArray;
291 if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
292 if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
293 else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd);
294 const int nt = tracklist->GetEntries();
298 // qvector full event
299 fQVector = new TVector2(GetQ(esdEP,tracklist));
300 fEventplaneQ = fQVector->Phi()/2;
301 // q vector subevents
302 GetQsub(qq1, qq2, tracklist, esdEP);
303 fQsub1 = new TVector2(qq1);
304 fQsub2 = new TVector2(qq2);
305 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
307 esdEP->SetQVector(fQVector);
308 esdEP->SetEventplaneQ(fEventplaneQ);
309 esdEP->SetQsub(fQsub1,fQsub2);
310 esdEP->SetQsubRes(fQsubRes);
312 fHOutEventplaneQ->Fill(fEventplaneQ);
313 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
314 fHOutNTEPRes->Fill(nt,fQsubRes);
316 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
318 for (int iter = 0; iter<nt;iter++){
319 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
321 float delta = track->Phi()-fEventplaneQ;
322 while (delta < 0) delta += TMath::Pi();
323 while (delta > TMath::Pi()) delta -= TMath::Pi();
324 fHOutPTPsi->Fill(track->Pt(),delta);
325 fHOutPhi->Fill(track->Phi());
326 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
330 AliESDtrack* trmax = esd->GetTrack(0);
331 for (int iter = 1; iter<nt;iter++){
332 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
333 if (track && (track->Pt() > trmax->Pt())) trmax = track;
335 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
343 else if (fAnalysisInput.CompareTo("AOD")==0){
344 AliVEvent* event = InputEvent();
345 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
349 // get centrality of the event
350 AliAODHeader *header=dynamic_cast<AliAODHeader*>(aod->GetHeader());
351 if(!header) AliFatal("Not a standard AOD");
352 AliCentrality *centrality=header->GetCentralityP();
353 if(!centrality) AliError(Form("No AliCentrality attached to AOD header"));
354 fCentrality = centrality->GetCentralityPercentile("V0M");
356 if (!(fRunNumber == aod->GetRunNumber())) {
357 fRunNumber = aod->GetRunNumber();
358 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
365 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
366 if (headerH) fRP = headerH->GetReactionPlaneAngle();
369 esdEP = header->GetEventplaneP();
370 if(!esdEP) return; // protection against missing EP branch (nanoAODs)
374 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
376 if (fSaveTrackContribution) {
377 esdEP->GetQContributionXArray()->Set(maxID+1);
378 esdEP->GetQContributionYArray()->Set(maxID+1);
379 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
380 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
381 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
382 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
385 const int NT = tracklist->GetEntries();
389 // qvector full event
390 fQVector = new TVector2(GetQ(esdEP,tracklist));
391 fEventplaneQ = fQVector->Phi()/2;
392 // q vector subevents
393 GetQsub(qq1, qq2, tracklist, esdEP);
394 fQsub1 = new TVector2(qq1);
395 fQsub2 = new TVector2(qq2);
396 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
398 esdEP->SetQVector(fQVector);
399 esdEP->SetEventplaneQ(fEventplaneQ);
400 esdEP->SetQsub(fQsub1,fQsub2);
401 esdEP->SetQsubRes(fQsubRes);
403 fHOutEventplaneQ->Fill(fEventplaneQ);
404 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
405 fHOutNTEPRes->Fill(NT,fQsubRes);
407 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
409 for (int iter = 0; iter<NT;iter++){
410 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
412 float delta = track->Phi()-fEventplaneQ;
413 while (delta < 0) delta += TMath::Pi();
414 while (delta > TMath::Pi()) delta -= TMath::Pi();
415 fHOutPTPsi->Fill(track->Pt(),delta);
416 fHOutPhi->Fill(track->Phi());
417 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
421 AliAODTrack* trmax = dynamic_cast<AliAODTrack*>(aod->GetTrack(0));
422 if(!trmax) AliFatal("Not a standard AOD");
423 for (int iter = 1; iter<NT;iter++){
424 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
425 if (track && (track->Pt() > trmax->Pt())) trmax = track;
427 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
438 printf(" Analysis Input not known!\n\n ");
441 PostData(1, fOutputList);
444 //________________________________________________________________________
445 void AliEPSelectionTask::Terminate(Option_t */*option*/)
447 // Terminate analysis
450 //__________________________________________________________________________
451 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
459 // get recentering values
460 Double_t mean[2], rms[2];
464 int nt = tracklist->GetEntries();
466 for (int i=0; i<nt; i++){
468 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
470 weight = GetWeight(track);
471 if (fSaveTrackContribution){
472 idtemp = track->GetID();
473 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
474 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
475 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
477 mQx += (weight*cos(2*track->Phi())/rms[0]);
478 mQy += (weight*sin(2*track->Phi())/rms[1]);
481 mQ.Set(mQx-(mean[0]/rms[0]), mQy-(mean[1]/rms[1]));
485 //________________________________________________________________________
486 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
490 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
492 // get recentering values
493 Double_t mean[2], rms[2];
500 int nt = tracklist->GetEntries();
501 int trackcounter1=0, trackcounter2=0;
504 if (fSplitMethod == AliEPSelectionTask::kRandom){
506 for (Int_t i = 0; i < nt; i++) {
508 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
509 if (!track) continue;
510 weight = GetWeight(track);
511 idtemp = track->GetID();
512 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
514 // This loop splits the track set into 2 random subsets
515 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
516 float random = rn.Rndm();
518 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
519 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
520 if (fSaveTrackContribution){
521 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
522 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
527 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
528 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
529 if (fSaveTrackContribution){
530 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
531 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
536 else if( trackcounter1 >= int(nt/2.)){
537 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
538 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
539 if (fSaveTrackContribution){
540 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
541 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
546 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
547 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
548 if (fSaveTrackContribution){
549 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
550 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
555 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
557 for (Int_t i = 0; i < nt; i++) {
559 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
560 if (!track) continue;
561 weight = GetWeight(track);
562 Double_t eta = track->Eta();
563 idtemp = track->GetID();
564 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
566 if (eta > fEtaGap/2.) {
567 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
568 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
569 if (fSaveTrackContribution){
570 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
571 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
573 } else if (eta < -1.*fEtaGap/2.) {
574 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
575 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
576 if (fSaveTrackContribution){
577 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
578 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
582 } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
584 for (Int_t i = 0; i < nt; i++) {
586 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
587 if (!track) continue;
588 weight = GetWeight(track);
589 Short_t cha = track->Charge();
590 idtemp = track->GetID();
591 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
594 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
595 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
596 if (fSaveTrackContribution){
597 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
598 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
600 } else if (cha < 0) {
601 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
602 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
603 if (fSaveTrackContribution){
604 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
605 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
610 printf("plane resolution determination method not available!\n\n ");
613 // apply recenetering
614 mQ[0].Set(mQx1-(mean[0]/rms[0]), mQy1-(mean[1]/rms[1]));
615 mQ[1].Set(mQx2-(mean[0]/rms[0]), mQy2-(mean[1]/rms[1]));
620 //________________________________________________________________________
621 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
624 delete fESDtrackCuts;
627 if (fAnalysisInput.CompareTo("AOD")==0){
628 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
634 fESDtrackCuts = trackcuts;
637 //________________________________________________________________________
638 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
641 delete fESDtrackCuts;
644 if (fAnalysisInput.CompareTo("ESD")==0){
645 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
651 fESDtrackCuts = new AliESDtrackCuts();
652 fESDtrackCuts->SetPtRange(ptlow,ptup);
653 fESDtrackCuts->SetMinNClustersTPC(ntpc);
654 fESDtrackCuts->SetEtaRange(etalow,etaup);
655 fAODfilterbit = filterbit;
658 //_____________________________________________________________________________
660 void AliEPSelectionTask::SetTrackType(TString tracktype){
661 fTrackType = tracktype;
663 if (fTrackType.CompareTo("GLOBAL")==0){
664 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
667 if (fTrackType.CompareTo("TPC")==0){
668 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
671 fESDtrackCuts->SetPtRange(0.15,20.);
672 fESDtrackCuts->SetEtaRange(-0.8,0.8);
676 //________________________________________________________________________
677 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
680 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
681 if (fUsePtWeight && track) {
682 if (track->Pt()<2) ptweight=track->Pt();
685 return ptweight*GetPhiWeight(track);
688 //________________________________________________________________________
689 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
691 Double_t phiweight=1;
692 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
695 if(track) phiDist = SelectPhiDist(track);
697 if (fUsePhiWeight && phiDist && track) {
698 Double_t nParticles = phiDist->Integral();
699 Double_t nPhibins = phiDist->GetNbinsX();
701 Double_t Phi = track->Phi();
703 while (Phi<0) Phi += TMath::TwoPi();
704 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
706 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
708 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
713 //________________________________________________________________________
714 void AliEPSelectionTask::Recenter(Int_t var, Double_t * values)
717 if (fUseRecentering && fQDist[0] && fQDist[1] && fCentrality!=-1.) {
718 Int_t centbin = fQDist[0]->FindBin(fCentrality);
720 if(var==0) { // fill mean
721 values[0] = fQDist[0]->GetBinContent(centbin);
722 values[1] = fQDist[1]->GetBinContent(centbin);
724 else if(var==1) { // fill rms
725 values[0] = fQDist[0]->GetBinError(centbin);
726 values[1] = fQDist[1]->GetBinError(centbin);
727 // protection against division by zero
728 if(values[0]==0.0) values[0]=1.0;
729 if(values[1]==0.0) values[1]=1.0;
732 else { //default (no recentering)
739 //__________________________________________________________________________
740 void AliEPSelectionTask::SetPhiDist()
742 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
744 if (fPeriod.CompareTo("LHC10h")==0)
746 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
747 else if(fPeriod.CompareTo("LHC11h")==0){
748 Int_t runbin=fHruns->FindBin(fRunNumber);
749 if (fHruns->GetBinContent(runbin) > 1){
750 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
752 else if(fHruns->GetBinContent(runbin) < 2){
753 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
754 AliInfo("Using integrated Phi-weights for this run");
756 for (Int_t i = 0; i<4 ;i++)
763 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
764 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
766 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
767 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
769 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
770 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
772 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
773 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
774 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
775 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
776 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
777 fSparseDist->GetAxis(2)->SetRange(1,2);
779 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
782 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
787 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
794 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
795 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
800 cout << "empty bins - rebinning!" << endl;
801 fPhiDist[0]->Rebin();
808 AliError("After Maximum of rebinning still empty Phi-bins!!!");
811 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
812 AliInfo("No Phi-weights available. All Phi weights set to 1");
813 SetUsePhiWeight(kFALSE);
817 //__________________________________________________________________________
818 void AliEPSelectionTask::SetQvectorDist()
820 if(!fUseRecentering) return;
821 AliInfo(Form("Setting q vector distributions"));
822 fQDist[0] = (TProfile*) fQxContainer->GetObject(fRunNumber, "Default");
823 fQDist[1] = (TProfile*) fQyContainer->GetObject(fRunNumber, "Default");
825 if (!fQDist[0] || !fQDist[1]) {
826 AliError(Form("Cannot find OADB q-vector distributions for run %d. Using default values (mean=0,rms=1).", fRunNumber));
836 for (int i=1; i<=fQDist[0]->GetNbinsX()-2; i++){
837 if (!((fQDist[0]->GetBinEntries(i))>0)) {
842 cout << "empty bins - rebinning!" << endl;
851 AliError("After Maximum of rebinning still empty Qxy-bins!!!");
855 //__________________________________________________________________________
856 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
859 fUserphidist = kTRUE;
862 TObject* list = f.Get(listname);
863 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
864 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
869 //_________________________________________________________________________
870 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
872 TObjArray *acctracks = new TObjArray();
876 Int_t maxidtemp = -1;
881 fESDtrackCuts->GetPtRange(ptlow,ptup);
882 fESDtrackCuts->GetEtaRange(etalow,etaup);
883 Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
885 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
886 tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
887 if(!tr) AliFatal("Not a standard AOD");
888 maxidtemp = tr->GetID();
889 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
890 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
891 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
892 if (maxidtemp > maxid1) maxid1 = maxidtemp;
893 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
904 //_________________________________________________________________________
905 void AliEPSelectionTask::SetOADBandPeriod()
907 TString oadbfilename;
909 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
911 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
913 if (fAnalysisInput.CompareTo("AOD")==0){
914 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
915 } else if (fAnalysisInput.CompareTo("ESD")==0){
916 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
919 TFile foadb(oadbfilename);
920 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
922 AliInfo("Using Standard OADB");
923 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
924 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
929 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
932 // if it's already set and custom class is required, we use the one provided by the user
934 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
935 TFile *foadb = TFile::Open(oadbfilename);
936 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
938 AliInfo("Using Standard OADB");
939 fSparseDist = (THnSparse*) foadb->Get("Default");
940 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
943 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
944 fHruns->SetName("runsHisto");
948 if(fUseRecentering) {
949 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/eprecentering.root", AliAnalysisManager::GetOADBPath()));
950 TFile *foadb = TFile::Open(oadbfilename);
951 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
953 AliInfo("Using Standard OADB");
954 fQxContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qx");
955 fQyContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qy");
956 if (!fQxContainer || !fQyContainer) AliFatal("Cannot fetch OADB container for EP recentering");
963 //_________________________________________________________________________
964 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
966 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
967 else if(fPeriod.CompareTo("LHC11h")==0)
969 if (track->Charge() < 0)
971 if(track->Eta() < 0.) return fPhiDist[0];
972 else if (track->Eta() > 0.) return fPhiDist[2];
974 else if (track->Charge() > 0)
976 if(track->Eta() < 0.) return fPhiDist[1];
977 else if (track->Eta() > 0.) return fPhiDist[3];
984 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
986 // 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
987 TObjArray *acctracks = new TObjArray();
988 acctracks->SetOwner(kTRUE);
990 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
996 fESDtrackCuts->GetPtRange(ptlow,ptup);
997 fESDtrackCuts->GetEtaRange(etalow,etaup);
999 Double_t pDCA[3] = { 0. }; // momentum at DCA
1000 Double_t rDCA[3] = { 0. }; // position at DCA
1001 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
1002 Float_t cDCA[3] = {0.}; // covariance of impact parameters
1006 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
1008 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
1011 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
1013 // create a tpc only tracl
1014 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
1015 if(!track) continue;
1019 // only constrain tracks above threshold
1020 AliExternalTrackParam exParam;
1021 // take the B-field from the ESD, no 3D fieldMap available at this point
1022 Bool_t relate = false;
1023 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
1028 // fetch the track parameters at the DCA (unconstraint)
1029 if(track->GetTPCInnerParam()){
1030 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1031 track->GetTPCInnerParam()->GetXYZ(rDCA);
1033 // get the DCA to the vertex:
1034 track->GetImpactParametersTPC(dDCA,cDCA);
1035 // set the constrained parameters to the track
1036 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1040 Float_t eta = track->Eta();
1041 Float_t pT = track->Pt();
1043 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
1048 acctracks->Add(track);