#include <TList.h>
#include <TH1F.h>
#include <TH2F.h>
+#include <THnSparse.h>
#include <TProfile.h>
#include <TFile.h>
#include <TObjString.h>
#include "AliVTrack.h"
#include "AliEventplane.h"
+using std::cout;
+using std::endl;
ClassImp(AliEPSelectionTask)
//________________________________________________________________________
AliAnalysisTaskSE(),
fAnalysisInput("ESD"),
fTrackType("TPC"),
+ fPeriod(""),
fUseMCRP(kFALSE),
fUsePhiWeight(kFALSE),
fUsePtWeight(kFALSE),
fSplitMethod(0),
fESDtrackCuts(0),
fEPContainer(0),
- fPhiDist(0),
+ fSparseDist(0),
+ fHruns(0),
fQVector(0),
fQContributionX(0),
fQContributionY(0),
{
// Default constructor
AliInfo("Event Plane Selection enabled.");
+ for(Int_t i = 0; i < 4; ++i) {
+ fPhiDist[i] = 0;
+ }
}
//________________________________________________________________________
AliAnalysisTaskSE(name),
fAnalysisInput("ESD"),
fTrackType("TPC"),
+ fPeriod(""),
fUseMCRP(kFALSE),
fUsePhiWeight(kFALSE),
fUsePtWeight(kFALSE),
fSplitMethod(0),
fESDtrackCuts(0),
fEPContainer(0),
- fPhiDist(0),
+ fSparseDist(0),
+ fHruns(0),
fQVector(0),
fQContributionX(0),
fQContributionY(0),
// Default constructor
AliInfo("Event Plane Selection enabled.");
DefineOutput(1, TList::Class());
+ for(Int_t i = 0; i < 4; i++) {
+ fPhiDist[i] = 0;
+ }
}
//________________________________________________________________________
fESDtrackCuts = 0;
}
if (fUserphidist) {
- if (fPhiDist) {
- delete fPhiDist;
- fPhiDist = 0;
+ 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;
+ }
}
//________________________________________________________________________
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;
-
- 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();
- }
}
AliEventplane *esdEP;
TVector2 qq1;
TVector2 qq2;
- Double_t fRP = 0.; // the monte carlo reaction plane angle
+ Double_t fRP = 0.; // monte carlo reaction plane angle
if (fAnalysisInput.CompareTo("ESD")==0){
if (esd){
if (!(fRunNumber == esd->GetRunNumber())) {
fRunNumber = esd->GetRunNumber();
- SetPhiDist();
+ AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
+ SetOADBandPeriod();
+ SetPhiDist();
}
TObjArray* tracklist = new TObjArray;
if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
- if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
+ 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){
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);
}
fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
}
+ tracklist->Clear();
delete tracklist;
tracklist = 0;
}
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<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
- if (headerH) fRP = headerH->GetReactionPlaneAngle();
- }
+ if (fUseMCRP) {
+
AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+ if (headerH) fRP = headerH->GetReactionPlaneAngle();
+ }
- if (aod){
esdEP = aod->GetHeader()->GetEventplaneP();
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);
- }
+ 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();
+ 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.);
+ 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);
+ 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);
- 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<NT;iter++){
}
}
}
+ } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
+
+ for (Int_t i = 0; i < nt; i++) {
+ weight = 1;
+ track = dynamic_cast<AliVTrack*> (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;
}
//________________________________________________________________________
-void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
+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;
fUsercuts = kTRUE;
fESDtrackCuts = new AliESDtrackCuts();
fESDtrackCuts->SetPtRange(ptlow,ptup);
+ fESDtrackCuts->SetMinNClustersTPC(ntpc);
fESDtrackCuts->SetEtaRange(etalow,etaup);
fAODfilterbit = filterbit;
}
{
Double_t phiweight=1;
AliVTrack* track = dynamic_cast<AliVTrack*>(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;
}
//__________________________________________________________________________
void AliEPSelectionTask::SetPhiDist()
{
- 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));
+ 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
+ }
+
+ 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; i<fPhiDist->GetNbinsX(); 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; i<fPhiDist[0]->GetNbinsX(); 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);
}
}
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();
}
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);
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){
+ if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
acctracks->Add(tr);
}
}
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<AliESDEvent*>(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;
+
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff