1 ////////////////////////////////////////////////////////////////////////////////
3 // AliFemtoEventReaderAOD - the reader class for the Alice AOD //
4 // Reads in AOD information and converts it into internal AliFemtoEvent //
5 // Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl //
6 // Adam Kisiel kisiel@mps.ohio-state.edu //
8 ////////////////////////////////////////////////////////////////////////////////
10 #include "AliFemtoEventReaderAOD.h"
14 #include "AliAODEvent.h"
15 #include "AliAODTrack.h"
16 #include "AliAODVertex.h"
17 #include "AliAODMCHeader.h"
18 #include "AliESDtrack.h"
20 #include "AliFmPhysicalHelixD.h"
21 #include "AliFmThreeVectorF.h"
23 #include "SystemOfUnits.h"
25 #include "AliFemtoEvent.h"
26 #include "AliFemtoModelHiddenInfo.h"
27 #include "AliFemtoModelGlobalHiddenInfo.h"
30 #include "AliAODpidUtil.h"
33 ClassImp(AliFemtoEventReaderAOD)
35 #if !(ST_NO_NAMESPACES)
36 using namespace units;
43 //____________________________
44 //constructor with 0 parameters , look at default settings
45 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
52 // fPWG2AODTracks(0x0),
64 // default constructor
65 fAllTrue.ResetAllBits(kTRUE);
66 fAllFalse.ResetAllBits(kFALSE);
71 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
72 AliFemtoEventReader(),
79 // fPWG2AODTracks(0x0),
92 fReadMC = aReader.fReadMC;
93 fReadV0 = aReader.fReadV0;
94 fInputFile = aReader.fInputFile;
95 fFileName = aReader.fFileName;
96 fNumberofEvent = aReader.fNumberofEvent;
97 fCurEvent = aReader.fCurEvent;
98 fEvent = new AliAODEvent();
99 fAodFile = new TFile(aReader.fAodFile->GetName());
100 fAllTrue.ResetAllBits(kTRUE);
101 fAllFalse.ResetAllBits(kFALSE);
102 fFilterBit = aReader.fFilterBit;
103 // fPWG2AODTracks = aReader.fPWG2AODTracks;
104 fAODpidUtil = aReader.fAODpidUtil;
105 fCentRange[0] = aReader.fCentRange[0];
106 fCentRange[1] = aReader.fCentRange[1];
107 fNoCentrality = aReader.fNoCentrality;
108 fUsePreCent = aReader.fUsePreCent;
112 AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
118 // if (fPWG2AODTracks) {
119 // fPWG2AODTracks->Delete();
120 // delete fPWG2AODTracks;
125 AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
127 // assignment operator
128 if (this == &aReader)
131 fInputFile = aReader.fInputFile;
132 fFileName = aReader.fFileName;
133 fNumberofEvent = aReader.fNumberofEvent;
134 fCurEvent = aReader.fCurEvent;
135 if (fTree) delete fTree;
136 if (fEvent) delete fEvent;
137 fEvent = new AliAODEvent();
138 if (fAodFile) delete fAodFile;
139 fAodFile = new TFile(aReader.fAodFile->GetName());
140 fAllTrue.ResetAllBits(kTRUE);
141 fAllFalse.ResetAllBits(kFALSE);
142 fFilterBit = aReader.fFilterBit;
143 // fPWG2AODTracks = aReader.fPWG2AODTracks;
144 fAODpidUtil = aReader.fAODpidUtil;
145 fCentRange[0] = aReader.fCentRange[0];
146 fCentRange[1] = aReader.fCentRange[1];
147 fUsePreCent = aReader.fUsePreCent;
148 fNoCentrality = aReader.fNoCentrality;
153 AliFemtoString AliFemtoEventReaderAOD::Report()
155 // create reader report
156 AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
161 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
163 //setting the name of file where names of AOD file are written
164 //it takes only this files which have good trees
166 fInputFile=string(inputFile);
167 ifstream infile(inputFile);
169 fTree = new TChain("aodTree");
171 if(infile.good()==true)
173 //checking if all give files have good tree inside
174 while (infile.eof()==false)
176 infile.getline(buffer,256);
177 TFile *aodFile=TFile::Open(buffer,"READ");
180 TTree* tree = (TTree*) aodFile->Get("aodTree");
183 // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
184 fTree->AddFile(buffer);
194 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
196 // read in a next hbt event from the chain
197 // convert it to AliFemtoEvent and return
198 // for further analysis
199 AliFemtoEvent *hbtEvent = 0;
200 cout<<"reader"<<endl;
201 if (fCurEvent==fNumberofEvent)//open next file
203 if(fNumberofEvent==0)
205 fEvent=new AliAODEvent();
206 fEvent->ReadFromTree(fTree);
208 // Check for the existence of the additional information
209 // fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
211 // if (fPWG2AODTracks) {
212 // cout << "Found additional PWG2 specific information in the AOD!" << endl;
213 // cout << "Reading only tracks with the additional information" << endl;
216 fNumberofEvent=fTree->GetEntries();
217 // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
220 else //no more data to read
222 cout<<"no more files "<<hbtEvent<<endl;
228 cout<<"starting to read event "<<fCurEvent<<endl;
229 fTree->GetEvent(fCurEvent);//getting next event
230 // cout << "Read event " << fEvent << " from file " << fTree << endl;
232 hbtEvent = new AliFemtoEvent;
234 CopyAODtoFemtoEvent(hbtEvent);
241 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
244 // A function that reads in the AOD event
245 // and transfers the neccessary information into
246 // the internal AliFemtoEvent
248 // setting global event characteristics
249 tEvent->SetRunNumber(fEvent->GetRunNumber());
250 tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
251 tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
252 tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
253 tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
254 tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
255 tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
256 tEvent->SetZDCParticipants(0);
257 tEvent->SetTriggerMask(fEvent->GetTriggerMask());
258 tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
260 // Attempt to access MC header
264 mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
266 cout << "AOD MC information requested, but no header found!" << endl;
269 mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
271 cout << "AOD MC information requested, but no particle array found!" << endl;
275 tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
279 motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
280 for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
282 // Read in mother ids
283 AliAODMCParticle *motherpart;
284 for (int ip=0; ip<mcP->GetEntries(); ip++) {
285 motherpart = (AliAODMCParticle *) mcP->At(ip);
286 if (motherpart->GetDaughter(0) > 0)
287 motherids[motherpart->GetDaughter(0)] = ip;
288 if (motherpart->GetDaughter(1) > 0)
289 motherids[motherpart->GetDaughter(1)] = ip;
293 // Primary Vertex position
295 fEvent->GetPrimaryVertex()->GetPosition(fV1);
297 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
298 tEvent->SetPrimVertPos(vertex);
300 //starting to reading tracks
301 int nofTracks=0; //number of reconstructed tracks in event
303 // Check to see whether the additional info exists
304 // if (fPWG2AODTracks)
305 // nofTracks=fPWG2AODTracks->GetEntries();
307 nofTracks=fEvent->GetNumberOfTracks();
309 AliCentrality *cent = fEvent->GetCentrality();
311 if (!fNoCentrality && cent && fUsePreCent) {
312 if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
313 (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
315 cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
321 int realnofTracks=0; // number of track which we use in a analysis
325 for (int il=0; il<20000; il++) labels[il] = -1;
327 // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
328 for (int i=0;i<nofTracks;i++) {
329 const AliAODTrack *aodtrack=fEvent->GetTrack(i);
330 if (!aodtrack->TestFilterBit(fFilterBit)) {
331 if(aodtrack->GetID() < 0) continue;
332 labels[aodtrack->GetID()] = i;
337 for (int i=0;i<nofTracks;i++)
339 AliFemtoTrack* trackCopy = new AliFemtoTrack();
341 // if (fPWG2AODTracks) {
342 // // Read tracks from the additional pwg2 specific AOD part
344 // // Note that in that case all the AOD tracks without the
345 // // additional information will be ignored !
346 // AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
348 // // Getting the AOD track through the ref of the additional info
349 // AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
350 // if (!aodtrack->TestFilterBit(fFilterBit)) {
356 // if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
357 // if (aodtrack->Chi2perNDF() < 6.0)
358 // if (aodtrack->Eta() < 0.9)
362 // CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
365 // // Fill the hidden information with the simulated data
366 // // Int_t pLabel = aodtrack->GetLabel();
367 // AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
369 // // Check the mother information
371 // // Using the new way of storing the freeze-out information
372 // // Final state particle is stored twice on the stack
373 // // one copy (mother) is stored with original freeze-out information
374 // // and is not tracked
375 // // the other one (daughter) is stored with primary vertex position
378 // // Freeze-out coordinates
379 // double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
380 // fpx = tPart->Xv() - fV1[0];
381 // fpy = tPart->Yv() - fV1[1];
382 // fpz = tPart->Zv() - fV1[2];
385 // AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
386 // tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
393 // // cout << "Looking for mother ids " << endl;
394 // if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
395 // // cout << "Got mother id" << endl;
396 // AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
397 // // Check if this is the same particle stored twice on the stack
398 // if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
399 // // It is the same particle
400 // // Read in the original freeze-out information
401 // // and convert it from to [fm]
404 // // fpx = mother->Xv()*1e13*0.197327;
405 // // fpy = mother->Yv()*1e13*0.197327;
406 // // fpz = mother->Zv()*1e13*0.197327;
407 // // fpt = mother->T() *1e13*0.197327*0.5;
410 // // Therminator style
411 // fpx = mother->Xv()*1e13;
412 // fpy = mother->Yv()*1e13;
413 // fpz = mother->Zv()*1e13;
414 // fpt = mother->T() *1e13*3e10;
419 // // if (fRotateToEventPlane) {
420 // // double tPhi = TMath::ATan2(fpy, fpx);
421 // // double tRad = TMath::Hypot(fpx, fpy);
423 // // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
424 // // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
427 // tInfo->SetPDGPid(tPart->GetPdgCode());
429 // // if (fRotateToEventPlane) {
430 // // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
431 // // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
433 // // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
434 // // tRad*TMath::Sin(tPhi - tReactionPlane),
438 // tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
439 // Double_t mass2 = (tPart->E() *tPart->E() -
440 // tPart->Px()*tPart->Px() -
441 // tPart->Py()*tPart->Py() -
442 // tPart->Pz()*tPart->Pz());
444 // tInfo->SetMass(TMath::Sqrt(mass2));
446 // tInfo->SetMass(0.0);
448 // tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
449 // trackCopy->SetHiddenInfo(tInfo);
454 // aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
455 // const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
456 // // Check the sanity of the tracks - reject zero momentum tracks
457 // if (ktP.Mag() == 0) {
463 // No additional information exists
464 // Read in the normal AliAODTracks
466 // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
467 AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
469 if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
471 if (!aodtrack->TestFilterBit(fFilterBit)) {
476 //counting particles to set multiplicity
479 if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
480 if(impact[0]<2.4 && TMath::Abs(impact[1]+fV1[2])<3.2)
481 if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
482 if (aodtrack->Chi2perNDF() < 4.0)
483 if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
484 if (aodtrack->GetTPCNcls() > 50)
485 if (aodtrack->Eta() < 0.8)
490 CopyAODtoFemtoTrack(aodtrack, trackCopy);
492 // copying PID information from the correspondent track
493 // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
496 AliAODTrack *aodtrackpid;
497 if(fFilterBit == (1 << (7))) //for TPC Only tracks we have to copy PID information from corresponding global tracks
498 aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
500 aodtrackpid = fEvent->GetTrack(i);
501 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
504 // Fill the hidden information with the simulated data
505 // Int_t pLabel = aodtrack->GetLabel();
506 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
508 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
509 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
514 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
516 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
517 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
521 // Check the mother information
523 // Using the new way of storing the freeze-out information
524 // Final state particle is stored twice on the stack
525 // one copy (mother) is stored with original freeze-out information
526 // and is not tracked
527 // the other one (daughter) is stored with primary vertex position
530 // Freeze-out coordinates
531 fpx = tPart->Xv() - fV1[0];
532 fpy = tPart->Yv() - fV1[1];
533 fpz = tPart->Zv() - fV1[2];
536 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
543 // cout << "Looking for mother ids " << endl;
544 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
545 // cout << "Got mother id" << endl;
546 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
547 // Check if this is the same particle stored twice on the stack
549 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
550 // It is the same particle
551 // Read in the original freeze-out information
552 // and convert it from to [fm]
555 // fpx = mother->Xv()*1e13*0.197327;
556 // fpy = mother->Yv()*1e13*0.197327;
557 // fpz = mother->Zv()*1e13*0.197327;
558 // fpt = mother->T() *1e13*0.197327*0.5;
562 fpx = mother->Xv()*1e13;
563 fpy = mother->Yv()*1e13;
564 fpz = mother->Zv()*1e13;
565 // fpt = mother->T() *1e13*3e10;
571 // if (fRotateToEventPlane) {
572 // double tPhi = TMath::ATan2(fpy, fpx);
573 // double tRad = TMath::Hypot(fpx, fpy);
575 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
576 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
579 tInfo->SetPDGPid(tPart->GetPdgCode());
581 // if (fRotateToEventPlane) {
582 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
583 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
585 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
586 // tRad*TMath::Sin(tPhi - tReactionPlane),
590 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
591 Double_t mass2 = (tPart->E() *tPart->E() -
592 tPart->Px()*tPart->Px() -
593 tPart->Py()*tPart->Py() -
594 tPart->Pz()*tPart->Pz());
596 tInfo->SetMass(TMath::Sqrt(mass2));
600 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
602 trackCopy->SetHiddenInfo(tInfo);
606 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
607 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
608 // Check the sanity of the tracks - reject zero momentum tracks
609 if (ktP.Mag() == 0) {
616 tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
617 realnofTracks++;//real number of tracks
620 tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
621 tEvent->SetNormalizedMult(tracksPrim);
623 if (!fNoCentrality) {
624 // AliCentrality *cent = fEvent->GetCentrality();
625 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
626 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
629 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
630 // tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
631 tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
632 // tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
636 tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
639 if (mcP) delete [] motherids;
641 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
646 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
647 AliAODv0* aodv0 = fEvent->GetV0(i);
648 if (!aodv0) continue;
649 if(aodv0->GetNDaughters()>2) continue;
650 if(aodv0->GetNProngs()>2) continue;
651 if(aodv0->GetCharge()!=0) continue;
652 if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
653 if(aodv0->CosPointingAngle(fV1)<0.998) continue;
654 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
656 CopyAODtoFemtoV0(aodv0, trackCopyV0);
657 tEvent->V0Collection()->push_back(trackCopyV0);
658 //cout<<"Pushback v0 to v0collection"<<endl;
664 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
665 AliFemtoTrack *tFemtoTrack
666 // AliPWG2AODTrack *tPWG2AODTrack
669 // Copy the track information from the AOD into the internal AliFemtoTrack
670 // If it exists, use the additional information from the PWG2 AOD
672 // Primary Vertex position
674 fEvent->GetPrimaryVertex()->GetPosition(fV1);
675 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
677 tFemtoTrack->SetCharge(tAodTrack->Charge());
680 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
681 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
682 tFemtoTrack->SetP(v);//setting momentum
683 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
684 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
685 //setting track helix
686 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
687 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
688 tFemtoTrack->SetHelix(helix);
691 tFemtoTrack->SetTrackId(tAodTrack->GetID());
692 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
693 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
695 // Track quality information
697 tAodTrack->GetCovMatrix(covmat);
699 // ! DCA information is done in CopyPIDtoFemtoTrack()
702 // double covimpact[3];
704 // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
705 // //cout << "sth went wrong with dca propagation" << endl;
706 // tFemtoTrack->SetImpactD(-1000.0);
707 // tFemtoTrack->SetImpactZ(-1000.0);
711 // tFemtoTrack->SetImpactD(impact[0]);
712 // tFemtoTrack->SetImpactZ(impact[1]+fV1[2]);
715 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
716 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
718 // tFemtoTrack->SetImpactD(0.0);
719 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
721 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
724 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
725 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
729 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
730 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
731 // << tAodTrack->Yv() - fV1[1]
732 // << "xv = " << tAodTrack->Xv() << endl
733 // << "fv1[0] = " << fV1[0] << endl
734 // << "yv = " << tAodTrack->Yv() << endl
735 // << "fv1[1] = " << fV1[1] << endl
736 // << "zv = " << tAodTrack->Zv() << endl
737 // << "fv1[2] = " << fV1[2] << endl
738 // << "impact[0] = " << impact[0] << endl
739 // << "impact[1] = " << impact[1] << endl
742 tFemtoTrack->SetCdd(covmat[0]);
743 tFemtoTrack->SetCdz(covmat[1]);
744 tFemtoTrack->SetCzz(covmat[2]);
745 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
746 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
747 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
748 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
749 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
750 tFemtoTrack->SetTPCsignalN(1);
751 tFemtoTrack->SetTPCsignalS(1);
752 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
754 // if (tPWG2AODTrack) {
755 // // Copy the PWG2 specific information if it exists
756 // tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
757 // tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
759 // double xtpc[3] = {0,0,0};
760 // tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
761 // tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
762 // tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
763 // tFemtoTrack->SetNominalTPCExitPoint(xtpc);
766 // If not use dummy values
767 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
768 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
771 float globalPositionsAtRadii[9][3];
772 float bfield = 5*fMagFieldSign;
773 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
774 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
775 double **tpcPositions;
776 tpcPositions = new double*[9];
778 tpcPositions[i] = new double[3];
779 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
782 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
783 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
784 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
786 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
787 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
788 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
792 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
796 for (int ik=0; ik<3; ik++) {
799 tFemtoTrack->SetKinkIndexes(indexes);
802 void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
804 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
805 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
806 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
807 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
808 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
809 tFemtoV0->SetdecayVertexV0(decayvertex);
810 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
811 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
812 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
813 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
814 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
815 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
816 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
817 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
818 tFemtoV0->SetmomPos(mompos);
819 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
820 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
821 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
822 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
823 tFemtoV0->SetmomNeg(momneg);
825 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
826 //void SettpcHitsPos(const int& i);
827 //void SettpcHitsNeg(const int& i);
829 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
830 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
832 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
833 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
834 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
835 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
836 tFemtoV0->SetmomV0(momv0);
837 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
838 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
839 tFemtoV0->SeteLambda(tAODv0->ELambda());
840 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
841 tFemtoV0->SetePosProton(tAODv0->EPosProton());
842 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
843 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
844 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
845 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
846 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
847 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
849 //void SetcTauLambda( float x);
850 //void SetcTauK0Short( float x);
852 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
853 tFemtoV0->SetptV0(tAODv0->Pt());
854 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
855 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
856 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
857 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
858 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
860 tFemtoV0->SetidNeg(tAODv0->GetNegID());
861 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
862 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
863 tFemtoV0->SetidPos(tAODv0->GetPosID());
865 tFemtoV0->SetEtaV0(tAODv0->Eta());
866 tFemtoV0->SetPhiV0(tAODv0->Phi());
867 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
868 //tFemtoV0->SetYV0(tAODv0->Y());
870 //void SetdedxNeg(float x);
871 //void SeterrdedxNeg(float x);//Gael 04Fev2002
872 //void SetlendedxNeg(float x);//Gael 04Fev2002
873 //void SetdedxPos(float x);
874 //void SeterrdedxPos(float x);//Gael 04Fev2002
875 //void SetlendedxPos(float x);//Gael 04Fev2002
877 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
878 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
880 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
881 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
883 if(trackpos && trackneg)
885 tFemtoV0->SetEtaPos(trackpos->Eta());
886 tFemtoV0->SetEtaNeg(trackneg->Eta());
887 tFemtoV0->SetptotPos(tAODv0->PProng(0));
888 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
889 tFemtoV0->SetptPos(trackpos->Pt());
890 tFemtoV0->SetptNeg(trackneg->Pt());
892 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
893 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
894 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
895 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
896 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
897 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
898 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
899 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
900 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
901 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
902 tFemtoV0->SetStatusPos(trackpos->GetStatus());
903 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
905 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
906 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
907 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
908 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
909 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
910 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
913 float bfield = 5*fMagFieldSign;
914 float globalPositionsAtRadiiPos[9][3];
915 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
916 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
917 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
919 float globalPositionsAtRadiiNeg[9][3];
920 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
921 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
922 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
924 AliFemtoThreeVector tmpVec;
925 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
926 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
928 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
929 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
931 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
932 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
934 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
935 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
937 AliFemtoThreeVector vecTpcPos[9];
938 AliFemtoThreeVector vecTpcNeg[9];
941 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
942 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
944 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
945 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
947 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
948 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
950 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
951 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
953 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFmismatch)>0)
955 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFmismatch)>0)
957 tFemtoV0->SetPosNSigmaTOFK(-1000);
958 tFemtoV0->SetNegNSigmaTOFK(-1000);
959 tFemtoV0->SetPosNSigmaTOFP(-1000);
960 tFemtoV0->SetNegNSigmaTOFP(-1000);
961 tFemtoV0->SetPosNSigmaTOFPi(-1000);
962 tFemtoV0->SetNegNSigmaTOFPi(-1000);
964 tFemtoV0->SetTOFProtonTimePos(-1000);
965 tFemtoV0->SetTOFPionTimePos(-1000);
966 tFemtoV0->SetTOFKaonTimePos(-1000);
967 tFemtoV0->SetTOFProtonTimeNeg(-1000);
968 tFemtoV0->SetTOFPionTimeNeg(-1000);
969 tFemtoV0->SetTOFKaonTimeNeg(-1000);
974 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
975 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
976 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
977 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
978 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
979 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
981 double TOFSignalPos = trackpos->GetTOFsignal();
982 double TOFSignalNeg = trackneg->GetTOFsignal();
985 trackpos->GetIntegratedTimes(pidPos);
986 trackneg->GetIntegratedTimes(pidNeg);
988 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
989 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
990 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
991 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
992 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
993 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
998 tFemtoV0->SetStatusPos(999);
999 tFemtoV0->SetStatusNeg(999);
1001 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1004 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1006 fFilterBit = (1 << (ibit));
1009 void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1015 void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1020 AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1022 if (aLabel < 0) return 0;
1023 AliAODMCParticle *aodP;
1025 if (aLabel > mcP->GetEntries())
1026 posstack = mcP->GetEntries();
1030 aodP = (AliAODMCParticle *) mcP->At(posstack);
1031 if (aodP->GetLabel() > posstack) {
1033 aodP = (AliAODMCParticle *) mcP->At(posstack);
1034 if (aodP->GetLabel() == aLabel) return aodP;
1037 while (posstack > 0);
1041 aodP = (AliAODMCParticle *) mcP->At(posstack);
1042 if (aodP->GetLabel() == aLabel) return aodP;
1045 while (posstack < mcP->GetEntries());
1051 void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1052 AliFemtoTrack *tFemtoTrack)
1055 // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
1058 double covimpact[3];
1060 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1061 //cout << "sth went wrong with dca propagation" << endl;
1062 tFemtoTrack->SetImpactD(-1000.0);
1063 tFemtoTrack->SetImpactZ(-1000.0);
1067 tFemtoTrack->SetImpactD(impact[0]);
1068 tFemtoTrack->SetImpactZ(impact[1]);
1072 tAodTrack->GetPID(aodpid);
1073 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1074 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1075 tFemtoTrack->SetPidProbPion(aodpid[2]);
1076 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1077 tFemtoTrack->SetPidProbProton(aodpid[4]);
1079 aodpid[0] = -100000.0;
1080 aodpid[1] = -100000.0;
1081 aodpid[2] = -100000.0;
1082 aodpid[3] = -100000.0;
1083 aodpid[4] = -100000.0;
1087 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1088 tTOF = tAodTrack->GetTOFsignal();
1089 tAodTrack->GetIntegratedTimes(aodpid);
1092 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
1094 ////// TPC ////////////////////////////////////////////
1096 float nsigmaTPCK=-1000.;
1097 float nsigmaTPCPi=-1000.;
1098 float nsigmaTPCP=-1000.;
1100 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1102 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1103 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1104 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1105 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1108 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1109 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1110 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
1112 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1113 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1114 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1116 tFemtoTrack->SetTPCsignalN(1);
1117 tFemtoTrack->SetTPCsignalS(1);
1118 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
1120 ///////TOF//////////////////////
1123 float nsigmaTOFPi=-1000.;
1124 float nsigmaTOFK=-1000.;
1125 float nsigmaTOFP=-1000.;
1127 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1128 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
1129 (tAodTrack->GetStatus() & AliESDtrack::kTIME) && //AliESDtrack::kTIME=0x80000000
1130 !(tAodTrack->GetStatus() & AliESDtrack::kTOFmismatch)) //AliESDtrack::kTOFmismatch=0x100000
1132 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
1135 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1136 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1137 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1139 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1140 Double_t tof=tAodTrack->GetTOFsignal();
1141 if(tof > 0.) vp=len/tof/0.03;
1144 tFemtoTrack->SetVTOF(vp);
1145 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1146 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1147 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1150 //////////////////////////////////////
1154 void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1156 fCentRange[0] = min; fCentRange[1] = max;
1161 void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1163 fNoCentrality = anocent;
1164 if (fNoCentrality) fUsePreCent = 0;
1167 void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1169 fAODpidUtil = aAODpidUtil;
1170 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1174 void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1185 void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1187 // Gets the global position of the track at nine different radii in the TPC
1188 // track is the track you want to propagate
1189 // bfield is the magnetic field of your event
1190 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1192 // Initialize the array to something indicating there was no propagation
1193 for(Int_t i=0;i<9;i++){
1194 for(Int_t j=0;j<3;j++){
1195 globalPositionsAtRadii[i][j]=-9999.;
1199 // Make a copy of the track to not change parameters of the track
1200 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1201 //printf("\nAfter CopyFromVTrack\n");
1204 // The global position of the the track
1205 Double_t xyz[3]={-9999.,-9999.,-9999.};
1207 // Counter for which radius we want
1209 // The radii at which we get the global positions
1210 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1211 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1212 // The global radius we are at
1213 Float_t globalRadius=0;
1215 // Propagation is done in local x of the track
1216 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1217 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1218 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1219 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1220 // adds to the global radius
1222 // Stop if the propagation was not succesful. This can happen for low pt tracks
1223 // that don't reach outer radii
1224 if(!etp.PropagateTo(x,bfield))break;
1225 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1226 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1228 // Roughly reached the radius we want
1229 if(globalRadius > Rwanted[iR]){
1231 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1232 while (globalRadius>Rwanted[iR]){
1234 // printf("propagating to x %5.2f\n",x);
1235 if(!etp.PropagateTo(x,bfield))break;
1236 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1237 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1239 //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]);
1240 globalPositionsAtRadii[iR][0]=xyz[0];
1241 globalPositionsAtRadii[iR][1]=xyz[1];
1242 globalPositionsAtRadii[iR][2]=xyz[2];
1243 // Indicate we want the next radius