#include "AliPhysicsSelection.h"
#include "AliBackgroundSelection.h"
#include "AliESDUtils.h"
+#include "AliOADBContainer.h"
#include "AliEventplane.h"
//________________________________________________________________________
AliEPSelectionTask::AliEPSelectionTask():
AliAnalysisTaskSE(),
- fDebug(0),
fAnalysisInput("ESD"),
- fStatus("TPC"),
+ fTrackType("TPC"),
fUseMCRP(kFALSE),
fUsePhiWeight(kFALSE),
fUsePtWeight(kFALSE),
fSaveTrackContribution(kFALSE),
+ fUserphidist(kFALSE),
+ fUsercuts(kFALSE),
+ fRunNumber(-15),
fESDtrackCuts(0),
- ftracklist(0),
+ fEPContainer(0),
fPhiDist(0),
fQVector(0),
fQContributionX(0),
//________________________________________________________________________
AliEPSelectionTask::AliEPSelectionTask(const char *name):
AliAnalysisTaskSE(name),
- fDebug(0),
fAnalysisInput("ESD"),
- fStatus("TPC"),
+ fTrackType("TPC"),
fUseMCRP(kFALSE),
fUsePhiWeight(kFALSE),
fUsePtWeight(kFALSE),
fSaveTrackContribution(kFALSE),
+ fUserphidist(kFALSE),
+ fUsercuts(kFALSE),
+ fRunNumber(-15),
fESDtrackCuts(0),
- ftracklist(0),
+ fEPContainer(0),
fPhiDist(0),
fQVector(0),
fQContributionX(0),
delete fESDtrackCuts;
fESDtrackCuts = 0;
}
- if (fQVector){
- delete fQVector;
- fQVector = 0;
- }
- if (fQsub1){
- delete fQsub1;
- fQsub1 = 0;
+ if (fUserphidist) {
+ if (fPhiDist) {
+ delete fPhiDist;
+ fPhiDist = 0;
+ }
}
- if (fQsub2){
- delete fQsub2;
- fQsub2 = 0;
+ if (fEPContainer){
+ delete fEPContainer;
+ fEPContainer = 0;
}
}
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();
+ }
}
//________________________________________________________________________
{
// Execute analysis for current event:
if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
+
+// fRunNumber = -15;
AliEventplane* esdEP = 0;
- TVector2 QQ1;
- TVector2 QQ2;
+ TVector2 qq1;
+ TVector2 qq2;
Double_t fRP = 0.; // the monte carlo reaction plane angle
if (fAnalysisInput.CompareTo("ESD")==0){
-
+
AliVEvent* event = InputEvent();
AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
-
- if (fUseMCRP) {
- AliMCEvent* mcEvent = MCEvent();
- if (mcEvent && mcEvent->GenEventHeader()) {
- AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
- fRP = headerH->ReactionPlaneAngle();
+ if (esd){
+ if (!(fRunNumber == esd->GetRunNumber())) {
+ fRunNumber = esd->GetRunNumber();
+ SetPhiDist();
}
- }
- if (esd){
+
+
+ if (fUseMCRP) {
+ AliMCEvent* mcEvent = MCEvent();
+ if (mcEvent && mcEvent->GenEventHeader()) {
+ AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
+ if (headerH) fRP = headerH->ReactionPlaneAngle();
+ }
+ }
+
esdEP = esd->GetEventplane();
if (fSaveTrackContribution) {
esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
}
- if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
- if (fStatus.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) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
+ 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);
+ 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);
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<NT;iter++){
- AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
- 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));
+ for (int iter = 0; iter<nt;iter++){
+ AliESDtrack* track = dynamic_cast<AliESDtrack*> (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));
+ }
}
AliESDtrack* trmax = esd->GetTrack(0);
- for (int iter = 1; iter<NT;iter++){
- AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
- if (track->Pt() > trmax->Pt()) trmax = track;
+ for (int iter = 1; iter<nt;iter++){
+ AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
+ if (track && (track->Pt() > trmax->Pt())) trmax = track;
}
fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
}
- delete ftracklist;
+ delete tracklist;
+ tracklist = 0;
}
}
//__________________________________________________________________________
TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
{
+// Get the Q vector
TVector2 mQ;
float mQx=0, mQy=0;
AliESDtrack* track;
Double_t weight;
- int NT = tracklist->GetEntries();
+ int nt = tracklist->GetEntries();
- for (int i=0; i<NT; i++){
+ for (int i=0; i<nt; i++){
weight = 1;
track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
- 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());
+ 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()));
}
- mQx += (weight*cos(2*track->Phi()));
- mQy += (weight*sin(2*track->Phi()));
}
mQ.Set(mQx,mQy);
return mQ;
//________________________________________________________________________
void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
{
+// Get Qsub
TVector2 mQ[2];
float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
Double_t weight;
AliESDtrack* track;
- TRandom2 RN = 0;
+ TRandom2 rn = 0;
- int NT = tracklist->GetEntries();
+ int nt = tracklist->GetEntries();
int trackcounter1=0, trackcounter2=0;
- for (Int_t i = 0; i < NT; i++) {
+ for (Int_t i = 0; i < nt; i++) {
weight = 1;
track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
+ if (!track) continue;
weight = GetWeight(track);
// This loop splits the track set into 2 random subsets
- if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
- float random = RN.Rndm();
+ 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()));
trackcounter2++;
}
}
- else if( trackcounter1 >= int(NT/2.)){
+ else if( trackcounter1 >= int(nt/2.)){
mQx2 += (weight*cos(2*track->Phi()));
mQy2 += (weight*sin(2*track->Phi()));
trackcounter2++;
}
//________________________________________________________________________
-void AliEPSelectionTask::SetESDtrackCuts(TString status){
+void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
- fStatus = status;
- if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
- if (fStatus.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+ if(fESDtrackCuts){
+ delete fESDtrackCuts;
+ fESDtrackCuts = 0;
+ }
+
+ fUsercuts = kTRUE;
+ fESDtrackCuts = trackcuts;
+}
+
+//_____________________________________________________________________________
+void AliEPSelectionTask::SetTrackType(TString tracktype){
+// Set the track type
+ 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);
fESDtrackCuts->SetEtaRange(-0.8,0.8);
+ }
}
-//__________________________________________________________________________
-void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
+//________________________________________________________________________
+Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
{
- TFile f(infilename);
- TObject* list = f.Get(listname);
- fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
- if (!fPhiDist) cout << "Phi Distribution not found!!!" << endl;
+// Get weight for track
+ Double_t ptweight=1;
+
+ if (fUsePtWeight) {
+ if (track->Pt()<2) ptweight=track->Pt();
+ else ptweight=2;
+ }
+ return ptweight*GetPhiWeight(track);
+}
+
+//________________________________________________________________________
+Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
+{
+// Get phi weight for track
+ Double_t phiweight=1;
+
+ if (fUsePhiWeight && fPhiDist && track) {
+ Double_t nParticles = fPhiDist->Integral();
+ Double_t nPhibins = fPhiDist->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()));
+
+ if (phiDistValue > 0) phiweight = nParticles/nPhibins/phiDistValue;
+ }
+ return phiweight;
+}
+
+//__________________________________________________________________________
+void AliEPSelectionTask::SetPhiDist()
+{
+// Set the phi distribution
+ if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
+
+ fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
+ if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
+
+ }
+ else {
+ AliInfo("Using Custom Phi Distribution");
+ }
+
Bool_t emptybins;
int iter = 0;
if (emptybins) {
AliError("After Maximum of rebinning still empty Phi-bins!!!");
}
- f.Close();
}
-//________________________________________________________________________
-Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
-{
- Double_t ptweight=1;
-
- if (fUsePtWeight) {
- if (track->Pt()<2) ptweight=track->Pt();
- else ptweight=2;
- }
- return ptweight*GetPhiWeight(track);
-}
-
-//________________________________________________________________________
-Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
+//__________________________________________________________________________
+void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
{
- Double_t phiweight=1;
-
- if (fUsePhiWeight) {
- Double_t nParticles = fPhiDist->Integral();
- Double_t nPhibins = fPhiDist->GetNbinsX();
+ // Set a personal phi distribution
+ fUserphidist = kTRUE;
- 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()));
-
- if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
- }
- return phiweight;
-}
+ TFile f(infilename);
+ TObject* list = f.Get(listname);
+ fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
+ if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
+
+ f.Close();
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff