From a68b227d8c642908055829169f1921cbf987c1c4 Mon Sep 17 00:00:00 2001 From: maszyman Date: Mon, 3 Mar 2014 17:57:46 +0100 Subject: [PATCH] fix in reading TPC signal in MC analysis (AliFemtoEventReaderESDChainKine class) --- .../AliFemtoEventReaderESDChainKine.cxx | 1386 +++++++++-------- 1 file changed, 694 insertions(+), 692 deletions(-) diff --git a/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx b/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx index e21360f602c..ef2a23d195e 100644 --- a/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx +++ b/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx @@ -39,7 +39,7 @@ ClassImp(AliFemtoEventReaderESDChainKine) #if !(ST_NO_NAMESPACES) - using namespace units; +using namespace units; #endif using namespace std; @@ -226,7 +226,7 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD"))) hbtEvent->SetTriggerCluster(1); else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) || - (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD"))) + (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD"))) hbtEvent->SetTriggerCluster(2); else hbtEvent->SetTriggerCluster(0); @@ -259,17 +259,17 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cdh) { TList *tGenHeaders = cdh->GetHeaders(); for (int ihead = 0; iheadGetEntries(); ihead++) { - hdh = dynamic_cast (fGenHeader); - if (hdh) break; + hdh = dynamic_cast (fGenHeader); + if (hdh) break; } } } if (hdh) - { - tReactionPlane = hdh->ReactionPlaneAngle(); - cout << "Got reaction plane " << tReactionPlane << endl; - } + { + tReactionPlane = hdh->ReactionPlaneAngle(); + cout << "Got reaction plane " << tReactionPlane << endl; + } hbtEvent->SetReactionPlaneAngle(tReactionPlane); @@ -315,9 +315,9 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() for (int ip=0; ipGetNtrack(); ip++) { motherpart = fStack->Particle(ip); if (motherpart->GetDaughter(0) > 0) - motherids[motherpart->GetDaughter(0)] = ip; + motherids[motherpart->GetDaughter(0)] = ip; if (motherpart->GetDaughter(1) > 0) - motherids[motherpart->GetDaughter(1)] = ip; + motherids[motherpart->GetDaughter(1)] = ip; // if (motherpart->GetPdgCode() == 211) { // cout << "Mother " << ip << " has daughters " @@ -338,125 +338,125 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } 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 (TMath::Abs(esdtrack->Eta()) < 1.2) { - esdtrack->GetImpactParameters(b,bCov); - if ((b[0]<0.2) && (b[1] < 0.25)) { - tNormMult++; - tTotalPt += esdtrack->Pt(); - } - } - } - } - else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) { - if (esdtrack->GetTPCNcls() > 100) - if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) { - if (TMath::Abs(esdtrack->Eta()) < 1.2) { - esdtrack->GetImpactParameters(b,bCov); - if ((b[0]<2.4) && (b[1] < 3.2)) { - tNormMult++; - tTotalPt += esdtrack->Pt(); - } - } - } - } + { + // cout << "Reading track " << i << endl; + bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track + + + const AliESDtrack *esdtrack=fEvent->GetTrack(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 (TMath::Abs(esdtrack->Eta()) < 1.2) { + esdtrack->GetImpactParameters(b,bCov); + if ((b[0]<0.2) && (b[1] < 0.25)) { + tNormMult++; + tTotalPt += esdtrack->Pt(); + } + } + } + } + else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) { + if (esdtrack->GetTPCNcls() > 100) + if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) { + if (TMath::Abs(esdtrack->Eta()) < 1.2) { + esdtrack->GetImpactParameters(b,bCov); + if ((b[0]<2.4) && (b[1] < 3.2)) { + 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; - } + 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()); + 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); - trackCopy->SetPidProbElectron(esdpid[0]); - trackCopy->SetPidProbMuon(esdpid[1]); - trackCopy->SetPidProbPion(esdpid[2]); - trackCopy->SetPidProbKaon(esdpid[3]); - trackCopy->SetPidProbProton(esdpid[4]); + //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); + 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; + esdpid[0] = -100000.0; + esdpid[1] = -100000.0; + esdpid[2] = -100000.0; + esdpid[3] = -100000.0; + esdpid[4] = -100000.0; - double tTOF = 0.0; + double tTOF = 0.0; - if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) { - tTOF = esdtrack->GetTOFsignal(); - esdtrack->GetIntegratedTimes(esdpid); - } + if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) { + tTOF = esdtrack->GetTOFsignal(); + esdtrack->GetIntegratedTimes(esdpid); + } - trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]); + trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]); - ////// TPC //////////////////////////////////////////// - float nsigmaTPCK=-1000.; - float nsigmaTPCPi=-1000.; - float nsigmaTPCP=-1000.; + ////// 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); + 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); + } + trackCopy->SetNSigmaTPCPi(nsigmaTPCPi); + trackCopy->SetNSigmaTPCK(nsigmaTPCK); + trackCopy->SetNSigmaTPCP(nsigmaTPCP); - ///// TOF /////////////////////////////////////////////// + ///// TOF /////////////////////////////////////////////// - float vp=-1000.; - float nsigmaTOFPi=-1000.; - float nsigmaTOFK=-1000.; - float nsigmaTOFP=-1000.; + 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)) { //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) && @@ -464,7 +464,7 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() //(esdtrack->GetStatus()&AliESDtrack::kTIME)){ // collect info from ESDpid class - if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) { + if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) { double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P()); @@ -479,360 +479,362 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } } - trackCopy->SetVTOF(vp); - trackCopy->SetNSigmaTOFPi(nsigmaTOFPi); - trackCopy->SetNSigmaTOFK(nsigmaTOFK); - trackCopy->SetNSigmaTOFP(nsigmaTOFP); + trackCopy->SetVTOF(vp); + trackCopy->SetNSigmaTOFPi(nsigmaTOFPi); + trackCopy->SetNSigmaTOFK(nsigmaTOFK); + trackCopy->SetNSigmaTOFP(nsigmaTOFP); - double pxyz[3]; - double rxyz[3]; - double impact[2]; - double covimpact[3]; + double pxyz[3]; + double rxyz[3]; + double impact[2]; + double covimpact[3]; - if (fUseTPCOnly) { - if (!esdtrack->GetTPCInnerParam()) { - cout << "No TPC inner param !" << endl; - 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); - } + if (fUseTPCOnly) { + if (!esdtrack->GetTPCInnerParam()) { + cout << "No TPC inner param !" << endl; + delete trackCopy; + continue; + } - if (fRotateToEventPlane) { - double tPhi = TMath::ATan2(pxyz[1], pxyz[0]); - double tRad = TMath::Hypot(pxyz[0], pxyz[1]); - - pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane); - pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane); - } - - AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); - if (v.Mag() < 0.0001) { - // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl; - delete trackCopy; - continue; - } - trackCopy->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; + 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); } - else { - if (fReadInner == true) { - - if (esdtrack->GetTPCInnerParam()) { - AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam()); - 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 (fRotateToEventPlane) { + double tPhi = TMath::ATan2(pxyz[1], pxyz[0]); + double tRad = TMath::Hypot(pxyz[0], pxyz[1]); - 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 - } + pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane); + pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane); + } - if (fRotateToEventPlane) { - double tPhi = TMath::ATan2(pxyz[1], pxyz[0]); - double tRad = TMath::Hypot(pxyz[0], pxyz[1]); - - pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane); - pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane); - } - - AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); - if (v.Mag() < 0.0001) { - // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl; - delete trackCopy; - continue; - } - - trackCopy->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]; - 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)); + AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); + if (v.Mag() < 0.0001) { + // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl; + delete trackCopy; + continue; + } + trackCopy->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->GetTPCInnerParam()); + 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); + } } - trackCopy->SetTrackId(esdtrack->GetID()); - trackCopy->SetFlags(esdtrack->GetStatus()); - trackCopy->SetLabel(esdtrack->GetLabel()); - trackCopy->SetITSchi2(esdtrack->GetITSchi2()); - trackCopy->SetITSncls(esdtrack->GetNcls(0)); - trackCopy->SetTPCchi2(esdtrack->GetTPCchi2()); - trackCopy->SetTPCncls(esdtrack->GetTPCNcls()); - trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF()); - trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class - trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); + 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 + } + + if (fRotateToEventPlane) { + double tPhi = TMath::ATan2(pxyz[1], pxyz[0]); + double tRad = TMath::Hypot(pxyz[0], pxyz[1]); + pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane); + pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane); + } - trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap()); - trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap()); + AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); + if (v.Mag() < 0.0001) { + // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl; + delete trackCopy; + continue; + } - double xtpc[3]; - esdtrack->GetInnerXYZ(xtpc); - xtpc[2] -= fV1[2]; - trackCopy->SetNominalTPCEntrancePoint(xtpc); + trackCopy->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]; + 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)); + } - esdtrack->GetOuterXYZ(xtpc); - xtpc[2] -= fV1[2]; - trackCopy->SetNominalTPCExitPoint(xtpc); + trackCopy->SetTrackId(esdtrack->GetID()); + trackCopy->SetFlags(esdtrack->GetStatus()); + trackCopy->SetLabel(esdtrack->GetLabel()); - int indexes[3]; - for (int ik=0; ik<3; ik++) { - indexes[ik] = esdtrack->GetKinkIndex(ik); - } - trackCopy->SetKinkIndexes(indexes); + trackCopy->SetITSchi2(esdtrack->GetITSchi2()); + trackCopy->SetITSncls(esdtrack->GetNcls(0)); + trackCopy->SetTPCchi2(esdtrack->GetTPCchi2()); + trackCopy->SetTPCncls(esdtrack->GetTPCNcls()); + trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF()); + trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class + trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); + trackCopy->SetTPCsignal(esdtrack->GetTPCsignal()); - for (int ii=0; ii<6; ii++){ - trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii)); - } + cout << "TPC signal = " << esdtrack->GetTPCsignal() << endl; + trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap()); + trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap()); - // Fill the hidden information with the simulated data - if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) { - TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel())); + double xtpc[3]; + esdtrack->GetInnerXYZ(xtpc); + xtpc[2] -= fV1[2]; + trackCopy->SetNominalTPCEntrancePoint(xtpc); - // Check the mother information + esdtrack->GetOuterXYZ(xtpc); + xtpc[2] -= fV1[2]; + trackCopy->SetNominalTPCExitPoint(xtpc); - // Using the new way of storing the freeze-out information - // Final state particle is stored twice on the stack - // one copy (mother) is stored with original freeze-out information - // and is not tracked - // the other one (daughter) is stored with primary vertex position - // and is tracked + int indexes[3]; + for (int ik=0; ik<3; ik++) { + indexes[ik] = esdtrack->GetKinkIndex(ik); + } + trackCopy->SetKinkIndexes(indexes); - // Freeze-out coordinates - double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; - fpx = tPart->Vx() - fV1[0]; - fpy = tPart->Vy() - fV1[1]; - fpz = tPart->Vz() - fV1[2]; - fpt = tPart->T(); + for (int ii=0; ii<6; ii++){ + trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii)); + } - AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); - tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); - trackCopy->SetGlobalEmissionPoint(fpx, fpy, fpz); - fpx *= 1e13; - fpy *= 1e13; - fpz *= 1e13; - fpt *= 1e13; + // Fill the hidden information with the simulated data + if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) { + TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel())); + // Check the mother information - if (fOnlyPrimaries && !fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){ - delete trackCopy; - continue; - } + // Using the new way of storing the freeze-out information + // Final state particle is stored twice on the stack + // one copy (mother) is stored with original freeze-out information + // and is not tracked + // the other one (daughter) is stored with primary vertex position + // and is tracked - // fillDCA + // Freeze-out coordinates + double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; + fpx = tPart->Vx() - fV1[0]; + fpy = tPart->Vy() - fV1[1]; + fpz = tPart->Vz() - fV1[2]; + fpt = tPart->T(); - if (TMath::Abs(impact[0]) > 0.001) { - if (fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){ - trackCopy->SetImpactDprim(impact[0]); - //cout << "prim" << endl; + AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); + tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); + trackCopy->SetGlobalEmissionPoint(fpx, fpy, fpz); - } - else if (fStack->IsSecondaryFromWeakDecay(TMath::Abs(esdtrack->GetLabel()))) { - trackCopy->SetImpactDweak(impact[0]); - //cout << "wea" << endl; - } - else if (fStack->IsSecondaryFromMaterial(TMath::Abs(esdtrack->GetLabel()))) { - trackCopy->SetImpactDmat(impact[0]); - //cout << "mat" << endl; - } - } - // end fillDCA - - // cout << "Looking for mother ids " << endl; - if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) { - // cout << "Got mother id" << endl; - TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]); - // Check if this is the same particle stored twice on the stack - if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) { - // It is the same particle - // Read in the original freeze-out information - // and convert it from to [fm] - - // EPOS style - fpx = mother->Vx()*1e13*0.197327; - fpy = mother->Vy()*1e13*0.197327; - fpz = mother->Vz()*1e13*0.197327; - fpt = mother->T() *1e13*0.197327; - - - // Therminator style + fpx *= 1e13; + fpy *= 1e13; + fpz *= 1e13; + fpt *= 1e13; + + + if (fOnlyPrimaries && !fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){ + delete trackCopy; + continue; + } + + // fillDCA + + if (TMath::Abs(impact[0]) > 0.001) { + if (fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){ + trackCopy->SetImpactDprim(impact[0]); + //cout << "prim" << endl; + + } + else if (fStack->IsSecondaryFromWeakDecay(TMath::Abs(esdtrack->GetLabel()))) { + trackCopy->SetImpactDweak(impact[0]); + //cout << "wea" << endl; + } + else if (fStack->IsSecondaryFromMaterial(TMath::Abs(esdtrack->GetLabel()))) { + trackCopy->SetImpactDmat(impact[0]); + //cout << "mat" << endl; + } + } + // end fillDCA + + // cout << "Looking for mother ids " << endl; + if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) { + // cout << "Got mother id" << endl; + TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]); + // Check if this is the same particle stored twice on the stack + if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) { + // It is the same particle + // Read in the original freeze-out information + // and convert it from to [fm] + + // EPOS style + fpx = mother->Vx()*1e13*0.197327; + fpy = mother->Vy()*1e13*0.197327; + fpz = mother->Vz()*1e13*0.197327; + fpt = mother->T() *1e13*0.197327; + + + // Therminator style // fpx = mother->Vx()*1e13; // fpy = mother->Vy()*1e13; // fpz = mother->Vz()*1e13; // fpt = mother->T() *1e13*3e10; - } - } - - if (fRotateToEventPlane) { - double tPhi = TMath::ATan2(fpy, fpx); - double tRad = TMath::Hypot(fpx, fpy); - - fpx = tRad*TMath::Cos(tPhi - tReactionPlane); - fpy = tRad*TMath::Sin(tPhi - tReactionPlane); - } - - tInfo->SetPDGPid(tPart->GetPdgCode()); - trackCopy->SetPDGPid(tPart->GetPdgCode()); - - if (fRotateToEventPlane) { - double tPhi = TMath::ATan2(tPart->Py(), tPart->Px()); - double tRad = TMath::Hypot(tPart->Px(), tPart->Py()); - - tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane), - tRad*TMath::Sin(tPhi - tReactionPlane), - tPart->Pz()); - trackCopy->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane), - tRad*TMath::Sin(tPhi - tReactionPlane), - tPart->Pz()); - } - else - { - tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); - trackCopy->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); - } - Double_t mass2 = (tPart->Energy() *tPart->Energy() - - tPart->Px()*tPart->Px() - - tPart->Py()*tPart->Py() - - tPart->Pz()*tPart->Pz()); - if (mass2>0.0) - { - tInfo->SetMass(TMath::Sqrt(mass2)); - trackCopy->SetMass(TMath::Sqrt(mass2)); - } - else - { - tInfo->SetMass(0.0); - trackCopy->SetMass(0.0); - } - - tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); - trackCopy->SetEmissionPoint(fpx, fpy, fpz, fpt); - trackCopy->SetHiddenInfo(tInfo); - } - else { - AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); - tInfo->SetMass(0.0); - double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; - fpx = fV1[0]*1e13; - fpy = fV1[1]*1e13; - fpz = fV1[2]*1e13; - fpt = 0.0; - tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); - - tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]); - - trackCopy->SetHiddenInfo(tInfo); + } } - // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl; - //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 (fRotateToEventPlane) { + double tPhi = TMath::ATan2(fpy, fpx); + double tRad = TMath::Hypot(fpx, fpy); - bool trackAccept = true; - if (isKaonAnalysis == true && trackCopy->GetPDGPid() != 321) { - trackAccept = false; + fpx = tRad*TMath::Cos(tPhi - tReactionPlane); + fpy = tRad*TMath::Sin(tPhi - tReactionPlane); } - if (tGoodMomentum==true && trackAccept == true) - { + tInfo->SetPDGPid(tPart->GetPdgCode()); + trackCopy->SetPDGPid(tPart->GetPdgCode()); - hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis - realnofTracks++;//real number of tracks - // delete trackCopy; - } + if (fRotateToEventPlane) { + double tPhi = TMath::ATan2(tPart->Py(), tPart->Px()); + double tRad = TMath::Hypot(tPart->Px(), tPart->Py()); + + tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane), + tRad*TMath::Sin(tPhi - tReactionPlane), + tPart->Pz()); + trackCopy->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane), + tRad*TMath::Sin(tPhi - tReactionPlane), + tPart->Pz()); + } + else + { + tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); + trackCopy->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz()); + } + Double_t mass2 = (tPart->Energy() *tPart->Energy() - + tPart->Px()*tPart->Px() - + tPart->Py()*tPart->Py() - + tPart->Pz()*tPart->Pz()); + if (mass2>0.0) + { + tInfo->SetMass(TMath::Sqrt(mass2)); + trackCopy->SetMass(TMath::Sqrt(mass2)); + } else - { - delete trackCopy; - } + { + tInfo->SetMass(0.0); + trackCopy->SetMass(0.0); + } + + tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); + trackCopy->SetEmissionPoint(fpx, fpy, fpz, fpt); + trackCopy->SetHiddenInfo(tInfo); + } + else { + AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); + tInfo->SetMass(0.0); + double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; + fpx = fV1[0]*1e13; + fpy = fV1[1]*1e13; + fpz = fV1[2]*1e13; + fpt = 0.0; + tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); + + tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]); + + trackCopy->SetHiddenInfo(tInfo); + } + // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl; + + //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 + + bool trackAccept = true; + if (isKaonAnalysis == true && trackCopy->GetPDGPid() != 321) { + trackAccept = false; + } + if (tGoodMomentum==true && trackAccept == true) + { + + hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis + realnofTracks++;//real number of tracks + // delete trackCopy; + } + else + { + delete trackCopy; } + } + delete [] motherids; hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event @@ -867,22 +869,22 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() 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); - } + { + 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); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M"))); + 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"))); + 10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M"))); } } else if (fEstEventMult == kCentralityV0A) { @@ -890,11 +892,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("V0A") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A"))); + 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"))); + 10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A"))); } } else if (fEstEventMult == kCentralityV0C) { @@ -902,11 +904,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("V0C") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C"))); + 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"))); + 10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C"))); } } else if (fEstEventMult == kCentralityZNA) { @@ -914,11 +916,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("ZNA") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA"))); + 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"))); + 10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA"))); } } else if (fEstEventMult == kCentralityZNC) { @@ -926,11 +928,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("ZNC") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC"))); + 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"))); + 10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC"))); } } else if (fEstEventMult == kCentralityCL1) { @@ -938,11 +940,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("CL1") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1"))); + 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"))); + 10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1"))); } } else if (fEstEventMult == kCentralityCL0) { @@ -950,11 +952,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("CL0") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0"))); + 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"))); + 10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0"))); } } else if (fEstEventMult == kCentralityTRK) { @@ -962,11 +964,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("TRK") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK"))); + 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"))); + 10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK"))); } } else if (fEstEventMult == kCentralityTKL) { @@ -974,11 +976,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("TKL") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL"))); + 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"))); + 10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL"))); } } else if (fEstEventMult == kCentralityNPA) { @@ -986,11 +988,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("NPA") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA"))); + 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"))); + 10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA"))); } } else if (fEstEventMult == kCentralityCND) { @@ -998,11 +1000,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("CND") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND"))); + 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"))); + 10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND"))); } } else if (fEstEventMult == kCentralityFMD) { @@ -1010,11 +1012,11 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (cent) { if (cent->GetCentralityPercentile("FMD") < 0.00001) - hbtEvent->SetNormalizedMult(-1); + hbtEvent->SetNormalizedMult(-1); else - hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD"))); + 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"))); + 10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD"))); } } @@ -1031,24 +1033,24 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } 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"<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"<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[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][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[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][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++) { - 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[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][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[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][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()); + 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(); - 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); - } + 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) + 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); @@ -1296,128 +1298,128 @@ void AliFemtoEventReaderESDChainKine::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemt 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); - } - } - - if ((TMath::Abs(trackpos->GetLabel()) < fStack->GetNtrack()) && (TMath::Abs(trackneg->GetLabel()) < fStack->GetNtrack())) { - - //cout<<"tESDv0->GetPdgCode(): "<GetPdgCode()<GetLabel()<<" "<GetLabel()<Particle(tESDv0->GetLabel()); //zle + } + 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); + } + } - int labelpos = TMath::Abs(trackpos->GetLabel()); - int labelneg = TMath::Abs(trackneg->GetLabel()); - TParticle *tPartPos = fStack->Particle(labelpos); - TParticle *tPartNeg = fStack->Particle(labelneg); + if ((TMath::Abs(trackpos->GetLabel()) < fStack->GetNtrack()) && (TMath::Abs(trackneg->GetLabel()) < fStack->GetNtrack())) { + //cout<<"tESDv0->GetPdgCode(): "<GetPdgCode()<GetLabel()<<" "<GetLabel()<Particle(tESDv0->GetLabel()); //zle - double impactpos[2]; - double impactneg[2]; - double covimpact[3]; + int labelpos = TMath::Abs(trackpos->GetLabel()); + int labelneg = TMath::Abs(trackneg->GetLabel()); + TParticle *tPartPos = fStack->Particle(labelpos); + TParticle *tPartNeg = fStack->Particle(labelneg); - AliExternalTrackParam *parampos = new AliExternalTrackParam(*trackpos->GetTPCInnerParam()); - parampos->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactpos, covimpact); - AliExternalTrackParam *paramneg = new AliExternalTrackParam(*trackneg->GetTPCInnerParam()); - paramneg->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactneg, covimpact); + double impactpos[2]; + double impactneg[2]; + double covimpact[3]; - // fillDCA - if (TMath::Abs(impactpos[0]) > 0.001) { - if (fStack->IsPhysicalPrimary(labelpos)){ - tFemtoV0->SetImpactDprimPos(impactpos[0]); - } - else if (fStack->IsSecondaryFromWeakDecay(labelpos)) { - tFemtoV0->SetImpactDweakPos(impactpos[0]); - //cout << "wea" << endl; - } - else if (fStack->IsSecondaryFromMaterial(labelpos)) { - tFemtoV0->SetImpactDmatPos(impactpos[0]); - //cout << "mat" << endl; - } - } - if (TMath::Abs(impactneg[0]) > 0.001) { - if (fStack->IsPhysicalPrimary(labelneg)){ - tFemtoV0->SetImpactDprimNeg(impactneg[0]); - //cout << "prim" << endl; - } - else if (fStack->IsSecondaryFromWeakDecay(labelneg)) { - tFemtoV0->SetImpactDweakNeg(impactneg[0]); - //cout << "wea" << endl; - } - else if (fStack->IsSecondaryFromMaterial(labelneg)) { - tFemtoV0->SetImpactDmatNeg(impactneg[0]); - //cout << "mat" << endl; - } + AliExternalTrackParam *parampos = new AliExternalTrackParam(*trackpos->GetTPCInnerParam()); + parampos->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactpos, covimpact); + AliExternalTrackParam *paramneg = new AliExternalTrackParam(*trackneg->GetTPCInnerParam()); + paramneg->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactneg, covimpact); + + + // fillDCA + if (TMath::Abs(impactpos[0]) > 0.001) { + if (fStack->IsPhysicalPrimary(labelpos)){ + tFemtoV0->SetImpactDprimPos(impactpos[0]); + } + else if (fStack->IsSecondaryFromWeakDecay(labelpos)) { + tFemtoV0->SetImpactDweakPos(impactpos[0]); + //cout << "wea" << endl; + } + else if (fStack->IsSecondaryFromMaterial(labelpos)) { + tFemtoV0->SetImpactDmatPos(impactpos[0]); + //cout << "mat" << endl; + } + } + if (TMath::Abs(impactneg[0]) > 0.001) { + if (fStack->IsPhysicalPrimary(labelneg)){ + tFemtoV0->SetImpactDprimNeg(impactneg[0]); + //cout << "prim" << endl; + } + else if (fStack->IsSecondaryFromWeakDecay(labelneg)) { + tFemtoV0->SetImpactDweakNeg(impactneg[0]); + //cout << "wea" << endl; + } + else if (fStack->IsSecondaryFromMaterial(labelneg)) { + tFemtoV0->SetImpactDmatNeg(impactneg[0]); + //cout << "mat" << endl; + } - } - // end fillDCA + } + // end fillDCA - //tInfo->SetPDGPid(); - //tInfo->SetMass(); - //tInfo->SetTrueMomentum(); - //tInfo->SetEmissionPoint(); + //tInfo->SetPDGPid(); + //tInfo->SetMass(); + //tInfo->SetTrueMomentum(); + //tInfo->SetEmissionPoint(); - // Freeze-out coordinates - double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; + // Freeze-out coordinates + double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; - fpx = tPartPos->Vx() - fPrimaryVtxPosition[0]; - fpy = tPartPos->Vy() - fPrimaryVtxPosition[1]; - fpz = tPartPos->Vz() - fPrimaryVtxPosition[2]; - fpt = tPartPos->T(); + fpx = tPartPos->Vx() - fPrimaryVtxPosition[0]; + fpy = tPartPos->Vy() - fPrimaryVtxPosition[1]; + fpz = tPartPos->Vz() - fPrimaryVtxPosition[2]; + fpt = tPartPos->T(); - fpx *= 1e13; - fpy *= 1e13; - fpz *= 1e13; - fpt *= 1e13; + fpx *= 1e13; + fpy *= 1e13; + fpz *= 1e13; + fpt *= 1e13; - tInfo->SetPDGPidPos(tPartPos->GetPdgCode()); - tInfo->SetMassPos(tPartPos->GetMass()); - tInfo->SetTrueMomentumPos(tPartPos->Px(),tPartPos->Py(),tPartPos->Pz()); - tInfo->SetEmissionPointPos(fpx,fpy,fpz,fpt); + tInfo->SetPDGPidPos(tPartPos->GetPdgCode()); + tInfo->SetMassPos(tPartPos->GetMass()); + tInfo->SetTrueMomentumPos(tPartPos->Px(),tPartPos->Py(),tPartPos->Pz()); + tInfo->SetEmissionPointPos(fpx,fpy,fpz,fpt); - fpx = tPartNeg->Vx() - fPrimaryVtxPosition[0]; - fpy = tPartNeg->Vy() - fPrimaryVtxPosition[1]; - fpz = tPartNeg->Vz() - fPrimaryVtxPosition[2]; - fpt = tPartNeg->T(); + fpx = tPartNeg->Vx() - fPrimaryVtxPosition[0]; + fpy = tPartNeg->Vy() - fPrimaryVtxPosition[1]; + fpz = tPartNeg->Vz() - fPrimaryVtxPosition[2]; + fpt = tPartNeg->T(); - fpx *= 1e13; - fpy *= 1e13; - fpz *= 1e13; - fpt *= 1e13; + fpx *= 1e13; + fpy *= 1e13; + fpz *= 1e13; + fpt *= 1e13; - tInfo->SetPDGPidNeg(tPartNeg->GetPdgCode()); - tInfo->SetMassNeg(tPartNeg->GetMass()); - tInfo->SetTrueMomentumNeg(tPartNeg->Px(),tPartNeg->Py(),tPartNeg->Pz()); - tInfo->SetEmissionPointNeg(fpx,fpy,fpz,fpt); + tInfo->SetPDGPidNeg(tPartNeg->GetPdgCode()); + tInfo->SetMassNeg(tPartNeg->GetMass()); + tInfo->SetTrueMomentumNeg(tPartNeg->Px(),tPartNeg->Py(),tPartNeg->Pz()); + tInfo->SetEmissionPointNeg(fpx,fpy,fpz,fpt); - tFemtoV0->SetHiddenInfo(tInfo); - } + tFemtoV0->SetHiddenInfo(tInfo); } + } else - { - tFemtoV0->SetStatusPos(999); - tFemtoV0->SetStatusNeg(999); - } + { + tFemtoV0->SetStatusPos(999); + tFemtoV0->SetStatusNeg(999); + } tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus()); } @@ -1479,11 +1481,11 @@ void AliFemtoEventReaderESDChainKine::GetGlobalPositionAtGlobalRadiiThroughTPC(A // 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 + 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]; -- 2.43.0