X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=ANALYSIS%2FAliEPSelectionTask.cxx;h=7d0d235c61d48c8aa98adea8557926064ab37dbd;hb=a1ce031d91be9bcb03e7de62e385e8ffbc01797f;hp=722ded872d8c729634099917c45ea08df0b9d932;hpb=5b53b816da5ffcbf851b82d59507c5df936f62d1;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliEPSelectionTask.cxx b/ANALYSIS/AliEPSelectionTask.cxx index 722ded872d8..7d0d235c61d 100644 --- a/ANALYSIS/AliEPSelectionTask.cxx +++ b/ANALYSIS/AliEPSelectionTask.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -61,6 +62,8 @@ #include "AliVTrack.h" #include "AliEventplane.h" +using std::cout; +using std::endl; ClassImp(AliEPSelectionTask) //________________________________________________________________________ @@ -68,6 +71,7 @@ AliEPSelectionTask::AliEPSelectionTask(): AliAnalysisTaskSE(), fAnalysisInput("ESD"), fTrackType("TPC"), + fPeriod(""), fUseMCRP(kFALSE), fUsePhiWeight(kFALSE), fUsePtWeight(kFALSE), @@ -80,7 +84,8 @@ AliAnalysisTaskSE(), fSplitMethod(0), fESDtrackCuts(0), fEPContainer(0), - fPhiDist(0), + fSparseDist(0), + fHruns(0), fQVector(0), fQContributionX(0), fQContributionY(0), @@ -100,6 +105,9 @@ AliAnalysisTaskSE(), { // Default constructor AliInfo("Event Plane Selection enabled."); + for(Int_t i = 0; i < 4; ++i) { + fPhiDist[i] = 0; + } } //________________________________________________________________________ @@ -107,6 +115,7 @@ AliEPSelectionTask::AliEPSelectionTask(const char *name): AliAnalysisTaskSE(name), fAnalysisInput("ESD"), fTrackType("TPC"), + fPeriod(""), fUseMCRP(kFALSE), fUsePhiWeight(kFALSE), fUsePtWeight(kFALSE), @@ -119,7 +128,8 @@ AliEPSelectionTask::AliEPSelectionTask(const char *name): fSplitMethod(0), fESDtrackCuts(0), fEPContainer(0), - fPhiDist(0), + fSparseDist(0), + fHruns(0), fQVector(0), fQContributionX(0), fQContributionY(0), @@ -140,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; + } } //________________________________________________________________________ @@ -155,15 +168,24 @@ AliEPSelectionTask::~AliEPSelectionTask() 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; + } } //________________________________________________________________________ @@ -194,26 +216,7 @@ 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; - - 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(); - } } @@ -228,7 +231,7 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) 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){ @@ -237,7 +240,9 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) if (esd){ if (!(fRunNumber == esd->GetRunNumber())) { fRunNumber = esd->GetRunNumber(); - SetPhiDist(); + AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber)); + SetOADBandPeriod(); + SetPhiDist(); } @@ -261,7 +266,8 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) 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){ @@ -291,8 +297,8 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) 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); @@ -302,6 +308,7 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) } fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ); } + tracklist->Clear(); delete tracklist; tracklist = 0; } @@ -311,51 +318,54 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/) AliVEvent* event = InputEvent(); AliAODEvent* aod = dynamic_cast(event); - if (!(fRunNumber == aod->GetRunNumber())) { - fRunNumber = aod->GetRunNumber(); - SetPhiDist(); - } - - if (fUseMCRP) { - AliAODMCHeader *headerH = dynamic_cast(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName())); - if (headerH) fRP = headerH->GetReactionPlaneAngle(); - } - 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); - } + 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 (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; @@ -551,7 +588,7 @@ void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){ } //________________________________________________________________________ -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; @@ -566,6 +603,7 @@ void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalo fUsercuts = kTRUE; fESDtrackCuts = new AliESDtrackCuts(); fESDtrackCuts->SetPtRange(ptlow,ptup); + fESDtrackCuts->SetMinNClustersTPC(ntpc); fESDtrackCuts->SetEtaRange(etalow,etaup); fAODfilterbit = filterbit; } @@ -605,17 +643,20 @@ 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; } @@ -625,37 +666,78 @@ Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1) //__________________________________________________________________________ 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; 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); } } @@ -667,8 +749,8 @@ void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char 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(); } @@ -688,6 +770,7 @@ TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& max 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); @@ -696,7 +779,7 @@ TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& max 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); } } @@ -706,3 +789,143 @@ TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& max 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; + +} +