1 ////////////////////////////////////////////////////////////////////////////////
3 // AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and //
4 // the model Kinematics information tailored for the Task framework and the //
5 // Reads in AliESDfriend to create shared hit/quality information //
6 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
8 ////////////////////////////////////////////////////////////////////////////////
9 #include "AliFemtoEventReaderESDChainKine.h"
14 #include "AliESDEvent.h"
15 #include "AliESDtrack.h"
16 #include "AliESDVertex.h"
18 #include "AliMultiplicity.h"
19 #include "AliCentrality.h"
20 #include "AliEventplane.h"
22 #include "AliFmPhysicalHelixD.h"
23 #include "AliFmThreeVectorF.h"
25 #include "SystemOfUnits.h"
27 #include "AliFemtoEvent.h"
29 #include "TParticle.h"
30 #include "AliFemtoModelHiddenInfo.h"
31 #include "AliFemtoModelGlobalHiddenInfo.h"
32 #include "AliGenHijingEventHeader.h"
33 #include "AliGenCocktailEventHeader.h"
35 #include "AliVertexerTracks.h"
39 ClassImp(AliFemtoEventReaderESDChainKine)
41 #if !(ST_NO_NAMESPACES)
42 using namespace units;
46 //____________________________
47 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
59 fEstEventMult(kV0Centrality),
60 fRotateToEventPlane(0),
65 //constructor with 0 parameters , look at default settings
69 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
70 AliFemtoEventReader(aReader),
82 fEstEventMult(kV0Centrality),
83 fRotateToEventPlane(0),
88 fConstrained = aReader.fConstrained;
89 fReadInner = aReader.fReadInner;
90 fUseTPCOnly = aReader.fUseTPCOnly;
91 fNumberofEvent = aReader.fNumberofEvent;
92 fCurEvent = aReader.fCurEvent;
93 fCurFile = aReader.fCurFile;
94 fEvent = new AliESDEvent();
95 fStack = aReader.fStack;
96 fTrackType = aReader.fTrackType;
97 fEstEventMult = aReader.fEstEventMult;
98 fRotateToEventPlane = aReader.fRotateToEventPlane;
99 fESDpid = aReader.fESDpid;
100 fIsPidOwner = aReader.fIsPidOwner;
103 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
110 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
112 // Assignment operator
113 if (this == &aReader)
116 fConstrained = aReader.fConstrained;
117 fUseTPCOnly = aReader.fUseTPCOnly;
118 fNumberofEvent = aReader.fNumberofEvent;
119 fCurEvent = aReader.fCurEvent;
120 fCurFile = aReader.fCurFile;
121 if (fEvent) delete fEvent;
122 fEvent = new AliESDEvent();
123 fStack = aReader.fStack;
124 fGenHeader = aReader.fGenHeader;
125 fRotateToEventPlane = aReader.fRotateToEventPlane;
126 fESDpid = aReader.fESDpid;
127 fIsPidOwner = aReader.fIsPidOwner;
133 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
135 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
140 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
142 // Select whether to read constrained or not constrained momentum
143 fConstrained=constrained;
146 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
148 // Check whether we read constrained or not constrained momentum
152 void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
154 fReadInner=readinner;
157 bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
163 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
165 fUseTPCOnly=usetpconly;
168 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
172 void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
174 fEstEventMult = aType;
177 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
179 // Get the event, read all the relevant information
180 // and fill the AliFemtoEvent class
181 // Returns a valid AliFemtoEvent
182 AliFemtoEvent *hbtEvent = 0;
183 string tFriendFileName;
185 // Get the friend information
186 cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
187 // fEvent->SetESDfriend(fEventFriend);
189 hbtEvent = new AliFemtoEvent;
190 //setting basic things
191 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
192 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
193 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
194 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
195 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
196 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
197 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
198 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
199 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
200 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
201 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
202 //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
204 if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
205 (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
206 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
207 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
208 hbtEvent->SetTriggerCluster(1);
209 else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
210 (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
211 hbtEvent->SetTriggerCluster(2);
213 hbtEvent->SetTriggerCluster(0);
219 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
220 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
221 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
225 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
226 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
227 if (!fEvent->GetPrimaryVertex()->GetStatus())
231 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
232 hbtEvent->SetPrimVertPos(vertex);
233 hbtEvent->SetPrimVertCov(fVCov);
235 Double_t tReactionPlane = 0;
237 AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
239 AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
241 TList *tGenHeaders = cdh->GetHeaders();
242 for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
243 hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
251 tReactionPlane = hdh->ReactionPlaneAngle();
252 cout << "Got reaction plane " << tReactionPlane << endl;
255 hbtEvent->SetReactionPlaneAngle(tReactionPlane);
257 Int_t spdetaonecount = 0;
259 for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
260 if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
263 // hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
264 hbtEvent->SetSPDMult(spdetaonecount);
266 //starting to reading tracks
267 int nofTracks=0; //number of reconstructed tracks in event
268 nofTracks=fEvent->GetNumberOfTracks();
269 int realnofTracks=0;//number of track which we use ina analysis
272 int tNormMultPos = 0;
273 int tNormMultNeg = 0;
275 Float_t tTotalPt = 0.0;
280 Int_t tTracklet=0, tITSTPC=0, tITSPure=0;
282 fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
284 hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
285 hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
286 // hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
287 hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
291 motherids = new Int_t[fStack->GetNtrack()];
292 for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
294 // Read in mother ids
295 TParticle *motherpart;
296 for (int ip=0; ip<fStack->GetNtrack(); ip++) {
297 motherpart = fStack->Particle(ip);
298 if (motherpart->GetDaughter(0) > 0)
299 motherids[motherpart->GetDaughter(0)] = ip;
300 if (motherpart->GetDaughter(1) > 0)
301 motherids[motherpart->GetDaughter(1)] = ip;
303 // if (motherpart->GetPdgCode() == 211) {
304 // cout << "Mother " << ip << " has daughters "
305 // << motherpart->GetDaughter(0) << " "
306 // << motherpart->GetDaughter(1) << " "
307 // << motherpart->Vx() << " "
308 // << motherpart->Vy() << " "
309 // << motherpart->Vz() << " "
316 printf ("No Stack ???\n");
321 for (int i=0;i<nofTracks;i++)
323 // cout << "Reading track " << i << endl;
324 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
327 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
328 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
331 if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
332 (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
333 if (esdtrack->GetTPCNcls() > 70)
334 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
335 if (TMath::Abs(esdtrack->Eta()) < 1.2) {
336 esdtrack->GetImpactParameters(b,bCov);
337 if ((b[0]<0.2) && (b[1] < 0.25)) {
339 tTotalPt += esdtrack->Pt();
344 else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
345 if (esdtrack->GetTPCNcls() > 100)
346 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
347 if (TMath::Abs(esdtrack->Eta()) < 1.2) {
348 esdtrack->GetImpactParameters(b,bCov);
349 if ((b[0]<2.4) && (b[1] < 3.2)) {
351 tTotalPt += esdtrack->Pt();
357 hbtEvent->SetZDCEMEnergy(tTotalPt);
358 // if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
359 // if (esdtrack->GetTPCNcls() > 80)
360 // if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
361 // if (esdtrack->GetConstrainedParam())
362 // if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
363 // if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
364 // if (esdtrack->GetSign() > 0)
366 // else if (esdtrack->GetSign() < 0)
370 // If reading ITS-only tracks, reject all with TPC
371 if (fTrackType == kITSOnly) {
372 if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
373 if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
374 if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
375 UChar_t iclm = esdtrack->GetITSClusterMap();
377 for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
379 if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
386 AliFemtoTrack* trackCopy = new AliFemtoTrack();
387 trackCopy->SetCharge((short)esdtrack->GetSign());
389 //in aliroot we have AliPID
390 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
391 //we use only 5 first
393 esdtrack->GetESDpid(esdpid);
394 trackCopy->SetPidProbElectron(esdpid[0]);
395 trackCopy->SetPidProbMuon(esdpid[1]);
396 trackCopy->SetPidProbPion(esdpid[2]);
397 trackCopy->SetPidProbKaon(esdpid[3]);
398 trackCopy->SetPidProbProton(esdpid[4]);
400 esdpid[0] = -100000.0;
401 esdpid[1] = -100000.0;
402 esdpid[2] = -100000.0;
403 esdpid[3] = -100000.0;
404 esdpid[4] = -100000.0;
408 if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
409 tTOF = esdtrack->GetTOFsignal();
410 esdtrack->GetIntegratedTimes(esdpid);
413 trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
416 ////// TPC ////////////////////////////////////////////
417 float nsigmaTPCK=-1000.;
418 float nsigmaTPCPi=-1000.;
419 float nsigmaTPCP=-1000.;
421 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
422 nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
423 nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
424 nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
427 trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
428 trackCopy->SetNSigmaTPCK(nsigmaTPCK);
429 trackCopy->SetNSigmaTPCP(nsigmaTPCP);
431 ///// TOF ///////////////////////////////////////////////
434 float nsigmaTOFPi=-1000.;
435 float nsigmaTOFK=-1000.;
436 float nsigmaTOFP=-1000.;
438 if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
439 (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
440 (esdtrack->GetStatus()&AliESDtrack::kTIME) &&
441 !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch))
444 //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
445 //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
446 //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
447 // collect info from ESDpid class
449 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
452 double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
454 nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
455 nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
456 nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
458 Double_t len=esdtrack->GetIntegratedLength();
459 Double_t tof=esdtrack->GetTOFsignal();
460 if(tof > 0.) vp=len/tof/0.03;
464 trackCopy->SetVTOF(vp);
465 trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
466 trackCopy->SetNSigmaTOFK(nsigmaTOFK);
467 trackCopy->SetNSigmaTOFP(nsigmaTOFP);
476 if (!esdtrack->GetTPCInnerParam()) {
477 cout << "No TPC inner param !" << endl;
482 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
484 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
485 param->GetPxPyPz(pxyz);//reading noconstarined momentum
487 if (fReadInner == true) {
488 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
489 tInfo->SetPDGPid(211);
490 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
491 tInfo->SetMass(0.13957);
492 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
493 // tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
494 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
495 trackCopy->SetHiddenInfo(tInfo);
498 if (fRotateToEventPlane) {
499 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
500 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
502 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
503 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
506 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
507 if (v.Mag() < 0.0001) {
508 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
512 trackCopy->SetP(v);//setting momentum
513 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
515 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
516 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
517 //setting helix I do not if it is ok
518 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
519 trackCopy->SetHelix(helix);
521 //some stuff which could be useful
522 trackCopy->SetImpactD(impact[0]);
523 trackCopy->SetImpactZ(impact[1]);
524 trackCopy->SetCdd(covimpact[0]);
525 trackCopy->SetCdz(covimpact[1]);
526 trackCopy->SetCzz(covimpact[2]);
527 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
532 if (fReadInner == true) {
534 if (esdtrack->GetTPCInnerParam()) {
535 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
537 // param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
538 param->GetPxPyPz(pxyz);//reading noconstarined momentum
541 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
542 tInfo->SetPDGPid(211);
543 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
544 tInfo->SetMass(0.13957);
545 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
546 //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
547 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
548 trackCopy->SetHiddenInfo(tInfo);
553 if(fTrackType == kGlobal) {
554 if (fConstrained==true)
555 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
557 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
559 else if (fTrackType == kTPCOnly) {
560 if (esdtrack->GetTPCInnerParam())
561 esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
567 else if (fTrackType == kITSOnly) {
568 if (fConstrained==true)
569 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
571 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
574 if (fRotateToEventPlane) {
575 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
576 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
578 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
579 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
582 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
583 if (v.Mag() < 0.0001) {
584 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
589 trackCopy->SetP(v);//setting momentum
590 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
591 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
592 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
593 //setting helix I do not if it is ok
594 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
595 trackCopy->SetHelix(helix);
597 //some stuff which could be useful
600 esdtrack->GetImpactParameters(imp,cim);
604 covimpact[0] = cim[0];
605 covimpact[1] = cim[1];
606 covimpact[2] = cim[2];
608 trackCopy->SetImpactD(impact[0]);
609 trackCopy->SetImpactZ(impact[1]);
610 trackCopy->SetCdd(covimpact[0]);
611 trackCopy->SetCdz(covimpact[1]);
612 trackCopy->SetCzz(covimpact[2]);
613 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
616 trackCopy->SetTrackId(esdtrack->GetID());
617 trackCopy->SetFlags(esdtrack->GetStatus());
618 trackCopy->SetLabel(esdtrack->GetLabel());
620 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
621 trackCopy->SetITSncls(esdtrack->GetNcls(0));
622 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
623 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
624 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
625 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
626 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
629 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
630 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
633 esdtrack->GetInnerXYZ(xtpc);
635 trackCopy->SetNominalTPCEntrancePoint(xtpc);
637 esdtrack->GetOuterXYZ(xtpc);
639 trackCopy->SetNominalTPCExitPoint(xtpc);
642 for (int ik=0; ik<3; ik++) {
643 indexes[ik] = esdtrack->GetKinkIndex(ik);
645 trackCopy->SetKinkIndexes(indexes);
648 // Fill the hidden information with the simulated data
649 if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
650 TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
652 // Check the mother information
654 // Using the new way of storing the freeze-out information
655 // Final state particle is stored twice on the stack
656 // one copy (mother) is stored with original freeze-out information
657 // and is not tracked
658 // the other one (daughter) is stored with primary vertex position
661 // Freeze-out coordinates
662 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
663 fpx = tPart->Vx() - fV1[0];
664 fpy = tPart->Vy() - fV1[1];
665 fpz = tPart->Vz() - fV1[2];
668 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
669 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
676 // cout << "Looking for mother ids " << endl;
677 if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
678 // cout << "Got mother id" << endl;
679 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
680 // Check if this is the same particle stored twice on the stack
681 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
682 // It is the same particle
683 // Read in the original freeze-out information
684 // and convert it from to [fm]
687 fpx = mother->Vx()*1e13*0.197327;
688 fpy = mother->Vy()*1e13*0.197327;
689 fpz = mother->Vz()*1e13*0.197327;
690 fpt = mother->T() *1e13*0.197327;
694 // fpx = mother->Vx()*1e13;
695 // fpy = mother->Vy()*1e13;
696 // fpz = mother->Vz()*1e13;
697 // fpt = mother->T() *1e13*3e10;
702 if (fRotateToEventPlane) {
703 double tPhi = TMath::ATan2(fpy, fpx);
704 double tRad = TMath::Hypot(fpx, fpy);
706 fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
707 fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
710 tInfo->SetPDGPid(tPart->GetPdgCode());
712 if (fRotateToEventPlane) {
713 double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
714 double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
716 tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
717 tRad*TMath::Sin(tPhi - tReactionPlane),
721 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
722 Double_t mass2 = (tPart->Energy() *tPart->Energy() -
723 tPart->Px()*tPart->Px() -
724 tPart->Py()*tPart->Py() -
725 tPart->Pz()*tPart->Pz());
727 tInfo->SetMass(TMath::Sqrt(mass2));
731 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
732 trackCopy->SetHiddenInfo(tInfo);
735 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
737 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
742 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
744 tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
746 trackCopy->SetHiddenInfo(tInfo);
748 // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl;
750 //decision if we want this track
751 //if we using diffrent labels we want that this label was use for first time
752 //if we use hidden info we want to have match between sim data and ESD
753 if (tGoodMomentum==true)
755 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
756 realnofTracks++;//real number of tracks
768 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
770 AliCentrality *cent = fEvent->GetCentrality();
772 hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
773 // hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
774 hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
775 // hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
777 if (Debug()>1) printf(" FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
780 if (fEstEventMult == kGlobalCount)
781 hbtEvent->SetNormalizedMult(tNormMult);
782 else if (fEstEventMult == kTracklet)
783 hbtEvent->SetNormalizedMult(tTracklet);
784 else if (fEstEventMult == kITSTPC)
785 hbtEvent->SetNormalizedMult(tITSTPC);
786 else if (fEstEventMult == kITSPure)
787 hbtEvent->SetNormalizedMult(tITSPure);
788 else if (fEstEventMult == kSPDLayer1)
789 hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
790 else if (fEstEventMult == kV0Centrality) {
791 // centrality between 0 (central) and 1 (very peripheral)
794 if (cent->GetCentralityPercentile("V0M") < 0.00001)
795 hbtEvent->SetNormalizedMult(-1);
797 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
798 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
799 10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
803 if (tNormMultPos > tNormMultNeg)
804 hbtEvent->SetZDCParticipants(tNormMultPos);
806 hbtEvent->SetZDCParticipants(tNormMultNeg);
808 AliEventplane* ep = fEvent->GetEventplane();
811 hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
816 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
819 //___________________
820 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
822 // The chain loads the ESD for us
823 // You must provide the address where it can be found
826 //___________________
827 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
829 // The chain loads the stack for us
830 // You must provide the address where it can be found
833 //___________________
834 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
836 // The chain loads the generator event header for us
837 // You must provide the address where it can be found
838 fGenHeader = aGenHeader;
840 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
842 fRotateToEventPlane=dorotate;
844 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
846 // Calculates the number of sigma to the vertex.
858 bRes[0] = TMath::Sqrt(bCov[0]);
859 bRes[1] = TMath::Sqrt(bCov[2]);
861 // -----------------------------------
862 // How to get to a n-sigma cut?
864 // The accumulated statistics from 0 to d is
866 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
867 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
869 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
870 // Can this be expressed in a different way?
872 if (bRes[0] == 0 || bRes[1] ==0)
875 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
877 // stupid rounding problem screws up everything:
878 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
879 if (TMath::Exp(-d * d / 2) < 1e-10)
882 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
886 void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)