//////////////////////////////////////////////////////////////////////////////// /// /// /// AliFemtoEventReaderESDChain - the reader class for the Alice ESD /// /// tailored for the Task framework /// /// Reads in AliESDfriend to create shared hit/quality information /// /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu /// /// /// //////////////////////////////////////////////////////////////////////////////// #include "AliFemtoEventReaderESDChain.h" #include "TFile.h" #include "TTree.h" #include "AliESDEvent.h" #include "AliESDtrack.h" #include "AliESDVertex.h" #include "AliMultiplicity.h" #include "AliCentrality.h" #include "AliEventplane.h" #include "AliESDVZERO.h" #include "AliFmPhysicalHelixD.h" #include "AliFmThreeVectorF.h" #include "SystemOfUnits.h" #include "AliFemtoEvent.h" #include "AliFemtoModelHiddenInfo.h" #include "AliPID.h" #include "AliAnalysisUtils.h" ClassImp(AliFemtoEventReaderESDChain) #if !(ST_NO_NAMESPACES) using namespace units; #endif using namespace std; //____________________________ AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(): fFileName(" "), fConstrained(true), fReadInner(false), fUseTPCOnly(false), fNumberofEvent(0), fCurEvent(0), fCurFile(0), fEvent(0x0), fTrackType(kGlobal), fEstEventMult(kReferenceITSTPC), fEventTrig(AliVEvent::kMB), //trigger fESDpid(0), fIsPidOwner(0), fReadV0(0), fMagFieldSign(0), fpA2013(kFALSE), fisPileUp(kFALSE), fMVPlp(kFALSE), fMinVtxContr(0), fMinPlpContribMV(0), fMinPlpContribSPD(0) { //constructor with 0 parameters , look at default settings // fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // } // fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // } } //__________________ AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader): AliFemtoEventReader(aReader), fFileName(" "), fConstrained(true), fReadInner(false), fUseTPCOnly(false), fNumberofEvent(0), fCurEvent(0), fCurFile(0), fEvent(0x0), fTrackType(kGlobal), fEstEventMult(kReferenceITSTPC), fEventTrig(AliVEvent::kMB), //trigger fESDpid(0), fIsPidOwner(0), fReadV0(0), fMagFieldSign(0), fpA2013(kFALSE), fisPileUp(kFALSE), fMVPlp(kFALSE), fMinVtxContr(0), fMinPlpContribMV(0), fMinPlpContribSPD(0) { // Copy constructor fConstrained = aReader.fConstrained; fReadInner = aReader.fReadInner; fUseTPCOnly = aReader.fUseTPCOnly; fNumberofEvent = aReader.fNumberofEvent; fCurEvent = aReader.fCurEvent; fCurFile = aReader.fCurFile; // fEvent = new AliESD(*aReader.fEvent); fEvent = new AliESDEvent(); fTrackType = aReader.fTrackType; fEstEventMult = aReader.fEstEventMult; fEventTrig = aReader.fEventTrig; //trigger fReadV0 = aReader.fReadV0; fMagFieldSign = aReader.fMagFieldSign; fpA2013 = aReader.fpA2013; fisPileUp = aReader.fisPileUp; fMVPlp = aReader.fMVPlp; fMinVtxContr = aReader.fMinVtxContr; fMinPlpContribMV = aReader.fMinPlpContribMV; fMinPlpContribSPD = aReader.fMinPlpContribSPD; // fEventFriend = aReader.fEventFriend; // fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // list::iterator iter; // for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) { // fClusterPerPadrow[tPad]->push_back(*iter); // } // } // fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // list::iterator iter; // for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) { // fSharedList[tPad]->push_back(*iter); // } // } } //__________________ AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain() { //Destructor delete fEvent; // for (int tPad=0; tPadclear(); // delete fClusterPerPadrow[tPad]; // } // delete [] fClusterPerPadrow; // for (int tPad=0; tPadclear(); // delete fSharedList[tPad]; // } // delete [] fSharedList; } //__________________ AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader) { // Assignment operator if (this == &aReader) return *this; fConstrained = aReader.fConstrained; fReadInner = aReader.fReadInner; fUseTPCOnly = aReader.fUseTPCOnly; fNumberofEvent = aReader.fNumberofEvent; fCurEvent = aReader.fCurEvent; fCurFile = aReader.fCurFile; if (fEvent) delete fEvent; fEvent = new AliESDEvent(); fTrackType = aReader.fTrackType; fEstEventMult = aReader.fEstEventMult; fReadV0 = aReader.fReadV0; fMagFieldSign = aReader.fMagFieldSign; fpA2013 = aReader.fpA2013; fisPileUp = aReader.fisPileUp; fMVPlp = aReader.fMVPlp; // fEventFriend = aReader.fEventFriend; // if (fClusterPerPadrow) { // for (int tPad=0; tPadclear(); // delete fClusterPerPadrow[tPad]; // } // delete [] fClusterPerPadrow; // } // if (fSharedList) { // for (int tPad=0; tPadclear(); // delete fSharedList[tPad]; // } // delete [] fSharedList; // } // fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // list::iterator iter; // for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) { // fClusterPerPadrow[tPad]->push_back(*iter); // } // } // fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); // for (int tPad=0; tPad(); // list::iterator iter; // for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) { // fSharedList[tPad]->push_back(*iter); // } // } return *this; } //__________________ // Simple report AliFemtoString AliFemtoEventReaderESDChain::Report() { AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n"; return temp; } //__________________ void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained) { // Select whether to read constrained or not constrained momentum fConstrained=constrained; } //__________________ bool AliFemtoEventReaderESDChain::GetConstrained() const { // Check whether we read constrained or not constrained momentum return fConstrained; } //__________________ void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner) { fReadInner=readinner; } bool AliFemtoEventReaderESDChain::GetReadTPCInner() const { return fReadInner; } //__________________ void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly) { fUseTPCOnly=usetpconly; } bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const { return fUseTPCOnly; } void AliFemtoEventReaderESDChain::SetUseMultiplicity(EstEventMult aType) { fEstEventMult = aType; } AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent() { // Get the event, read all the relevant information // and fill the AliFemtoEvent class // Returns a valid AliFemtoEvent AliFemtoEvent *hbtEvent = 0; hbtEvent = new AliFemtoEvent; CopyESDtoFemtoEvent(hbtEvent); fCurEvent++; return hbtEvent; } void AliFemtoEventReaderESDChain::CopyESDtoFemtoEvent(AliFemtoEvent *hbtEvent) { //string tFriendFileName; // Get the friend information if (Debug()>1) cout<<"starting to read event "<SetESDfriend(fEventFriend); if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD(); //setting basic things // hbtEvent->SetEventNumber(fEvent->GetEventNumber()); hbtEvent->SetRunNumber(fEvent->GetRunNumber()); //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks()); hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy()); hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy()); hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy()); hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy()); hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy()); hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants()); hbtEvent->SetTriggerMask(fEvent->GetTriggerMask()); // hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster()); if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) || (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) || (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) || (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD"))) hbtEvent->SetTriggerCluster(1); else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) || (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD"))) hbtEvent->SetTriggerCluster(2); else hbtEvent->SetTriggerCluster(0); //Vertex double fV1[3]; double fVCov[6]; //AliAnalysisUtils if(fisPileUp||fpA2013) { AliAnalysisUtils *anaUtil=new AliAnalysisUtils(); if(fMinVtxContr) anaUtil->SetMinVtxContr(fMinVtxContr); if(fpA2013) if(anaUtil->IsVertexSelected2013pA(fEvent)==kFALSE) return; //Vertex rejection for pA analysis. if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE); else anaUtil->SetUseMVPlpSelection(kFALSE); if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV); if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD); if(fisPileUp) if(anaUtil->IsPileUpEvent(fEvent)) return; //Pile-up rejection. delete anaUtil; } if (fUseTPCOnly) { const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertexTPC(); if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event. esdvertex->GetXYZ(fV1); esdvertex->GetCovMatrix(fVCov); if (!esdvertex->GetStatus()) fVCov[4] = -1001.0; } else { const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertex(); if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event. esdvertex->GetXYZ(fV1); esdvertex->GetCovMatrix(fVCov); if (!esdvertex->GetStatus()) fVCov[4] = -1001.0; } AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]); hbtEvent->SetPrimVertPos(vertex); hbtEvent->SetPrimVertCov(fVCov); Int_t spdetaonecount = 0; for (int iter=0; iterGetMultiplicity()->GetNumberOfTracklets(); iter++) if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0) spdetaonecount++; // hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets()); hbtEvent->SetSPDMult(spdetaonecount); //starting to reading tracks int nofTracks=0; //number of reconstructed tracks in event nofTracks=fEvent->GetNumberOfTracks(); int realnofTracks=0;//number of track which we use ina analysis // // Clear the shared cluster list // for (int tPad=0; tPadclear(); // } // for (int tPad=0; tPadclear(); // } // for (int i=0;iGetTrack(i);//getting next track // list::iterator tClustIter; // Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster]; // Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices); // for (int tNcl=0; tNcl= 0) { // tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]); // if (tClustIter == fClusterPerPadrow[tNcl]->end()) { // fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]); // } // else { // fSharedList[tNcl]->push_back(tTrackIndices[tNcl]); // } // } // } // } int tNormMult = 0; int tNormMultPos = 0; int tNormMultNeg = 0; Float_t tTotalPt = 0.0; Float_t b[2]; Float_t bCov[3]; Int_t tTracklet=0, tITSTPC=0; //W-AliESDEvent::EstimateMultiplicity: This obsolete method will be eliminated soon. Use AliESDtrackCuts::GetReferenceMultiplicity //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2); hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC); hbtEvent->SetMultiplicityEstimateTracklets(tTracklet); // hbtEvent->SetMultiplicityEstimateITSPure(tITSPure); hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1)); for (int i=0;iGetTrack(i);//getting next track // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack(); if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) && (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) { if (esdtrack->GetTPCNcls() > 70) if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) { if (esdtrack->Pt() > 0.15 && esdtrack->Pt() < 20) if (TMath::Abs(esdtrack->Eta()) < 0.8) { esdtrack->GetImpactParameters(b,bCov); if ((b[0]<0.2) && (b[1] < 2.0)) { tNormMult++; tTotalPt += esdtrack->Pt(); } } } } hbtEvent->SetZDCEMEnergy(tTotalPt); // if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) // if (esdtrack->GetTPCNcls() > 80) // if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0) // if (esdtrack->GetConstrainedParam()) // if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5) // if (esdtrack->GetConstrainedParam()->Pt() < 1.0) { // if (esdtrack->GetSign() > 0) // tNormMultPos++; // else if (esdtrack->GetSign() < 0) // tNormMultNeg--; // } // If reading ITS-only tracks, reject all with TPC if (fTrackType == kITSOnly) { if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue; if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue; if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue; UChar_t iclm = esdtrack->GetITSClusterMap(); Int_t incls = 0; for (int iter=0; iter<6; iter++) if (iclm&(1<1) cout << "Rejecting track with " << incls << " clusters" << endl; continue; } } AliFemtoTrack* trackCopy = new AliFemtoTrack(); trackCopy->SetCharge((short)esdtrack->GetSign()); //in aliroot we have AliPID //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon //we use only 5 first double esdpid[5]; // esdtrack->GetESDpid(esdpid); esdtrack->GetTPCpid(esdpid); trackCopy->SetPidProbElectron(esdpid[0]); trackCopy->SetPidProbMuon(esdpid[1]); trackCopy->SetPidProbPion(esdpid[2]); trackCopy->SetPidProbKaon(esdpid[3]); trackCopy->SetPidProbProton(esdpid[4]); esdpid[0] = -100000.0; esdpid[1] = -100000.0; esdpid[2] = -100000.0; esdpid[3] = -100000.0; esdpid[4] = -100000.0; double tTOF = 0.0; if (esdtrack->GetStatus()&AliESDtrack::kTOFout&AliESDtrack::kTIME) { tTOF = esdtrack->GetTOFsignal(); esdtrack->GetIntegratedTimes(esdpid); } trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]); ////// TPC //////////////////////////////////////////// float nsigmaTPCK=-1000.; float nsigmaTPCPi=-1000.; float nsigmaTPCP=-1000.; if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){ nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon); nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion); nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton); } trackCopy->SetNSigmaTPCPi(nsigmaTPCPi); trackCopy->SetNSigmaTPCK(nsigmaTPCK); trackCopy->SetNSigmaTPCP(nsigmaTPCP); ///// TOF /////////////////////////////////////////////// float vp=-1000.; float nsigmaTOFPi=-1000.; float nsigmaTOFK=-1000.; float nsigmaTOFP=-1000.; if (// (esdtrack->GetStatus()&AliESDtrack::kTOFpid) && (esdtrack->GetStatus()&AliESDtrack::kTOFout) && (esdtrack->GetStatus()&AliESDtrack::kTIME)) { //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) && //(esdtrack->GetStatus()&AliESDtrack::kTOFout) && //(esdtrack->GetStatus()&AliESDtrack::kTIME)){ // collect info from ESDpid class if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME))) { double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P()); nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero); nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero); nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero); Double_t len=esdtrack->GetIntegratedLength(); Double_t tof=esdtrack->GetTOFsignal(); if(tof > 0.) vp=len/tof/0.03; } } trackCopy->SetVTOF(vp); trackCopy->SetNSigmaTOFPi(nsigmaTOFPi); trackCopy->SetNSigmaTOFK(nsigmaTOFK); trackCopy->SetNSigmaTOFP(nsigmaTOFP); double pxyz[3]; double rxyz[3]; double impact[2]; double covimpact[3]; if (fUseTPCOnly) { if (!esdtrack->GetTPCInnerParam()) { delete trackCopy; continue; } AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam()); param->GetXYZ(rxyz); param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact); param->GetPxPyPz(pxyz);//reading noconstarined momentum if (fReadInner == true) { AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo(); tInfo->SetPDGPid(211); tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]); tInfo->SetMass(0.13957); // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0); // tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0); tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0); trackCopy->SetHiddenInfo(tInfo); } AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); if (v.Mag() < 0.0001) { // cout << "Found 0 momentum ???? " <SetP(v);//setting momentum trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]); const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]); //setting helix I do not if it is ok AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); trackCopy->SetHelix(helix); //some stuff which could be useful trackCopy->SetImpactD(impact[0]); trackCopy->SetImpactZ(impact[1]); trackCopy->SetCdd(covimpact[0]); trackCopy->SetCdz(covimpact[1]); trackCopy->SetCzz(covimpact[2]); trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact)); delete param; } else { if (fReadInner == true) { if (esdtrack->GetTPCInnerParam()) { AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetInnerParam()); //trackCopy->SetInnerMomentum(param->P()); trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum()); param->GetXYZ(rxyz); // param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000); param->GetPxPyPz(pxyz);//reading noconstarined momentum delete param; AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo(); tInfo->SetPDGPid(211); tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]); tInfo->SetMass(0.13957); // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0); //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0); tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0); trackCopy->SetHiddenInfo(tInfo); } } if (fTrackType == kGlobal) { if (fConstrained==true) tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum else tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum } else if (fTrackType == kTPCOnly) { if (esdtrack->GetTPCInnerParam()) esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz); else { delete trackCopy; continue; } } else if (fTrackType == kITSOnly) { if (fConstrained==true) tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum else tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum } AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); if (v.Mag() < 0.0001) { // cout << "Found 0 momentum ???? " <SetP(v);//setting momentum trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]); const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]); //setting helix I do not if it is ok AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); trackCopy->SetHelix(helix); //some stuff which could be useful float imp[2]; float cim[3]; // if (fTrackType == kTPCOnly) { // esdtrack->GetTPCInnerParam()->GetImpactParameters(imp,cim); // } // else { esdtrack->GetImpactParameters(imp,cim); // } impact[0] = imp[0]; impact[1] = imp[1]; covimpact[0] = cim[0]; covimpact[1] = cim[1]; covimpact[2] = cim[2]; trackCopy->SetImpactD(impact[0]); trackCopy->SetImpactZ(impact[1]); trackCopy->SetCdd(covimpact[0]); trackCopy->SetCdz(covimpact[1]); trackCopy->SetCzz(covimpact[2]); trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact)); } trackCopy->SetTrackId(esdtrack->GetID()); trackCopy->SetFlags(esdtrack->GetStatus()); trackCopy->SetLabel(esdtrack->GetLabel()); trackCopy->SetITSchi2(esdtrack->GetITSchi2()); if (esdtrack->GetITSFakeFlag()) trackCopy->SetITSncls(-esdtrack->GetNcls(0)); else trackCopy->SetITSncls(esdtrack->GetNcls(0)); trackCopy->SetTPCchi2(esdtrack->GetTPCchi2()); trackCopy->SetTPCncls(esdtrack->GetTPCNcls()); trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF()); trackCopy->SetTPCsignal(esdtrack->GetTPCsignal()); trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap()); trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap()); double xtpc[3]; esdtrack->GetInnerXYZ(xtpc); xtpc[2] -= fV1[2]; trackCopy->SetNominalTPCEntrancePoint(xtpc); esdtrack->GetOuterXYZ(xtpc); xtpc[2] -= fV1[2]; trackCopy->SetNominalTPCExitPoint(xtpc); int indexes[3]; for (int ik=0; ik<3; ik++) { indexes[ik] = esdtrack->GetKinkIndex(ik); } trackCopy->SetKinkIndexes(indexes); for (int ii=0; ii<6; ii++){ trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii)); } //decision if we want this track //if we using diffrent labels we want that this label was use for first time //if we use hidden info we want to have match between sim data and ESD if (tGoodMomentum==true) { hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis realnofTracks++;//real number of tracks // delete trackCopy; } else { delete trackCopy; } } hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event AliCentrality *cent = fEvent->GetCentrality(); if (cent) { hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M")); hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A")); hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C")); hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA")); hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC")); hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1")); hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0")); hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL")); hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD")); hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA")); // hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1")); hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK")); hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND")); if (Debug()>1) printf(" FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0); } if (fEstEventMult == kGlobalCount) hbtEvent->SetNormalizedMult(tNormMult); else if (fEstEventMult == kReferenceITSTPC) hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.0)); else if(fEstEventMult == kReferenceITSSA) hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.0)); else if(fEstEventMult == kReferenceTracklets) hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.0)); else if (fEstEventMult == kSPDLayer1) hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1)); else if (fEstEventMult == kVZERO) { Float_t multV0 = 0; for (Int_t i=0; i<64; i++) multV0 += fEvent->GetVZEROData()->GetMultiplicity(i); hbtEvent->SetNormalizedMult(multV0); } else if (fEstEventMult == kCentrality) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("V0M") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M"))); } } else if (fEstEventMult == kCentralityV0A) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("V0A") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A"))); } } else if (fEstEventMult == kCentralityV0C) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("V0C") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C"))); } } else if (fEstEventMult == kCentralityZNA) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("ZNA") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA"))); } } else if (fEstEventMult == kCentralityZNC) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("ZNC") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC"))); } } else if (fEstEventMult == kCentralityCL1) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("CL1") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1"))); } } else if (fEstEventMult == kCentralityCL0) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("CL0") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0"))); } } else if (fEstEventMult == kCentralityTRK) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("TRK") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK"))); } } else if (fEstEventMult == kCentralityTKL) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("TKL") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL"))); } } else if (fEstEventMult == kCentralityNPA) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("NPA") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA"))); } } else if (fEstEventMult == kCentralityCND) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("CND") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND"))); } } else if (fEstEventMult == kCentralityFMD) { // centrality between 0 (central) and 1 (very peripheral) if (cent) { if (cent->GetCentralityPercentile("FMD") < 0.00001) hbtEvent->SetNormalizedMult(-1); else hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD"))); if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD"))); } } if (tNormMultPos > tNormMultNeg) hbtEvent->SetZDCParticipants(tNormMultPos); else hbtEvent->SetZDCParticipants(tNormMultNeg); AliEventplane* ep = fEvent->GetEventplane(); if (ep) { hbtEvent->SetEP(ep); hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q")); } //V0 if(fReadV0) { for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) { AliESDv0* esdv0 = fEvent->GetV0(i); if (!esdv0) continue; //if(esdv0->GetNDaughters()>2) continue; //if(esdv0->GetNProngs()>2) continue; if(esdv0->Charge()!=0) continue; AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex()); if(!trackPos) continue; AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex()); if(!trackNeg) continue; if(trackPos->Charge()==trackNeg->Charge()) continue; AliFemtoV0* trackCopyV0 = new AliFemtoV0(); CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent); hbtEvent->V0Collection()->push_back(trackCopyV0); //cout<<"Pushback v0 to v0collection"<1) cout<<"end of reading nt "< Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma) // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma) // // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2) // Can this be expressed in a different way? if (bRes[0] == 0 || bRes[1] ==0) return -1; Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); // stupid rounding problem screws up everything: // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :( if (TMath::Exp(-d * d / 2) < 1e-10) return 1000; d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); return d; } void AliFemtoEventReaderESDChain::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent) { double fPrimaryVtxPosition[3]; tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition); tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2])); //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X()); //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y()); //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0 //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z()); //tFemtoV0->SetdecayVertexV0(decayvertex); tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters()); tFemtoV0->SetmomV0X(tESDv0->Px()); tFemtoV0->SetmomV0Y(tESDv0->Py()); tFemtoV0->SetmomV0Z(tESDv0->Pz()); AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz()); tFemtoV0->SetmomV0(momv0); tFemtoV0->SetalphaV0(tESDv0->AlphaV0()); tFemtoV0->SetptArmV0(tESDv0->PtArmV0()); //tFemtoV0->SeteLambda(tESDv0->ELambda()); //tFemtoV0->SeteK0Short(tESDv0->EK0Short()); //tFemtoV0->SetePosProton(tESDv0->EPosProton()); //tFemtoV0->SeteNegProton(tESDv0->ENegProton()); tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2)); tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4)); tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2)); //tFemtoV0->SetrapLambda(tESDv0->RapLambda()); //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short()); tFemtoV0->SetptV0(tESDv0->Pt()); tFemtoV0->SetptotV0(tESDv0->P()); tFemtoV0->SetEtaV0(tESDv0->Eta()); tFemtoV0->SetPhiV0(tESDv0->Phi()); tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2])); tFemtoV0->SetYV0(tESDv0->Y()); AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0); AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1); if(trackpos && trackneg) { tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField()))); tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField()))); double MomPos[3]; trackpos->PxPyPz(MomPos); tFemtoV0->SetmomPosX(MomPos[0]); tFemtoV0->SetmomPosY(MomPos[1]); tFemtoV0->SetmomPosZ(MomPos[2]); AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]); tFemtoV0->SetmomPos(mompos); double MomNeg[3]; trackneg->PxPyPz(MomNeg); tFemtoV0->SetmomNegX(MomNeg[0]); tFemtoV0->SetmomNegY(MomNeg[1]); tFemtoV0->SetmomNegZ(MomNeg[2]); AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]); tFemtoV0->SetmomNeg(momneg); tFemtoV0->SetptPos(trackpos->Pt()); tFemtoV0->SetptotPos(trackpos->P()); tFemtoV0->SetptNeg(trackneg->Pt()); tFemtoV0->SetptotNeg(trackneg->P()); tFemtoV0->SetidNeg(trackneg->GetID()); //cout<<"tESDv0->GetNegID(): "<GetNegID()<IdNeg(): "<IdNeg()<SetidPos(trackpos->GetID()); tFemtoV0->SetEtaPos(trackpos->Eta()); tFemtoV0->SetEtaNeg(trackneg->Eta()); tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos() tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg() tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls()); tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls()); tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap()); tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap()); tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap()); tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap()); tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls()); tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls()); tFemtoV0->SetStatusPos(trackpos->GetStatus()); tFemtoV0->SetStatusNeg(trackneg->GetStatus()); float bfield = 5*fMagFieldSign; float globalPositionsAtRadiiPos[9][3]; GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos); double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]}; double tpcExitPos[3]={globalPositionsAtRadiiPos[7][0],globalPositionsAtRadiiPos[7][1],globalPositionsAtRadiiPos[7][2]}; float globalPositionsAtRadiiNeg[9][3]; GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg); double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]}; double tpcExitNeg[3]={globalPositionsAtRadiiNeg[7][0],globalPositionsAtRadiiNeg[7][1],globalPositionsAtRadiiNeg[7][2]}; AliFemtoThreeVector tmpVec; tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]); tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec); tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]); tFemtoV0->SetNominalTpcExitPointPos(tmpVec); tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]); tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec); tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]); tFemtoV0->SetNominalTpcExitPointNeg(tmpVec); AliFemtoThreeVector vecTpcPos[9]; AliFemtoThreeVector vecTpcNeg[9]; for(int i=0;i<9;i++) { vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]); vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]); } tFemtoV0->SetNominalTpcPointPos(vecTpcPos); tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg); tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum(); tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum(); tFemtoV0->SetdedxPos(trackpos->GetTPCsignal()); tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal()); if (fESDpid) { tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon)); tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon)); tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton)); tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton)); tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion)); tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion)); } else { tFemtoV0->SetPosNSigmaTPCK(-1000); tFemtoV0->SetNegNSigmaTPCK(-1000); tFemtoV0->SetPosNSigmaTPCP(-1000); tFemtoV0->SetNegNSigmaTPCP(-1000); tFemtoV0->SetPosNSigmaTPCPi(-1000); tFemtoV0->SetNegNSigmaTPCPi(-1000); } if(// (tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0) { if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0) { tFemtoV0->SetPosNSigmaTOFK(-1000); tFemtoV0->SetNegNSigmaTOFK(-1000); tFemtoV0->SetPosNSigmaTOFP(-1000); tFemtoV0->SetNegNSigmaTOFP(-1000); tFemtoV0->SetPosNSigmaTOFPi(-1000); tFemtoV0->SetNegNSigmaTOFPi(-1000); } } else { if (fESDpid) { tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon)); tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon)); tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton)); tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton)); tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion)); tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion)); } else { tFemtoV0->SetPosNSigmaTOFK(-1000); tFemtoV0->SetNegNSigmaTOFK(-1000); tFemtoV0->SetPosNSigmaTOFP(-1000); tFemtoV0->SetNegNSigmaTOFP(-1000); tFemtoV0->SetPosNSigmaTOFPi(-1000); tFemtoV0->SetNegNSigmaTOFPi(-1000); } } } else { tFemtoV0->SetStatusPos(999); tFemtoV0->SetStatusNeg(999); } tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus()); } void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType) { fTrackType = aType; } //trigger void AliFemtoEventReaderESDChain::SetEventTrigger(UInt_t eventtrig) { fEventTrig = eventtrig; } //V0 reading void AliFemtoEventReaderESDChain::SetReadV0(bool a) { fReadV0 = a; } void AliFemtoEventReaderESDChain::SetMagneticFieldSign(int s) { if(s>0) fMagFieldSign = 1; else if(s<0) fMagFieldSign = -1; else fMagFieldSign = 0; } void AliFemtoEventReaderESDChain::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3]) { // Gets the global position of the track at nine different radii in the TPC // track is the track you want to propagate // bfield is the magnetic field of your event // globalPositionsAtRadii is the array of global positions in the radii and xyz // Initialize the array to something indicating there was no propagation for(Int_t i=0;i<9;i++){ for(Int_t j=0;j<3;j++){ globalPositionsAtRadii[i][j]=-9999.; } } // Make a copy of the track to not change parameters of the track AliExternalTrackParam etp; etp.CopyFromVTrack(track); //printf("\nAfter CopyFromVTrack\n"); //etp.Print(); // The global position of the the track Double_t xyz[3]={-9999.,-9999.,-9999.}; // Counter for which radius we want Int_t iR=0; // The radii at which we get the global positions // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm) Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.}; // The global radius we are at Float_t globalRadius=0; // Propagation is done in local x of the track for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent. // If the track's momentum is smaller than infinite, it will develop a y-component, which // adds to the global radius // Stop if the propagation was not succesful. This can happen for low pt tracks // that don't reach outer radii if(!etp.PropagateTo(x,bfield))break; etp.GetXYZ(xyz); // GetXYZ returns global coordinates globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii // Roughly reached the radius we want if(globalRadius > Rwanted[iR]){ // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps. while (globalRadius>Rwanted[iR]){ x-=.1; // printf("propagating to x %5.2f\n",x); if(!etp.PropagateTo(x,bfield))break; etp.GetXYZ(xyz); // GetXYZ returns global coordinates globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii } //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]); globalPositionsAtRadii[iR][0]=xyz[0]; globalPositionsAtRadii[iR][1]=xyz[1]; globalPositionsAtRadii[iR][2]=xyz[2]; // Indicate we want the next radius iR+=1; } if(iR>=8){ // TPC edge reached return; } } } void AliFemtoEventReaderESDChain::SetpA2013(Bool_t pA2013) { fpA2013 = pA2013; } void AliFemtoEventReaderESDChain::SetUseMVPlpSelection(Bool_t mvplp) { fMVPlp = mvplp; } void AliFemtoEventReaderESDChain::SetIsPileUpEvent(Bool_t ispileup) { fisPileUp = ispileup; }