X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FAliEPSelectionTask.cxx;h=7d0d235c61d48c8aa98adea8557926064ab37dbd;hb=4a8b470ced4e27715bcd99e6b957fc08529b6418;hp=d04ef61c0a58bdacd37c9bb64d9759e688531b42;hpb=90267ca6770aecc93227a97e4a9c634d9f034632;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliEPSelectionTask.cxx b/ANALYSIS/AliEPSelectionTask.cxx index d04ef61c0a5..7d0d235c61d 100644 --- a/ANALYSIS/AliEPSelectionTask.cxx +++ b/ANALYSIS/AliEPSelectionTask.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -56,28 +57,35 @@ #include "AliBackgroundSelection.h" #include "AliESDUtils.h" #include "AliOADBContainer.h" - +#include "AliAODMCHeader.h" +#include "AliAODTrack.h" +#include "AliVTrack.h" #include "AliEventplane.h" +using std::cout; +using std::endl; ClassImp(AliEPSelectionTask) //________________________________________________________________________ AliEPSelectionTask::AliEPSelectionTask(): AliAnalysisTaskSE(), - fDebug(0), fAnalysisInput("ESD"), fTrackType("TPC"), + fPeriod(""), fUseMCRP(kFALSE), fUsePhiWeight(kFALSE), fUsePtWeight(kFALSE), fSaveTrackContribution(kFALSE), - fuserphidist(kFALSE), - fusercuts(kFALSE), - frunNumber(-15), + fUserphidist(kFALSE), + fUsercuts(kFALSE), + fRunNumber(-15), + fAODfilterbit(1), + fEtaGap(0.), + fSplitMethod(0), fESDtrackCuts(0), - ftracklist(0), fEPContainer(0), - fPhiDist(0), + fSparseDist(0), + fHruns(0), fQVector(0), fQContributionX(0), fQContributionY(0), @@ -97,25 +105,31 @@ AliAnalysisTaskSE(), { // Default constructor AliInfo("Event Plane Selection enabled."); + for(Int_t i = 0; i < 4; ++i) { + fPhiDist[i] = 0; + } } //________________________________________________________________________ AliEPSelectionTask::AliEPSelectionTask(const char *name): AliAnalysisTaskSE(name), - fDebug(0), fAnalysisInput("ESD"), fTrackType("TPC"), + fPeriod(""), fUseMCRP(kFALSE), fUsePhiWeight(kFALSE), fUsePtWeight(kFALSE), fSaveTrackContribution(kFALSE), - fuserphidist(kFALSE), - fusercuts(kFALSE), - frunNumber(-15), + fUserphidist(kFALSE), + fUsercuts(kFALSE), + fRunNumber(-15), + fAODfilterbit(1), + fEtaGap(0.), + fSplitMethod(0), fESDtrackCuts(0), - ftracklist(0), fEPContainer(0), - fPhiDist(0), + fSparseDist(0), + fHruns(0), fQVector(0), fQContributionX(0), fQContributionY(0), @@ -136,6 +150,9 @@ AliEPSelectionTask::AliEPSelectionTask(const char *name): // Default constructor AliInfo("Event Plane Selection enabled."); DefineOutput(1, TList::Class()); + for(Int_t i = 0; i < 4; i++) { + fPhiDist[i] = 0; + } } //________________________________________________________________________ @@ -150,30 +167,25 @@ AliEPSelectionTask::~AliEPSelectionTask() delete fESDtrackCuts; fESDtrackCuts = 0; } - if (fQVector){ - delete fQVector; - fQVector = 0; - } - if (fQsub1){ - delete fQsub1; - fQsub1 = 0; - } - if (fQsub2){ - delete fQsub2; - fQsub2 = 0; - } - if (fPhiDist){ - delete fPhiDist; - fPhiDist = 0; - } - if (ftracklist){ - delete ftracklist; - ftracklist = 0; + if (fUserphidist) { + if (fPhiDist[0]) { + delete fPhiDist[0]; + fPhiDist[0] = 0; + } } if (fEPContainer){ delete fEPContainer; fEPContainer = 0; } + if (fPeriod.CompareTo("LHC11h")==0){ + for(Int_t i = 0; i < 4; i++) { + if(fPhiDist[i]){ + delete fPhiDist[i]; + fPhiDist[i] = 0; + } + } + if(fHruns) delete fHruns; + } } //________________________________________________________________________ @@ -204,19 +216,8 @@ void AliEPSelectionTask::UserCreateOutputObjects() fOutputList->Add(fHOutDiff); PostData(1, fOutputList); - - if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user - TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath())); - TFile foadb(oadbfilename); - if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data())); - - AliInfo("Using Standard OADB"); - fEPContainer = (AliOADBContainer*) foadb.Get("epphidist"); - if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection"); - foadb.Close(); - } } //________________________________________________________________________ @@ -225,47 +226,56 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) // Execute analysis for current event: if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n"); -// frunNumber = -15; +// fRunNumber = -15; - AliEventplane* esdEP = 0; - TVector2 QQ1; - TVector2 QQ2; - Double_t fRP = 0.; // the monte carlo reaction plane angle + AliEventplane *esdEP; + TVector2 qq1; + TVector2 qq2; + Double_t fRP = 0.; // monte carlo reaction plane angle if (fAnalysisInput.CompareTo("ESD")==0){ - + AliVEvent* event = InputEvent(); AliESDEvent* esd = dynamic_cast(event); - if (!(frunNumber == esd->GetRunNumber())) { - frunNumber = esd->GetRunNumber(); - SetPhiDist(); - } - - - if (fUseMCRP) { - AliMCEvent* mcEvent = MCEvent(); - if (mcEvent && mcEvent->GenEventHeader()) { - AliGenHijingEventHeader* headerH = dynamic_cast(mcEvent->GenEventHeader()); - if (headerH) fRP = headerH->ReactionPlaneAngle(); + if (esd){ + if (!(fRunNumber == esd->GetRunNumber())) { + fRunNumber = esd->GetRunNumber(); + AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber)); + SetOADBandPeriod(); + SetPhiDist(); } - } - if (esd){ + + + if (fUseMCRP) { + AliMCEvent* mcEvent = MCEvent(); + if (mcEvent && mcEvent->GenEventHeader()) { + AliGenHijingEventHeader* headerH = dynamic_cast(mcEvent->GenEventHeader()); + if (headerH) fRP = headerH->ReactionPlaneAngle(); + } + } + esdEP = esd->GetEventplane(); if (fSaveTrackContribution) { esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks()); esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks()); + esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks()); + esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks()); + esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks()); + esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks()); } - if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE); - if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE); - const int NT = ftracklist->GetEntries(); + TObjArray* tracklist = new TObjArray; + if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE); + if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE); + else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd); + const int nt = tracklist->GetEntries(); - if (NT>4){ - fQVector = new TVector2(GetQ(esdEP,ftracklist)); + if (nt>4){ + fQVector = new TVector2(GetQ(esdEP,tracklist)); fEventplaneQ = fQVector->Phi()/2; - GetQsub(QQ1, QQ2, ftracklist); - fQsub1 = new TVector2(QQ1); - fQsub2 = new TVector2(QQ2); + GetQsub(qq1, qq2, tracklist, esdEP); + fQsub1 = new TVector2(qq1); + fQsub2 = new TVector2(qq2); fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2); esdEP->SetQVector(fQVector); @@ -275,39 +285,116 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) fHOutEventplaneQ->Fill(fEventplaneQ); fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2); - fHOutNTEPRes->Fill(NT,fQsubRes); + fHOutNTEPRes->Fill(nt,fQsubRes); if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP); - for (int iter = 0; iter (ftracklist->At(iter)); + for (int iter = 0; iter (tracklist->At(iter)); if (track) { float delta = track->Phi()-fEventplaneQ; while (delta < 0) delta += TMath::Pi(); while (delta > TMath::Pi()) delta -= TMath::Pi(); fHOutPTPsi->Fill(track->Pt(),delta); fHOutPhi->Fill(track->Phi()); - fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track)); - } + fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track)); + } } AliESDtrack* trmax = esd->GetTrack(0); - for (int iter = 1; iter (ftracklist->At(iter)); + for (int iter = 1; iter (tracklist->At(iter)); if (track && (track->Pt() > trmax->Pt())) trmax = track; } fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ); } - delete ftracklist; + tracklist->Clear(); + delete tracklist; + tracklist = 0; } } - else if (fAnalysisInput.CompareTo("AOD")==0){ - //AliAODEvent *aod = dynamic_cast (InputEvent()); - // to be implemented - printf(" AOD analysis not yet implemented!!!\n\n"); - return; + else if (fAnalysisInput.CompareTo("AOD")==0){ + AliVEvent* event = InputEvent(); + AliAODEvent* aod = dynamic_cast(event); + + if (aod){ + if (!(fRunNumber == aod->GetRunNumber())) { + fRunNumber = aod->GetRunNumber(); + AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber)); + SetOADBandPeriod(); + SetPhiDist(); + } + + if (fUseMCRP) { + AliAODMCHeader *headerH = dynamic_cast(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + if (headerH) fRP = headerH->GetReactionPlaneAngle(); + } + + esdEP = aod->GetHeader()->GetEventplaneP(); + if(!esdEP) return; // protection against missing EP branch (nanoAODs) + esdEP->Reset(); + + Int_t maxID = 0; + TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID); + + if (fSaveTrackContribution) { + esdEP->GetQContributionXArray()->Set(maxID+1); + esdEP->GetQContributionYArray()->Set(maxID+1); + esdEP->GetQContributionXArraysub1()->Set(maxID+1); + esdEP->GetQContributionYArraysub1()->Set(maxID+1); + esdEP->GetQContributionXArraysub2()->Set(maxID+1); + esdEP->GetQContributionYArraysub2()->Set(maxID+1); + } + + const int NT = tracklist->GetEntries(); + + if (NT>4){ + fQVector = new TVector2(GetQ(esdEP,tracklist)); + fEventplaneQ = fQVector->Phi()/2.; + GetQsub(qq1, qq2, tracklist, esdEP); + fQsub1 = new TVector2(qq1); + fQsub2 = new TVector2(qq2); + fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.); + + esdEP->SetQVector(fQVector); + esdEP->SetEventplaneQ(fEventplaneQ); + esdEP->SetQsub(fQsub1,fQsub2); + esdEP->SetQsubRes(fQsubRes); + + fHOutEventplaneQ->Fill(fEventplaneQ); + fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.); + fHOutNTEPRes->Fill(NT,fQsubRes); + + if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP); + + for (int iter = 0; iter (tracklist->At(iter)); + if (track) { + float delta = track->Phi()-fEventplaneQ; + while (delta < 0) delta += TMath::Pi(); + while (delta > TMath::Pi()) delta -= TMath::Pi(); + fHOutPTPsi->Fill(track->Pt(),delta); + fHOutPhi->Fill(track->Phi()); + fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track)); + } + } + + AliAODTrack* trmax = aod->GetTrack(0); + for (int iter = 1; iter (tracklist->At(iter)); + if (track && (track->Pt() > trmax->Pt())) trmax = track; + } + fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ); + } + delete tracklist; + tracklist = 0; + } + + } + + else { printf(" Analysis Input not known!\n\n "); return; @@ -324,24 +411,28 @@ void AliEPSelectionTask::Terminate(Option_t */*option*/) //__________________________________________________________________________ TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist) { +// Get the Q vector TVector2 mQ; float mQx=0, mQy=0; - AliESDtrack* track; + AliVTrack* track; Double_t weight; + Int_t idtemp = -1; - int NT = tracklist->GetEntries(); + int nt = tracklist->GetEntries(); - for (int i=0; i (tracklist->At(i)); + track = dynamic_cast (tracklist->At(i)); if (track) { weight = GetWeight(track); - if (fSaveTrackContribution){ - EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID()); - EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID()); - } - mQx += (weight*cos(2*track->Phi())); - mQy += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + idtemp = track->GetID(); + if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1; + EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + mQx += (weight*cos(2*track->Phi())); + mQy += (weight*sin(2*track->Phi())); } } mQ.Set(mQx,mQy); @@ -349,49 +440,130 @@ TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist) } //________________________________________________________________________ -void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist) +void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP) { +// Get Qsub TVector2 mQ[2]; float mQx1=0, mQy1=0, mQx2=0, mQy2=0; Double_t weight; - AliESDtrack* track; - TRandom2 RN = 0; + AliVTrack* track; + TRandom2 rn = 0; - int NT = tracklist->GetEntries(); + int nt = tracklist->GetEntries(); int trackcounter1=0, trackcounter2=0; - - for (Int_t i = 0; i < NT; i++) { - weight = 1; - track = dynamic_cast (tracklist->At(i)); - if (!track) continue; - weight = GetWeight(track); + int idtemp = 0; + + if (fSplitMethod == AliEPSelectionTask::kRandom){ + + for (Int_t i = 0; i < nt; i++) { + weight = 1; + track = dynamic_cast (tracklist->At(i)); + if (!track) continue; + weight = GetWeight(track); + idtemp = track->GetID(); + if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1; - // This loop splits the track set into 2 random subsets - if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){ - float random = RN.Rndm(); - if(random < .5){ + // This loop splits the track set into 2 random subsets + if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){ + float random = rn.Rndm(); + if(random < .5){ + mQx1 += (weight*cos(2*track->Phi())); + mQy1 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + trackcounter1++; + } + else { + mQx2 += (weight*cos(2*track->Phi())); + mQy2 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + trackcounter2++; + } + } + else if( trackcounter1 >= int(nt/2.)){ + mQx2 += (weight*cos(2*track->Phi())); + mQy2 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + trackcounter2++; + } + else { mQx1 += (weight*cos(2*track->Phi())); mQy1 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp); + } trackcounter1++; } - else { + } + } else if (fSplitMethod == AliEPSelectionTask::kEta) { + + for (Int_t i = 0; i < nt; i++) { + weight = 1; + track = dynamic_cast (tracklist->At(i)); + if (!track) continue; + weight = GetWeight(track); + Double_t eta = track->Eta(); + idtemp = track->GetID(); + if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1; + + if (eta > fEtaGap/2.) { + mQx1 += (weight*cos(2*track->Phi())); + mQy1 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + } else if (eta < -1.*fEtaGap/2.) { mQx2 += (weight*cos(2*track->Phi())); mQy2 += (weight*sin(2*track->Phi())); - trackcounter2++; + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp); + } } } - else if( trackcounter1 >= int(NT/2.)){ - mQx2 += (weight*cos(2*track->Phi())); - mQy2 += (weight*sin(2*track->Phi())); - trackcounter2++; - } - else { - mQx1 += (weight*cos(2*track->Phi())); - mQy1 += (weight*sin(2*track->Phi())); - trackcounter1++; + } else if (fSplitMethod == AliEPSelectionTask::kCharge) { + + for (Int_t i = 0; i < nt; i++) { + weight = 1; + track = dynamic_cast (tracklist->At(i)); + if (!track) continue; + weight = GetWeight(track); + Short_t cha = track->Charge(); + idtemp = track->GetID(); + if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1; + + if (cha > 0) { + mQx1 += (weight*cos(2*track->Phi())); + mQy1 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + } else if (cha < 0) { + mQx2 += (weight*cos(2*track->Phi())); + mQy2 += (weight*sin(2*track->Phi())); + if (fSaveTrackContribution){ + EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp); + EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp); + } + } } + } else { + printf("plane resolution determination method not available!\n\n "); + return; } + mQ[0].Set(mQx1,mQy1); mQ[1].Set(mQx2,mQy2); Q1 = mQ[0]; @@ -401,27 +573,65 @@ void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklis //________________________________________________________________________ void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){ - fusercuts = kTRUE; + if(fESDtrackCuts){ + delete fESDtrackCuts; + fESDtrackCuts = 0; + } + if (fAnalysisInput.CompareTo("AOD")==0){ + AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts"); + fUsercuts = kFALSE; + SetTrackType("TPC"); + return; + } + fUsercuts = kTRUE; fESDtrackCuts = trackcuts; } +//________________________________________________________________________ +void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){ + + if(fESDtrackCuts){ + delete fESDtrackCuts; + fESDtrackCuts = 0; + } + if (fAnalysisInput.CompareTo("ESD")==0){ + AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts"); + fUsercuts = kFALSE; + SetTrackType("TPC"); + return; + } + fUsercuts = kTRUE; + fESDtrackCuts = new AliESDtrackCuts(); + fESDtrackCuts->SetPtRange(ptlow,ptup); + fESDtrackCuts->SetMinNClustersTPC(ntpc); + fESDtrackCuts->SetEtaRange(etalow,etaup); + fAODfilterbit = filterbit; +} + //_____________________________________________________________________________ + void AliEPSelectionTask::SetTrackType(TString tracktype){ fTrackType = tracktype; - if (!fusercuts) { - if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); - if (fTrackType.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); - fESDtrackCuts->SetPtRange(0.15,20); + if (!fUsercuts) { + if (fTrackType.CompareTo("GLOBAL")==0){ + fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); + fAODfilterbit = 32; + } + if (fTrackType.CompareTo("TPC")==0){ + fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); + fAODfilterbit = 128; + } + fESDtrackCuts->SetPtRange(0.15,20.); fESDtrackCuts->SetEtaRange(-0.8,0.8); } } //________________________________________________________________________ -Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track) +Double_t AliEPSelectionTask::GetWeight(TObject* track1) { Double_t ptweight=1; - - if (fUsePtWeight) { + AliVTrack* track = dynamic_cast(track1); + if (fUsePtWeight && track) { if (track->Pt()<2) ptweight=track->Pt(); else ptweight=2; } @@ -429,20 +639,24 @@ Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track) } //________________________________________________________________________ -Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track) +Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1) { Double_t phiweight=1; + AliVTrack* track = dynamic_cast(track1); + + TH1F *phiDist = 0x0; + if(track) phiDist = SelectPhiDist(track); - if (fUsePhiWeight && fPhiDist && track) { - Double_t nParticles = fPhiDist->Integral(); - Double_t nPhibins = fPhiDist->GetNbinsX(); + if (fUsePhiWeight && phiDist && track) { + Double_t nParticles = phiDist->Integral(); + Double_t nPhibins = phiDist->GetNbinsX(); Double_t Phi = track->Phi(); while (Phi<0) Phi += TMath::TwoPi(); while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi(); - Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi())); + Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi())); if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue; } @@ -452,50 +666,266 @@ Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track) //__________________________________________________________________________ void AliEPSelectionTask::SetPhiDist() { - if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user + 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 + + if (fPeriod.CompareTo("LHC10h")==0) + { + fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");} + else if(fPeriod.CompareTo("LHC11h")==0){ + Int_t runbin=fHruns->FindBin(fRunNumber); + if (fHruns->GetBinContent(runbin) > 1){ + fSparseDist->GetAxis(0)->SetRange(runbin,runbin); + } + else if(fHruns->GetBinContent(runbin) < 2){ + fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights + AliInfo("Using integrated Phi-weights for this run"); + } + for (Int_t i = 0; i<4 ;i++) + { + if(fPhiDist[i]){ + delete fPhiDist[i]; + fPhiDist[i] = 0x0; + } + if(i == 0){ + fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge + fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta + if(i == 1){ + fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge + fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta + if(i == 2){ + fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge + fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta + if(i == 3){ + fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge + fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta + fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi + fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber)); + fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes + fSparseDist->GetAxis(2)->SetRange(1,2); + } + fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis + } - fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default"); - if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber)); + if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber)); } - else { - AliInfo("Using Custom Phi Distribution"); - } + - Bool_t emptybins; + if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){ + Bool_t emptybins; - int iter = 0; - while (iter<3){ - emptybins = kFALSE; + int iter = 0; + while (iter<3){ + emptybins = kFALSE; - for (int i=1; iGetNbinsX(); i++){ - if (!((fPhiDist->GetBinContent(i))>0)) { - emptybins = kTRUE; - } - } - if (emptybins) { - cout << "empty bins - rebinning!" << endl; - fPhiDist->Rebin(); - iter++; - } - else iter = 3; - } + for (int i=1; iGetNbinsX(); i++){ + if (!((fPhiDist[0]->GetBinContent(i))>0)) { + emptybins = kTRUE; + } + } + if (emptybins) { + cout << "empty bins - rebinning!" << endl; + fPhiDist[0]->Rebin(); + iter++; + } + else iter = 3; + } - if (emptybins) { - AliError("After Maximum of rebinning still empty Phi-bins!!!"); + if (emptybins) { + AliError("After Maximum of rebinning still empty Phi-bins!!!"); + } + } + if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){ + AliInfo("No Phi-weights available. All Phi weights set to 1"); + SetUsePhiWeight(kFALSE); } } //__________________________________________________________________________ -void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname) +void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname) { - fuserphidist = kTRUE; + fUserphidist = kTRUE; TFile f(infilename); TObject* list = f.Get(listname); - fPhiDist = (TH1F*)list->FindObject("fHOutPhi"); - if (!fPhiDist) AliFatal("Phi Distribution not found!!!"); + fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi"); + if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!"); f.Close(); } + + +//_________________________________________________________________________ +TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid) +{ + TObjArray *acctracks = new TObjArray(); + + AliAODTrack *tr = 0; + Int_t maxid1 = 0; + Int_t maxidtemp = -1; + Float_t ptlow = 0; + Float_t ptup = 0; + Float_t etalow = 0; + Float_t etaup = 0; + fESDtrackCuts->GetPtRange(ptlow,ptup); + fESDtrackCuts->GetEtaRange(etalow,etaup); + Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC(); + + for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){ + tr = aod->GetTrack(i); + maxidtemp = tr->GetID(); + if(maxidtemp < 0 && fAODfilterbit != 128) continue; + if(maxidtemp > -1 && fAODfilterbit == 128) continue; + if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1; + if (maxidtemp > maxid1) maxid1 = maxidtemp; + if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){ + acctracks->Add(tr); + } + } + + maxid = maxid1; + + return acctracks; + +} + +//_________________________________________________________________________ +void AliEPSelectionTask::SetOADBandPeriod() +{ + TString oadbfilename; + + if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h + {fPeriod = "LHC10h"; + if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user + + if (fAnalysisInput.CompareTo("AOD")==0){ + oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath())); + } else if (fAnalysisInput.CompareTo("ESD")==0){ + oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath())); + } + + TFile foadb(oadbfilename); + if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data())); + + AliInfo("Using Standard OADB"); + fEPContainer = (AliOADBContainer*) foadb.Get("epphidist"); + if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection"); + foadb.Close(); + } + } + + if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h + {fPeriod = "LHC11h"; + if (!fUserphidist) { + // if it's already set and custom class is required, we use the one provided by the user + + oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath())); + TFile *foadb = TFile::Open(oadbfilename); + if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data())); + + AliInfo("Using Standard OADB"); + fSparseDist = (THnSparse*) foadb->Get("Default"); + if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection"); + foadb->Close(); + if(!fHruns){ + fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis; + fHruns->SetName("runsHisto"); + } + } + } +} + +//_________________________________________________________________________ +TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track) +{ + if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0]; + else if(fPeriod.CompareTo("LHC11h")==0) + { + if (track->Charge() < 0) + { + if(track->Eta() < 0.) return fPhiDist[0]; + else if (track->Eta() > 0.) return fPhiDist[2]; + } + else if (track->Charge() > 0) + { + if(track->Eta() < 0.) return fPhiDist[1]; + else if (track->Eta() > 0.) return fPhiDist[3]; + } + + } + return 0; +} + +TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd) +{ + // 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 + TObjArray *acctracks = new TObjArray(); + acctracks->SetOwner(kTRUE); + + const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD(); + + Float_t ptlow = 0; + Float_t ptup = 0; + Float_t etalow = 0; + Float_t etaup = 0; + fESDtrackCuts->GetPtRange(ptlow,ptup); + fESDtrackCuts->GetEtaRange(etalow,etaup); + + Double_t pDCA[3] = { 0. }; // momentum at DCA + Double_t rDCA[3] = { 0. }; // position at DCA + Float_t dDCA[2] = {0.}; // DCA to the vertex d and z + Float_t cDCA[3] = {0.}; // covariance of impact parameters + + + + for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){ + + AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy + // + // Track selection + if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue; + + // create a tpc only tracl + AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast(esd),esdTrack->GetID()); + if(!track) continue; + + if(track->Pt()>0.) + { + // only constrain tracks above threshold + AliExternalTrackParam exParam; + // take the B-field from the ESD, no 3D fieldMap available at this point + Bool_t relate = false; + relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam); + if(!relate){ + delete track; + continue; + } + // fetch the track parameters at the DCA (unconstraint) + if(track->GetTPCInnerParam()){ + track->GetTPCInnerParam()->GetPxPyPz(pDCA); + track->GetTPCInnerParam()->GetXYZ(rDCA); + } + // get the DCA to the vertex: + track->GetImpactParametersTPC(dDCA,cDCA); + // set the constrained parameters to the track + track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance()); + } + + + Float_t eta = track->Eta(); + Float_t pT = track->Pt(); + + if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){ + delete track; + continue; + } + + acctracks->Add(track); + } + + + return acctracks; + +} +