]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix the bug of setting extrance and exit points of V0 daughter tracks
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderAOD.cxx
CommitLineData
76ce4b5b 1////////////////////////////////////////////////////////////////////////////////
2// //
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 //
7// //
8////////////////////////////////////////////////////////////////////////////////
9
10#include "AliFemtoEventReaderAOD.h"
11
12#include "TFile.h"
13#include "TTree.h"
14#include "AliAODEvent.h"
15#include "AliAODTrack.h"
16#include "AliAODVertex.h"
17#include "AliAODMCHeader.h"
18#include "AliESDtrack.h"
19
20#include "AliFmPhysicalHelixD.h"
21#include "AliFmThreeVectorF.h"
22
23#include "SystemOfUnits.h"
24
25#include "AliFemtoEvent.h"
26#include "AliFemtoModelHiddenInfo.h"
27#include "AliFemtoModelGlobalHiddenInfo.h"
28#include "AliPID.h"
29
30#include "AliAODpidUtil.h"
31
32ClassImp(AliFemtoEventReaderAOD)
33
34#if !(ST_NO_NAMESPACES)
87ede832 35using namespace units;
76ce4b5b 36#endif
37
38using namespace std;
973a91f8 39
40double fV1[3];
41
76ce4b5b 42//____________________________
87ede832 43//constructor with 0 parameters , look at default settings
76ce4b5b 44AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
45 fNumberofEvent(0),
46 fCurEvent(0),
47 fEvent(0x0),
48 fAllTrue(160),
49 fAllFalse(160),
50 fFilterBit(0),
64536eaf 51 fFilterMask(0),
76ce4b5b 52 // fPWG2AODTracks(0x0),
53 fReadMC(0),
973a91f8 54 fReadV0(0),
76ce4b5b 55 fUsePreCent(0),
1777446b 56 fEstEventMult(kCentrality),
76ce4b5b 57 fAODpidUtil(0),
1777446b 58 fAODheader(0),
98e67b33 59 fInputFile(""),
76ce4b5b 60 fTree(0x0),
ba3c23a4 61 fAodFile(0x0),
734c121a 62 fMagFieldSign(1),
63 fisEPVZ(kTRUE),
50bb1be0 64 fpA2013(kFALSE),
65 fDCAglobalTrack(kFALSE)
76ce4b5b 66{
67 // default constructor
68 fAllTrue.ResetAllBits(kTRUE);
69 fAllFalse.ResetAllBits(kFALSE);
70 fCentRange[0] = 0;
71 fCentRange[1] = 1000;
72}
73
74AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
75 AliFemtoEventReader(),
76 fNumberofEvent(0),
77 fCurEvent(0),
78 fEvent(0x0),
79 fAllTrue(160),
80 fAllFalse(160),
81 fFilterBit(0),
64536eaf 82 fFilterMask(0),
76ce4b5b 83 // fPWG2AODTracks(0x0),
84 fReadMC(0),
973a91f8 85 fReadV0(0),
76ce4b5b 86 fUsePreCent(0),
1777446b 87 fEstEventMult(kCentrality),
76ce4b5b 88 fAODpidUtil(0),
1777446b 89 fAODheader(0),
98e67b33 90 fInputFile(""),
76ce4b5b 91 fTree(0x0),
ba3c23a4 92 fAodFile(0x0),
734c121a 93 fMagFieldSign(1),
94 fisEPVZ(kTRUE),
50bb1be0 95 fpA2013(kFALSE),
96 fDCAglobalTrack(kFALSE)
76ce4b5b 97{
98 // copy constructor
973a91f8 99 fReadMC = aReader.fReadMC;
100 fReadV0 = aReader.fReadV0;
76ce4b5b 101 fInputFile = aReader.fInputFile;
76ce4b5b 102 fNumberofEvent = aReader.fNumberofEvent;
103 fCurEvent = aReader.fCurEvent;
104 fEvent = new AliAODEvent();
105 fAodFile = new TFile(aReader.fAodFile->GetName());
106 fAllTrue.ResetAllBits(kTRUE);
107 fAllFalse.ResetAllBits(kFALSE);
108 fFilterBit = aReader.fFilterBit;
109 // fPWG2AODTracks = aReader.fPWG2AODTracks;
110 fAODpidUtil = aReader.fAODpidUtil;
1777446b 111 fAODheader = aReader.fAODheader;
76ce4b5b 112 fCentRange[0] = aReader.fCentRange[0];
113 fCentRange[1] = aReader.fCentRange[1];
1777446b 114 fEstEventMult = aReader.fEstEventMult;
973a91f8 115 fUsePreCent = aReader.fUsePreCent;
734c121a 116 fpA2013 = aReader.fpA2013;
50bb1be0 117 fDCAglobalTrack = aReader.fDCAglobalTrack;
118
76ce4b5b 119}
120//__________________
121//Destructor
122AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
123{
124 // destructor
125 delete fTree;
126 delete fEvent;
127 delete fAodFile;
128// if (fPWG2AODTracks) {
129// fPWG2AODTracks->Delete();
130// delete fPWG2AODTracks;
131// }
132}
133
134//__________________
135AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
136{
137 // assignment operator
138 if (this == &aReader)
87ede832 139 return *this;
76ce4b5b 140
141 fInputFile = aReader.fInputFile;
76ce4b5b 142 fNumberofEvent = aReader.fNumberofEvent;
143 fCurEvent = aReader.fCurEvent;
144 if (fTree) delete fTree;
145 if (fEvent) delete fEvent;
146 fEvent = new AliAODEvent();
147 if (fAodFile) delete fAodFile;
148 fAodFile = new TFile(aReader.fAodFile->GetName());
149 fAllTrue.ResetAllBits(kTRUE);
150 fAllFalse.ResetAllBits(kFALSE);
151 fFilterBit = aReader.fFilterBit;
64536eaf 152 fFilterMask = aReader.fFilterMask;
76ce4b5b 153 // fPWG2AODTracks = aReader.fPWG2AODTracks;
154 fAODpidUtil = aReader.fAODpidUtil;
1777446b 155 fAODheader = aReader.fAODheader;
76ce4b5b 156 fCentRange[0] = aReader.fCentRange[0];
157 fCentRange[1] = aReader.fCentRange[1];
973a91f8 158 fUsePreCent = aReader.fUsePreCent;
1777446b 159 fEstEventMult = aReader.fEstEventMult;
734c121a 160 fpA2013 = aReader.fpA2013;
e81ec993 161 fDCAglobalTrack = aReader.fDCAglobalTrack;
76ce4b5b 162
163 return *this;
164}
165//__________________
166AliFemtoString AliFemtoEventReaderAOD::Report()
167{
168 // create reader report
169 AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
170 return temp;
171}
172
173//__________________
174void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
175{
87ede832 176 //setting the name of file where names of AOD file are written
76ce4b5b 177 //it takes only this files which have good trees
178 char buffer[256];
179 fInputFile=string(inputFile);
180 ifstream infile(inputFile);
181
182 fTree = new TChain("aodTree");
183
184 if(infile.good()==true)
87ede832 185 {
186 //checking if all give files have good tree inside
187 while (infile.eof()==false)
188 {
189 infile.getline(buffer,256);
190 TFile *aodFile=TFile::Open(buffer,"READ");
191 if (aodFile!=0x0)
192 {
76ce4b5b 193 TTree* tree = (TTree*) aodFile->Get("aodTree");
194 if (tree!=0x0)
87ede832 195 {
196 // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
197 fTree->AddFile(buffer);
198 delete tree;
199 }
200 aodFile->Close();
76ce4b5b 201 }
87ede832 202 delete aodFile;
76ce4b5b 203 }
87ede832 204 }
76ce4b5b 205}
206
207AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
208{
209 // read in a next hbt event from the chain
210 // convert it to AliFemtoEvent and return
211 // for further analysis
212 AliFemtoEvent *hbtEvent = 0;
87ede832 213 // cout<<"reader"<<endl;
214 if (fCurEvent==fNumberofEvent)//open next file
215 {
216 if(fNumberofEvent==0)
76ce4b5b 217 {
87ede832 218 fEvent=new AliAODEvent();
219 fEvent->ReadFromTree(fTree);
76ce4b5b 220
87ede832 221 // Check for the existence of the additional information
76ce4b5b 222// fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
223
224// if (fPWG2AODTracks) {
225// cout << "Found additional PWG2 specific information in the AOD!" << endl;
226// cout << "Reading only tracks with the additional information" << endl;
227// }
228
87ede832 229 fNumberofEvent=fTree->GetEntries();
230 // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
231 fCurEvent=0;
232 }
233 else //no more data to read
234 {
235 // cout<<"no more files "<<hbtEvent<<endl;
236 fReaderStatus=1;
237 return hbtEvent;
238 }
239 }
76ce4b5b 240
87ede832 241 // cout<<"starting to read event "<<fCurEvent<<endl;
76ce4b5b 242 fTree->GetEvent(fCurEvent);//getting next event
243 // cout << "Read event " << fEvent << " from file " << fTree << endl;
87ede832 244
76ce4b5b 245 hbtEvent = new AliFemtoEvent;
246
247 CopyAODtoFemtoEvent(hbtEvent);
248 fCurEvent++;
249
250
87ede832 251 return hbtEvent;
76ce4b5b 252}
253
254void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
255{
256
257 // A function that reads in the AOD event
258 // and transfers the neccessary information into
259 // the internal AliFemtoEvent
260
261 // setting global event characteristics
262 tEvent->SetRunNumber(fEvent->GetRunNumber());
263 tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
264 tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
265 tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
266 tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
267 tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
268 tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
269 tEvent->SetZDCParticipants(0);
270 tEvent->SetTriggerMask(fEvent->GetTriggerMask());
271 tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
87ede832 272
76ce4b5b 273 // Attempt to access MC header
274 AliAODMCHeader *mcH;
275 TClonesArray *mcP=0;
276 if (fReadMC) {
277 mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
278 if (!mcH) {
279 cout << "AOD MC information requested, but no header found!" << endl;
280 }
281
282 mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
283 if (!mcP) {
284 cout << "AOD MC information requested, but no particle array found!" << endl;
285 }
286 }
287
288 tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
289
290 Int_t *motherids=0;
291 if (mcP) {
292 motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
293 for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
294
295 // Read in mother ids
296 AliAODMCParticle *motherpart;
297 for (int ip=0; ip<mcP->GetEntries(); ip++) {
298 motherpart = (AliAODMCParticle *) mcP->At(ip);
299 if (motherpart->GetDaughter(0) > 0)
87ede832 300 motherids[motherpart->GetDaughter(0)] = ip;
76ce4b5b 301 if (motherpart->GetDaughter(1) > 0)
87ede832 302 motherids[motherpart->GetDaughter(1)] = ip;
76ce4b5b 303 }
304 }
305
306 // Primary Vertex position
973a91f8 307 // double fV1[3];
734c121a 308 if(fpA2013)
87ede832 309 {
310 const AliAODVertex* trkVtx = fEvent->GetPrimaryVertex();
311 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
312 TString vtxTtl = trkVtx->GetTitle();
313 if (!vtxTtl.Contains("VertexerTracks")) return;
314 Float_t zvtx = trkVtx->GetZ();
315 const AliAODVertex* spdVtx = fEvent->GetPrimaryVertexSPD();
316 if (spdVtx->GetNContributors()<=0) return;
317 TString vtxTyp = spdVtx->GetTitle();
318 Double_t cov[6]={0};
319 spdVtx->GetCovarianceMatrix(cov);
320 Double_t zRes = TMath::Sqrt(cov[5]);
321 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
322 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
323
324 if (TMath::Abs(zvtx) > 10) return;
325 }
76ce4b5b 326 fEvent->GetPrimaryVertex()->GetPosition(fV1);
327
328 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
329 tEvent->SetPrimVertPos(vertex);
87ede832 330
76ce4b5b 331 //starting to reading tracks
332 int nofTracks=0; //number of reconstructed tracks in event
333
334 // Check to see whether the additional info exists
1777446b 335 // if (fPWG2AODTracks)
336 // nofTracks=fPWG2AODTracks->GetEntries();
337 // else
338 nofTracks=fEvent->GetNumberOfTracks();
87ede832 339
1777446b 340 AliEventplane *ep = fEvent->GetEventplane();
341 if (ep) {
342 tEvent->SetEP(ep);
87ede832 343 if (fisEPVZ)
344 tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
345 else
346 tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
1777446b 347 }
76ce4b5b 348
36cfef5c 349 AliCentrality *cent = fEvent->GetCentrality();
87ede832 350
1777446b 351 if (!fEstEventMult && cent && fUsePreCent) {
76ce4b5b 352 if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
87ede832 353 (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
354 {
355 // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
356
357 return;
358 }
76ce4b5b 359 }
360
361 int realnofTracks=0; // number of track which we use in a analysis
87ede832 362 int tracksPrim=0;
76ce4b5b 363
87ede832 364 int labels[20000];
76ce4b5b 365 for (int il=0; il<20000; il++) labels[il] = -1;
366
367 // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
368 for (int i=0;i<nofTracks;i++) {
369 const AliAODTrack *aodtrack=fEvent->GetTrack(i);
370 if (!aodtrack->TestFilterBit(fFilterBit)) {
973a91f8 371 if(aodtrack->GetID() < 0) continue;
76ce4b5b 372 labels[aodtrack->GetID()] = i;
373 }
374 }
375
ba3c23a4 376 int tNormMult = 0;
76ce4b5b 377 for (int i=0;i<nofTracks;i++)
87ede832 378 {
379 AliFemtoTrack* trackCopy = new AliFemtoTrack();
76ce4b5b 380
381// if (fPWG2AODTracks) {
382// // Read tracks from the additional pwg2 specific AOD part
383// // if they exist
87ede832 384// // Note that in that case all the AOD tracks without the
76ce4b5b 385// // additional information will be ignored !
386// AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
387
388// // Getting the AOD track through the ref of the additional info
87ede832 389// AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
76ce4b5b 390// if (!aodtrack->TestFilterBit(fFilterBit)) {
391// delete trackCopy;
392// continue;
393// }
394
87ede832 395
76ce4b5b 396// if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
87ede832 397// if (aodtrack->Chi2perNDF() < 6.0)
76ce4b5b 398// if (aodtrack->Eta() < 0.9)
399// tNormMult++;
400
401
402// CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
87ede832 403
76ce4b5b 404// if (mcP) {
405// // Fill the hidden information with the simulated data
406// // Int_t pLabel = aodtrack->GetLabel();
407// AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
408
409// // Check the mother information
87ede832 410
76ce4b5b 411// // Using the new way of storing the freeze-out information
412// // Final state particle is stored twice on the stack
413// // one copy (mother) is stored with original freeze-out information
414// // and is not tracked
415// // the other one (daughter) is stored with primary vertex position
416// // and is tracked
87ede832 417
76ce4b5b 418// // Freeze-out coordinates
419// double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
420// fpx = tPart->Xv() - fV1[0];
421// fpy = tPart->Yv() - fV1[1];
422// fpz = tPart->Zv() - fV1[2];
423// fpt = tPart->T();
424
425// AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
426// tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
427
428// fpx *= 1e13;
429// fpy *= 1e13;
430// fpz *= 1e13;
431// fpt *= 1e13;
87ede832 432
76ce4b5b 433// // cout << "Looking for mother ids " << endl;
434// if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
435// // cout << "Got mother id" << endl;
436// AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
437// // Check if this is the same particle stored twice on the stack
438// if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
439// // It is the same particle
440// // Read in the original freeze-out information
441// // and convert it from to [fm]
87ede832 442
443// // EPOS style
76ce4b5b 444// // fpx = mother->Xv()*1e13*0.197327;
445// // fpy = mother->Yv()*1e13*0.197327;
446// // fpz = mother->Zv()*1e13*0.197327;
447// // fpt = mother->T() *1e13*0.197327*0.5;
87ede832 448
449
450// // Therminator style
76ce4b5b 451// fpx = mother->Xv()*1e13;
452// fpy = mother->Yv()*1e13;
453// fpz = mother->Zv()*1e13;
454// fpt = mother->T() *1e13*3e10;
87ede832 455
76ce4b5b 456// }
457// }
87ede832 458
76ce4b5b 459// // if (fRotateToEventPlane) {
460// // double tPhi = TMath::ATan2(fpy, fpx);
461// // double tRad = TMath::Hypot(fpx, fpy);
87ede832 462
76ce4b5b 463// // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
464// // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
465// // }
466
467// tInfo->SetPDGPid(tPart->GetPdgCode());
468
469// // if (fRotateToEventPlane) {
470// // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
471// // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
87ede832 472
76ce4b5b 473// // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
474// // tRad*TMath::Sin(tPhi - tReactionPlane),
475// // tPart->Pz());
476// // }
477// // else
478// tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
479// Double_t mass2 = (tPart->E() *tPart->E() -
480// tPart->Px()*tPart->Px() -
481// tPart->Py()*tPart->Py() -
482// tPart->Pz()*tPart->Pz());
483// if (mass2>0.0)
484// tInfo->SetMass(TMath::Sqrt(mass2));
87ede832 485// else
76ce4b5b 486// tInfo->SetMass(0.0);
87ede832 487
76ce4b5b 488// tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
489// trackCopy->SetHiddenInfo(tInfo);
490
491// }
492
493// double pxyz[3];
494// aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
495// const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
496// // Check the sanity of the tracks - reject zero momentum tracks
497// if (ktP.Mag() == 0) {
498// delete trackCopy;
499// continue;
500// }
501// }
502// else {
87ede832 503 // No additional information exists
504 // Read in the normal AliAODTracks
76ce4b5b 505
87ede832 506 // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
507 AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
4eac0b05 508
509
64536eaf 510
87ede832 511 if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
64536eaf 512
87ede832 513 if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
514 delete trackCopy;
515 continue;
516 }
76ce4b5b 517
87ede832 518 if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
519 delete trackCopy;
520 continue;
521 }
7eb9f5d0 522
87ede832 523 // cout << "Muon? " << aodtrack->IsMuonTrack() << endl;
524 // cout << "Type? " << aodtrack->GetType() << endl;
525
526 // if (aodtrack->IsMuonTrack()) {
527 // cout << "muon" << endl;
528 // delete trackCopy;
529 // continue;
530 // }
531
532 //counting particles to set multiplicity
c684d45e 533 if(fEstEventMult==kGlobalCount){
534 double impact[2];
535 double covimpact[3];
536 if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
537 if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
538 //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
539 if (aodtrack->Chi2perNDF() < 4.0)
540 if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
541 if (aodtrack->GetTPCNcls() > 70)
542 if (aodtrack->Eta() < 0.8)
543 tNormMult++;
544 }
290ffcb9 545 }
7eb9f5d0 546
87ede832 547 CopyAODtoFemtoTrack(aodtrack, trackCopy);
548
549 // copying PID information from the correspondent track
550 // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
551
552
553 AliAODTrack *aodtrackpid;
554 if((fFilterBit == (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
555 aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
556 else
557 aodtrackpid = fEvent->GetTrack(i);
558 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
559
560 if (mcP) {
561 // Fill the hidden information with the simulated data
562 // Int_t pLabel = aodtrack->GetLabel();
563 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
564
565 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
566 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
567 if (!tPart) {
568 fpx = fV1[0];
569 fpy = fV1[1];
570 fpz = fV1[2];
571 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
572 tInfo->SetPDGPid(0);
573 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
574 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
575 tInfo->SetMass(0);
576 }
577 else {
578 // Check the mother information
579
580 // Using the new way of storing the freeze-out information
581 // Final state particle is stored twice on the stack
582 // one copy (mother) is stored with original freeze-out information
583 // and is not tracked
584 // the other one (daughter) is stored with primary vertex position
585 // and is tracked
586
587 // Freeze-out coordinates
588 fpx = tPart->Xv() - fV1[0];
589 fpy = tPart->Yv() - fV1[1];
590 fpz = tPart->Zv() - fV1[2];
591 // fpt = tPart->T();
592
593 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
594
50bb1be0 595
87ede832 596 fpx *= 1e13;
597 fpy *= 1e13;
598 fpz *= 1e13;
599 // fpt *= 1e13;
600
601 // cout << "Looking for mother ids " << endl;
602 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
603 // cout << "Got mother id" << endl;
604 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
605 // Check if this is the same particle stored twice on the stack
606 if (mother) {
607 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
608 // It is the same particle
609 // Read in the original freeze-out information
610 // and convert it from to [fm]
611
612 // EPOS style
613 // fpx = mother->Xv()*1e13*0.197327;
614 // fpy = mother->Yv()*1e13*0.197327;
615 // fpz = mother->Zv()*1e13*0.197327;
616 // fpt = mother->T() *1e13*0.197327*0.5;
617
618
619 // Therminator style
620 fpx = mother->Xv()*1e13;
621 fpy = mother->Yv()*1e13;
622 fpz = mother->Zv()*1e13;
623 // fpt = mother->T() *1e13*3e10;
624
625 }
626 }
627 }
628
629 // if (fRotateToEventPlane) {
630 // double tPhi = TMath::ATan2(fpy, fpx);
631 // double tRad = TMath::Hypot(fpx, fpy);
632
633 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
634 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
635 // }
636
637 tInfo->SetPDGPid(tPart->GetPdgCode());
638
639 // if (fRotateToEventPlane) {
640 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
641 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
642
643 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
644 // tRad*TMath::Sin(tPhi - tReactionPlane),
645 // tPart->Pz());
646 // }
647 // else
648 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
649 Double_t mass2 = (tPart->E() *tPart->E() -
650 tPart->Px()*tPart->Px() -
651 tPart->Py()*tPart->Py() -
652 tPart->Pz()*tPart->Pz());
653 if (mass2>0.0)
654 tInfo->SetMass(TMath::Sqrt(mass2));
655 else
656 tInfo->SetMass(0.0);
657
658 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
50bb1be0 659
ab558586 660 // // fillDCA
661 // //if (TMath::Abs(impact[0]) > 0.001) {
662 // if (tPart->IsPhysicalPrimary()){
663 // tInfo->SetPartOrigin(0);
664 // // trackCopy->SetImpactDprim(impact[0]);
290ffcb9 665 // //cout << "Read prim" << endl;
ab558586 666 // }
667 // else if (tPart->IsSecondaryFromWeakDecay()) {
668 // tInfo->SetPartOrigin(1);
669 // // trackCopy->SetImpactDweak(impact[0]);
670 // //cout << "Read wea" << endl;
671 // }
672 // else if (tPart->IsSecondaryFromMaterial()) {
673 // tInfo->SetPartOrigin(2);
674 // // trackCopy->SetImpactDmat(impact[0]);
675 // //cout << "Read mat" << endl;
676 // }
677 // //}
678 // // end fillDCA
50bb1be0 679
680
87ede832 681 }
682 trackCopy->SetHiddenInfo(tInfo);
76ce4b5b 683 }
87ede832 684
685 double pxyz[3];
686
687 //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
688 trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
689
690 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
691 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
692 // Check the sanity of the tracks - reject zero momentum tracks
693 if (ktP.Mag() == 0) {
694 delete trackCopy;
695 continue;
696 }
697 // }
698
290ffcb9 699 tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
87ede832 700 realnofTracks++;//real number of tracks
701 }
702
703 tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
76ce4b5b 704 tEvent->SetNormalizedMult(tracksPrim);
705
2bcd4b96 706 if (cent) {
707 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
18757d69 708 tEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
709 tEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
2bcd4b96 710 tEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
18757d69 711 tEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
2bcd4b96 712 tEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
18757d69 713 tEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
714 tEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
715 tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
716 tEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
717 // tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
607aa748 718 tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
719 tEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
2bcd4b96 720 }
721
1777446b 722 if (fEstEventMult==kCentrality) {
2bcd4b96 723 //AliCentrality *cent = fEvent->GetCentrality();
4eac0b05 724 //cout<<"AliFemtoEventReaderAOD:"<<lrint(10*cent->GetCentralityPercentile("V0M"))<<endl;
973a91f8 725 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
726 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
c863b24a 727 }
18757d69 728 else if (fEstEventMult==kCentralityV0A) {
729 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0A")));
730 }
731 else if (fEstEventMult==kCentralityV0C) {
732 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0C")));
733 }
c863b24a 734 else if (fEstEventMult==kCentralityZNA) {
c863b24a 735 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA")));
c863b24a 736 }
87ede832 737 else if (fEstEventMult==kCentralityZNC) {
18757d69 738 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNC")));
739 }
c863b24a 740 else if (fEstEventMult==kCentralityCL1) {
c863b24a 741 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL1")));
76ce4b5b 742 }
18757d69 743 else if (fEstEventMult==kCentralityCL0) {
744 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL0")));
745 }
607aa748 746 else if (fEstEventMult==kCentralityTRK) {
747 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TRK")));
748 }
18757d69 749 else if (fEstEventMult==kCentralityTKL) {
750 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TKL")));
751 }
607aa748 752 else if (fEstEventMult==kCentralityCND) {
753 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CND")));
754 }
18757d69 755 else if (fEstEventMult==kCentralityNPA) {
756 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("NPA")));
757 }
87ede832 758 else if (fEstEventMult==kCentralityFMD) {
18757d69 759 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("FMD")));
760 }
1777446b 761 else if(fEstEventMult==kGlobalCount){
ba3c23a4 762 tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
763 }
1777446b 764 else if(fEstEventMult==kReference)
87ede832 765 {
766 tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
767 }
1777446b 768 else if(fEstEventMult==kTPCOnlyRef)
87ede832 769 {
770 tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
771 }
1938eadf 772 else if(fEstEventMult == kVZERO)
87ede832 773 {
774 Float_t multV0 = 0;
775 for (Int_t i=0; i<64; i++)
776 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
777 tEvent->SetNormalizedMult(multV0);
778 }
76ce4b5b 779
780 if (mcP) delete [] motherids;
781
87ede832 782 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
973a91f8 783
784 if(fReadV0)
87ede832 785 {
786 int count_pass = 0;
787 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
788 AliAODv0* aodv0 = fEvent->GetV0(i);
789 if (!aodv0) continue;
790 if(aodv0->GetNDaughters()>2) continue;
791 if(aodv0->GetNProngs()>2) continue;
792 if(aodv0->GetCharge()!=0) continue;
793 if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
794 if(aodv0->CosPointingAngle(fV1)<0.998) continue;
795 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
796 count_pass++;
797 CopyAODtoFemtoV0(aodv0, trackCopyV0);
798 tEvent->V0Collection()->push_back(trackCopyV0);
799 //cout<<"Pushback v0 to v0collection"<<endl;
973a91f8 800 }
87ede832 801 }
973a91f8 802
76ce4b5b 803}
804
87ede832 805void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
806 AliFemtoTrack *tFemtoTrack
807 // AliPWG2AODTrack *tPWG2AODTrack
808 )
76ce4b5b 809{
810 // Copy the track information from the AOD into the internal AliFemtoTrack
811 // If it exists, use the additional information from the PWG2 AOD
812
813 // Primary Vertex position
87ede832 814
76ce4b5b 815 fEvent->GetPrimaryVertex()->GetPosition(fV1);
816 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
817
818 tFemtoTrack->SetCharge(tAodTrack->Charge());
87ede832 819
76ce4b5b 820 double pxyz[3];
821 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
290ffcb9 822
76ce4b5b 823 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
824 tFemtoTrack->SetP(v);//setting momentum
825 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
826 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
87ede832 827 //setting track helix
76ce4b5b 828 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
87ede832 829 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
76ce4b5b 830 tFemtoTrack->SetHelix(helix);
87ede832 831
76ce4b5b 832 // Flags
833 tFemtoTrack->SetTrackId(tAodTrack->GetID());
834 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
835 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
87ede832 836
837 // Track quality information
76ce4b5b 838 float covmat[6];
87ede832 839 tAodTrack->GetCovMatrix(covmat);
840
50bb1be0 841 if (!fDCAglobalTrack) {
842 double impact[2];
843 double covimpact[3];
76ce4b5b 844
50bb1be0 845 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
846 //cout << "sth went wrong with dca propagation" << endl;
847 tFemtoTrack->SetImpactD(-1000.0);
848 tFemtoTrack->SetImpactZ(-1000.0);
76ce4b5b 849
50bb1be0 850 }
851 else {
852 tFemtoTrack->SetImpactD(impact[0]);
853 tFemtoTrack->SetImpactZ(impact[1]);
854 }
855 }
76ce4b5b 856
857 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
858 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
859 // else
860 // tFemtoTrack->SetImpactD(0.0);
861 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
87ede832 862
76ce4b5b 863 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
864
865
866 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
867 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
868
869
87ede832 870 // cout
871 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
872 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
873 // << tAodTrack->Yv() - fV1[1]
874// << "xv = " << tAodTrack->Xv() << endl
875// << "fv1[0] = " << fV1[0] << endl
876// << "yv = " << tAodTrack->Yv() << endl
877// << "fv1[1] = " << fV1[1] << endl
878// << "zv = " << tAodTrack->Zv() << endl
879// << "fv1[2] = " << fV1[2] << endl
880// << "impact[0] = " << impact[0] << endl
881// << "impact[1] = " << impact[1] << endl
76ce4b5b 882// << endl << endl ;
883
884 tFemtoTrack->SetCdd(covmat[0]);
885 tFemtoTrack->SetCdz(covmat[1]);
886 tFemtoTrack->SetCzz(covmat[2]);
973a91f8 887 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
888 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
889 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
890 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
891 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
87ede832 892 tFemtoTrack->SetTPCsignalN(1);
893 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 894 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
895
896// if (tPWG2AODTrack) {
897// // Copy the PWG2 specific information if it exists
898// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
899// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
87ede832 900
76ce4b5b 901// double xtpc[3] = {0,0,0};
902// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
903// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
904// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
905// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
906// }
907// else {
87ede832 908 // If not use dummy values
76ce4b5b 909 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
910 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
87ede832 911
ce7b3d98 912 float globalPositionsAtRadii[9][3];
913 float bfield = 5*fMagFieldSign;
914 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
915 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
916 double **tpcPositions;
917 tpcPositions = new double*[9];
918 for(int i=0;i<9;i++)
919 tpcPositions[i] = new double[3];
920 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
921 for(int i=0;i<9;i++)
87ede832 922 {
923 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
924 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
925 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
926 }
7db4e4cf 927
ce7b3d98 928 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
929 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
930 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
87ede832 931 for(int i=0;i<9;i++)
932 delete [] tpcPositions[i];
933 delete [] tpcPositions;
76ce4b5b 934 // }
87ede832 935
76ce4b5b 936 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
87ede832 937
938
76ce4b5b 939 int indexes[3];
940 for (int ik=0; ik<3; ik++) {
941 indexes[ik] = 0;
942 }
943 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 944
945
946 for (int ii=0; ii<6; ii++){
947 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
948 }
949
950
76ce4b5b 951}
952
973a91f8 953void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
954{
ce7b3d98 955 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 956 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
957 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
958 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
959 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
960 tFemtoV0->SetdecayVertexV0(decayvertex);
961 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
962 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
963 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
964 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
965 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
966 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
967 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
968 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
969 tFemtoV0->SetmomPos(mompos);
970 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
971 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
972 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
973 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
974 tFemtoV0->SetmomNeg(momneg);
975
976 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
87ede832 977 //void SettpcHitsPos(const int& i);
978 //void SettpcHitsNeg(const int& i);
973a91f8 979
980 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
981 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
982
983 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
984 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
985 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
986 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
987 tFemtoV0->SetmomV0(momv0);
988 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
989 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
990 tFemtoV0->SeteLambda(tAODv0->ELambda());
991 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
992 tFemtoV0->SetePosProton(tAODv0->EPosProton());
993 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
994 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
995 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
996 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
997 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
998 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
87ede832 999
1000 //void SetcTauLambda( float x);
1001 //void SetcTauK0Short( float x);
1002
ce7b3d98 1003 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
1004 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 1005 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 1006 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
1007 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
1008 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
1009 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
87ede832 1010
973a91f8 1011 tFemtoV0->SetidNeg(tAODv0->GetNegID());
1012 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
1013 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1014 tFemtoV0->SetidPos(tAODv0->GetPosID());
1015
1016 tFemtoV0->SetEtaV0(tAODv0->Eta());
1017 tFemtoV0->SetPhiV0(tAODv0->Phi());
1018 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
1019 //tFemtoV0->SetYV0(tAODv0->Y());
1020
1021 //void SetdedxNeg(float x);
1022 //void SeterrdedxNeg(float x);//Gael 04Fev2002
1023 //void SetlendedxNeg(float x);//Gael 04Fev2002
1024 //void SetdedxPos(float x);
1025 //void SeterrdedxPos(float x);//Gael 04Fev2002
1026 //void SetlendedxPos(float x);//Gael 04Fev2002
1027
ce7b3d98 1028 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
1029 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 1030
1031 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
1032 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
1033
1034 if(trackpos && trackneg)
87ede832 1035 {
1036 tFemtoV0->SetEtaPos(trackpos->Eta());
1037 tFemtoV0->SetEtaNeg(trackneg->Eta());
1038 tFemtoV0->SetptotPos(tAODv0->PProng(0));
1039 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
1040 tFemtoV0->SetptPos(trackpos->Pt());
1041 tFemtoV0->SetptNeg(trackneg->Pt());
1042
1043 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
1044 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
1045 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1046 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1047 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1048 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1049 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1050 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1051 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
1052 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
1053 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1054 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1055
1056 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1057 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1058 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1059 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1060 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1061 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1062
1063
1064 float bfield = 5*fMagFieldSign;
1065 float globalPositionsAtRadiiPos[9][3];
1066 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1067 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1068 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1069
1070 float globalPositionsAtRadiiNeg[9][3];
1071 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1072 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1073 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1074
1075 AliFemtoThreeVector tmpVec;
290ffcb9 1076 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetY(tpcEntrancePos[1]); tmpVec.SetZ(tpcEntrancePos[2]);
87ede832 1077 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1078
290ffcb9 1079 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetY(tpcExitPos[1]); tmpVec.SetZ(tpcExitPos[2]);
87ede832 1080 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1081
290ffcb9 1082 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetY(tpcEntranceNeg[1]); tmpVec.SetZ(tpcEntranceNeg[2]);
87ede832 1083 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1084
290ffcb9 1085 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
87ede832 1086 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1087
1088 AliFemtoThreeVector vecTpcPos[9];
1089 AliFemtoThreeVector vecTpcNeg[9];
1090 for(int i=0;i<9;i++)
973a91f8 1091 {
87ede832 1092 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1093 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1094 }
1095 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1096 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
ce7b3d98 1097
87ede832 1098 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
1099 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
ce7b3d98 1100
87ede832 1101 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1102 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
ce7b3d98 1103
87ede832 1104 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
1105 {
1106 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
973a91f8 1107 {
ce7b3d98 1108 tFemtoV0->SetPosNSigmaTOFK(-1000);
1109 tFemtoV0->SetNegNSigmaTOFK(-1000);
1110 tFemtoV0->SetPosNSigmaTOFP(-1000);
1111 tFemtoV0->SetNegNSigmaTOFP(-1000);
1112 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1113 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1114
1115 tFemtoV0->SetTOFProtonTimePos(-1000);
1116 tFemtoV0->SetTOFPionTimePos(-1000);
1117 tFemtoV0->SetTOFKaonTimePos(-1000);
1118 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1119 tFemtoV0->SetTOFPionTimeNeg(-1000);
1120 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1121 }
973a91f8 1122 }
87ede832 1123 else
973a91f8 1124 {
87ede832 1125 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1126 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1127 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1128 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1129 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1130 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1131
1132 double TOFSignalPos = trackpos->GetTOFsignal();
1133 double TOFSignalNeg = trackneg->GetTOFsignal();
1134 double pidPos[5];
1135 double pidNeg[5];
1136 trackpos->GetIntegratedTimes(pidPos);
1137 trackneg->GetIntegratedTimes(pidNeg);
1138
1139 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1140 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1141 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1142 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1143 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1144 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1145 }
87ede832 1146 }
1147 else
1148 {
1149 tFemtoV0->SetStatusPos(999);
1150 tFemtoV0->SetStatusNeg(999);
1151 }
973a91f8 1152 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1153}
1154
76ce4b5b 1155void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1156{
1157 fFilterBit = (1 << (ibit));
1158}
1159
64536eaf 1160
1161void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1162{
1163 fFilterMask = ibit;
1164}
1165
76ce4b5b 1166void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1167{
1168 fReadMC = a;
1169}
1170
973a91f8 1171
1172void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1173{
1174 fReadV0 = a;
1175}
1176
1777446b 1177void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1178{
1179 fEstEventMult = aType;
1180}
1181
76ce4b5b 1182AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1183{
1184 if (aLabel < 0) return 0;
1185 AliAODMCParticle *aodP;
1186 Int_t posstack = 0;
1187 if (aLabel > mcP->GetEntries())
1188 posstack = mcP->GetEntries();
1189 else
1190 posstack = aLabel;
1191
1192 aodP = (AliAODMCParticle *) mcP->At(posstack);
1193 if (aodP->GetLabel() > posstack) {
1194 do {
1195 aodP = (AliAODMCParticle *) mcP->At(posstack);
1196 if (aodP->GetLabel() == aLabel) return aodP;
1197 posstack--;
1198 }
1199 while (posstack > 0);
1200 }
1201 else {
1202 do {
1203 aodP = (AliAODMCParticle *) mcP->At(posstack);
1204 if (aodP->GetLabel() == aLabel) return aodP;
1205 posstack++;
1206 }
1207 while (posstack < mcP->GetEntries());
1208 }
87ede832 1209
76ce4b5b 1210 return 0;
1211}
1212
87ede832 1213void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1214 AliFemtoTrack *tFemtoTrack)
76ce4b5b 1215{
41cb6c1e 1216
50bb1be0 1217 if (fDCAglobalTrack) {
1218 // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
41cb6c1e 1219
50bb1be0 1220 double impact[2];
1221 double covimpact[3];
6691ea4f 1222
50bb1be0 1223 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1224 //cout << "sth went wrong with dca propagation" << endl;
1225 tFemtoTrack->SetImpactD(-1000.0);
1226 tFemtoTrack->SetImpactZ(-1000.0);
6691ea4f 1227
50bb1be0 1228 }
1229 else {
1230 tFemtoTrack->SetImpactD(impact[0]);
1231 tFemtoTrack->SetImpactZ(impact[1]);
1232 }
1233 }
41cb6c1e 1234
76ce4b5b 1235 double aodpid[10];
1236 tAodTrack->GetPID(aodpid);
1237 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1238 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1239 tFemtoTrack->SetPidProbPion(aodpid[2]);
1240 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1241 tFemtoTrack->SetPidProbProton(aodpid[4]);
1242
1243 aodpid[0] = -100000.0;
1244 aodpid[1] = -100000.0;
1245 aodpid[2] = -100000.0;
1246 aodpid[3] = -100000.0;
1247 aodpid[4] = -100000.0;
87ede832 1248
76ce4b5b 1249 double tTOF = 0.0;
87ede832 1250 Float_t probMis = 1.0;
76ce4b5b 1251
4eac0b05 1252 //what is that code? for what do we need it? nsigma values are not enough?
87ede832 1253 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1254 tTOF = tAodTrack->GetTOFsignal();
1255 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1256
87ede832 1257 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1258
1259 probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
1260 }
1261
1262
1263
1264 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1265
76ce4b5b 1266 ////// TPC ////////////////////////////////////////////
1267
87ede832 1268 float nsigmaTPCK=-1000.;
1269 float nsigmaTPCPi=-1000.;
1270 float nsigmaTPCP=-1000.;
1271
76ce4b5b 1272 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1273
1274 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1275 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1276 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1277 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1278 }
1279
1280 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1281 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1282 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
1283
87ede832 1284 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1285 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1286 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1287
1288 tFemtoTrack->SetTPCsignalN(1);
1289 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 1290 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
76ce4b5b 1291
87ede832 1292 ///////TOF//////////////////////
76ce4b5b 1293
87ede832 1294 float vp=-1000.;
1295 float nsigmaTOFPi=-1000.;
1296 float nsigmaTOFK=-1000.;
1297 float nsigmaTOFP=-1000.;
1298
1299 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1300 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
1301 (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
1302 && probMis < 0.01) // TOF mismatch probaility
1303 {
1304 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
76ce4b5b 1305 {
1306
1307 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1308 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1309 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1310
1311 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1312 Double_t tof=tAodTrack->GetTOFsignal();
1313 if(tof > 0.) vp=len/tof/0.03;
1314 }
87ede832 1315 }
1316 tFemtoTrack->SetVTOF(vp);
1317 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1318 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1319 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1320
76ce4b5b 1321
87ede832 1322 //////////////////////////////////////
76ce4b5b 1323
1324}
1325
1326void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1327{
1328 fCentRange[0] = min; fCentRange[1] = max;
87ede832 1329 fUsePreCent = 1;
1777446b 1330 fEstEventMult = kCentrality;
76ce4b5b 1331}
1332
973a91f8 1333void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1334{
1777446b 1335 if(anocent==false) {fEstEventMult=kCentrality;}
1336 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1337}
76ce4b5b 1338
1339void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1340{
1341 fAODpidUtil = aAODpidUtil;
1342 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1343}
1344
1777446b 1345void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1346{
1347 fAODheader = aAODheader;
1348}
1349
76ce4b5b 1350
ba3c23a4 1351void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1352{
1353 if(s>0)
1354 fMagFieldSign = 1;
1355 else if(s<0)
1356 fMagFieldSign = -1;
1357 else
1358 fMagFieldSign = 0;
1359}
1360
5e2038a9 1361void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1362{
87ede832 1363 fisEPVZ = iepvz;
5e2038a9 1364}
76ce4b5b 1365
ce7b3d98 1366void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1367{
1368 // Gets the global position of the track at nine different radii in the TPC
1369 // track is the track you want to propagate
1370 // bfield is the magnetic field of your event
1371 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1372
1373 // Initialize the array to something indicating there was no propagation
1374 for(Int_t i=0;i<9;i++){
1375 for(Int_t j=0;j<3;j++){
1376 globalPositionsAtRadii[i][j]=-9999.;
1377 }
1378 }
1379
1380 // Make a copy of the track to not change parameters of the track
1381 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1382 //printf("\nAfter CopyFromVTrack\n");
1383 //etp.Print();
1384
1385 // The global position of the the track
1386 Double_t xyz[3]={-9999.,-9999.,-9999.};
1387
1388 // Counter for which radius we want
1389 Int_t iR=0;
1390 // The radii at which we get the global positions
1391 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1392 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1393 // The global radius we are at
1394 Float_t globalRadius=0;
1395
1396 // Propagation is done in local x of the track
1397 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1398 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1399 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1400 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1401 // adds to the global radius
1402
1403 // Stop if the propagation was not succesful. This can happen for low pt tracks
1404 // that don't reach outer radii
1405 if(!etp.PropagateTo(x,bfield))break;
1406 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1407 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1408
1409 // Roughly reached the radius we want
1410 if(globalRadius > Rwanted[iR]){
1411
1412 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1413 while (globalRadius>Rwanted[iR]){
87ede832 1414 x-=.1;
1415 // printf("propagating to x %5.2f\n",x);
1416 if(!etp.PropagateTo(x,bfield))break;
1417 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1418 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1419 }
1420 //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]);
1421 globalPositionsAtRadii[iR][0]=xyz[0];
1422 globalPositionsAtRadii[iR][1]=xyz[1];
1423 globalPositionsAtRadii[iR][2]=xyz[2];
1424 // Indicate we want the next radius
1425 iR+=1;
1426 }
1427 if(iR>=8){
1428 // TPC edge reached
1429 return;
1430 }
1431 }
1432}
76ce4b5b 1433
734c121a 1434void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
1435{
1436 fpA2013 = pa2013;
1437}
50bb1be0 1438
1439void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
1440{
1441 fDCAglobalTrack = dcagt;
1442}
1443