X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG2%2FFEMTOSCOPY%2FAliFemto%2FAliFemtoEventReaderESDChainKine.cxx;h=a49f97db5d1467fbd66a6a5404a1d822948dbe63;hb=b84aaa8bf0ffd8228a9f1b6d35ba9851227dd05d;hp=b0ce8f04406784be45cb3c04ef06bd85c661f7d5;hpb=63a5982ace8ce90682ea7a8075d2c2b9830293a5;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx index b0ce8f04406..a49f97db5d1 100644 --- a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx +++ b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx @@ -10,10 +10,15 @@ #include "TFile.h" #include "TTree.h" +#include "TList.h" #include "AliESDEvent.h" #include "AliESDtrack.h" #include "AliESDVertex.h" +#include "AliMultiplicity.h" +#include "AliCentrality.h" +#include "AliEventplane.h" + #include "AliFmPhysicalHelixD.h" #include "AliFmThreeVectorF.h" @@ -23,7 +28,13 @@ #include "TParticle.h" #include "AliFemtoModelHiddenInfo.h" +#include "AliFemtoModelGlobalHiddenInfo.h" #include "AliGenHijingEventHeader.h" +#include "AliGenCocktailEventHeader.h" + +#include "AliVertexerTracks.h" + +#include "AliPID.h" ClassImp(AliFemtoEventReaderESDChainKine) @@ -36,13 +47,20 @@ using namespace std; AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(): fFileName(" "), fConstrained(true), + fReadInner(false), fUseTPCOnly(false), fNumberofEvent(0), fCurEvent(0), fCurFile(0), fEvent(0x0), fStack(0x0), - fGenHeader(0x0) + fGenHeader(0x0), + fTrackType(kGlobal), + fEstEventMult(kV0Centrality), + fRotateToEventPlane(0), + fESDpid(0), + fIsPidOwner(0) + { //constructor with 0 parameters , look at default settings } @@ -52,22 +70,34 @@ AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoE AliFemtoEventReader(aReader), fFileName(" "), fConstrained(true), + fReadInner(false), fUseTPCOnly(false), fNumberofEvent(0), fCurEvent(0), fCurFile(0), fEvent(0x0), fStack(0x0), - fGenHeader(0x0) + fGenHeader(0x0), + fTrackType(kGlobal), + fEstEventMult(kV0Centrality), + fRotateToEventPlane(0), + fESDpid(0), + fIsPidOwner(0) { // Copy constructor fConstrained = aReader.fConstrained; + fReadInner = aReader.fReadInner; fUseTPCOnly = aReader.fUseTPCOnly; fNumberofEvent = aReader.fNumberofEvent; fCurEvent = aReader.fCurEvent; fCurFile = aReader.fCurFile; fEvent = new AliESDEvent(); fStack = aReader.fStack; + fTrackType = aReader.fTrackType; + fEstEventMult = aReader.fEstEventMult; + fRotateToEventPlane = aReader.fRotateToEventPlane; + fESDpid = aReader.fESDpid; + fIsPidOwner = aReader.fIsPidOwner; } //__________________ AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine() @@ -92,6 +122,9 @@ AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(cons fEvent = new AliESDEvent(); fStack = aReader.fStack; fGenHeader = aReader.fGenHeader; + fRotateToEventPlane = aReader.fRotateToEventPlane; + fESDpid = aReader.fESDpid; + fIsPidOwner = aReader.fIsPidOwner; return *this; } @@ -116,6 +149,31 @@ bool AliFemtoEventReaderESDChainKine::GetConstrained() const return fConstrained; } //__________________ +void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner) +{ + fReadInner=readinner; +} + +bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const +{ + return fReadInner; +} + +//__________________ +void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly) +{ + fUseTPCOnly=usetpconly; +} + +bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const +{ + return fUseTPCOnly; +} +void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType) +{ + fEstEventMult = aType; +} +//__________________ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() { // Get the event, read all the relevant information @@ -125,7 +183,7 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() string tFriendFileName; // Get the friend information - cout<<"starting to read event "<SetESDfriend(fEventFriend); hbtEvent = new AliFemtoEvent; @@ -141,8 +199,19 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy()); hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants()); hbtEvent->SetTriggerMask(fEvent->GetTriggerMask()); - hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster()); - + //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]; @@ -160,35 +229,77 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]); - hbtEvent->SetPrimVertPos(vertex); hbtEvent->SetPrimVertCov(fVCov); - AliGenHijingEventHeader *hdh = dynamic_cast (fGenHeader); - Double_t tReactionPlane = 0; + + AliGenHijingEventHeader *hdh = dynamic_cast (fGenHeader); + if (!hdh) { + AliGenCocktailEventHeader *cdh = dynamic_cast (fGenHeader); + if (cdh) { + TList *tGenHeaders = cdh->GetHeaders(); + for (int ihead = 0; iheadGetEntries(); ihead++) { + hdh = dynamic_cast (fGenHeader); + if (hdh) break; + } + } + } + if (hdh) { tReactionPlane = hdh->ReactionPlaneAngle(); + cout << "Got reaction plane " << tReactionPlane << endl; } + + hbtEvent->SetReactionPlaneAngle(tReactionPlane); + + 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 - Int_t *motherids; - motherids = new Int_t[fStack->GetNtrack()]; - for (int ip=0; ipGetNtrack(); ip++) motherids[ip] = 0; - - // Read in mother ids - TParticle *motherpart; - for (int ip=0; ipGetNtrack(); ip++) { - motherpart = fStack->Particle(ip); - if (motherpart->GetDaughter(0) > 0) - motherids[motherpart->GetDaughter(0)] = ip; - if (motherpart->GetDaughter(1) > 0) - motherids[motherpart->GetDaughter(1)] = ip; + 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, tITSPure=0; + + fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2); + + hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC); + hbtEvent->SetMultiplicityEstimateTracklets(tTracklet); + // hbtEvent->SetMultiplicityEstimateITSPure(tITSPure); + hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1)); + Int_t *motherids; + if (fStack) { + motherids = new Int_t[fStack->GetNtrack()]; + for (int ip=0; ipGetNtrack(); ip++) motherids[ip] = 0; + + // Read in mother ids + TParticle *motherpart; + for (int ip=0; ipGetNtrack(); ip++) { + motherpart = fStack->Particle(ip); + if (motherpart->GetDaughter(0) > 0) + motherids[motherpart->GetDaughter(0)] = ip; + if (motherpart->GetDaughter(1) > 0) + motherids[motherpart->GetDaughter(1)] = ip; + // if (motherpart->GetPdgCode() == 211) { // cout << "Mother " << ip << " has daughters " // << motherpart->GetDaughter(0) << " " @@ -199,16 +310,80 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() // << endl; // } + } + } + else { + printf ("No Stack ???\n"); + delete hbtEvent; + return 0; } 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(); + } + } + } + } + + 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 @@ -221,6 +396,76 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() 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::kTOFpid) { + 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) && + !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch)) + { + + //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) && + //(esdtrack->GetStatus()&AliESDtrack::kTOFout) && + //(esdtrack->GetStatus()&AliESDtrack::kTIME)){ + // collect info from ESDpid class + + if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) { + + + 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]; @@ -229,6 +474,7 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() if (fUseTPCOnly) { if (!esdtrack->GetTPCInnerParam()) { + cout << "No TPC inner param !" << endl; delete trackCopy; continue; } @@ -238,9 +484,28 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() 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 (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 ???? " <GetTPCInnerParam()) { + AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam()); + 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 + } + + 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 ???? " <SetP(v);//setting momentum trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]); @@ -333,54 +644,108 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } trackCopy->SetKinkIndexes(indexes); - // Fill the hidden information with the simulated data - TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel())); - // Check the mother information + // Fill the hidden information with the simulated data + if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) { + TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel())); - // 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 + // Check the mother information + + // 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 + + // 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(); + + AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); + tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz); - // Freeze-out coordinates - double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0; - fpx = tPart->Vx(); - fpy = tPart->Vy(); - fpz = tPart->Vz(); - fpt = tPart->T(); - - if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) { - 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] - fpx = mother->Vx()*1e13; - fpy = mother->Vy()*1e13; - fpz = mother->Vz()*1e13; - fpt = mother->T()*1e13*3e10; + fpx *= 1e13; + fpy *= 1e13; + fpz *= 1e13; + fpt *= 1e13; + + // 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()); + + 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()); } + else + tInfo->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)); + else + tInfo->SetMass(0.0); + + tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); + trackCopy->SetHiddenInfo(tInfo); } - - AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo(); - tInfo->SetPDGPid(tPart->GetPdgCode()); - tInfo->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)); - else + else { + AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo(); tInfo->SetMass(0.0); - - tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt); - trackCopy->SetHiddenInfo(tInfo); + 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 @@ -397,10 +762,58 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent() } } + + delete [] motherids; hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event + + AliCentrality *cent = fEvent->GetCentrality(); + if (cent) { + hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M")); + // hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD")); + hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1")); + // hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK")); + + 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 == kTracklet) + hbtEvent->SetNormalizedMult(tTracklet); + else if (fEstEventMult == kITSTPC) + hbtEvent->SetNormalizedMult(tITSTPC); + else if (fEstEventMult == kITSPure) + hbtEvent->SetNormalizedMult(tITSPure); + else if (fEstEventMult == kSPDLayer1) + hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1)); + else if (fEstEventMult == kV0Centrality) { + // 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"))); + } + } + + if (tNormMultPos > tNormMultNeg) + hbtEvent->SetZDCParticipants(tNormMultPos); + else + hbtEvent->SetZDCParticipants(tNormMultNeg); + + AliEventplane* ep = fEvent->GetEventplane(); + if (ep) { + hbtEvent->SetEP(ep); + hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q")); + } + + fCurEvent++; - cout<<"end of reading nt "<