]> 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 in the average separation calculation
[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];
c710fd31 913
914 float s = fEvent->GetMagneticField();
915 if(s>0)
916 fMagFieldSign = 1;
917 else if(s<0)
918 fMagFieldSign = -1;
919 else
920 fMagFieldSign = 0;
921
ce7b3d98 922 float bfield = 5*fMagFieldSign;
c710fd31 923
ce7b3d98 924 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
925 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
926 double **tpcPositions;
927 tpcPositions = new double*[9];
928 for(int i=0;i<9;i++)
929 tpcPositions[i] = new double[3];
930 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
931 for(int i=0;i<9;i++)
87ede832 932 {
933 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
934 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
935 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
936 }
7db4e4cf 937
ce7b3d98 938 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
939 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
940 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
87ede832 941 for(int i=0;i<9;i++)
942 delete [] tpcPositions[i];
943 delete [] tpcPositions;
76ce4b5b 944 // }
87ede832 945
76ce4b5b 946 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
87ede832 947
948
76ce4b5b 949 int indexes[3];
950 for (int ik=0; ik<3; ik++) {
951 indexes[ik] = 0;
952 }
953 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 954
955
956 for (int ii=0; ii<6; ii++){
957 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
958 }
959
960
76ce4b5b 961}
962
973a91f8 963void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
964{
ce7b3d98 965 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 966 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
967 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
968 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
969 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
970 tFemtoV0->SetdecayVertexV0(decayvertex);
971 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
972 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
973 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
974 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
975 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
976 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
977 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
978 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
979 tFemtoV0->SetmomPos(mompos);
980 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
981 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
982 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
983 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
984 tFemtoV0->SetmomNeg(momneg);
985
986 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
87ede832 987 //void SettpcHitsPos(const int& i);
988 //void SettpcHitsNeg(const int& i);
973a91f8 989
990 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
991 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
992
993 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
994 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
995 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
996 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
997 tFemtoV0->SetmomV0(momv0);
998 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
999 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
1000 tFemtoV0->SeteLambda(tAODv0->ELambda());
1001 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
1002 tFemtoV0->SetePosProton(tAODv0->EPosProton());
1003 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
1004 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
1005 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
1006 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
1007 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
1008 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
87ede832 1009
1010 //void SetcTauLambda( float x);
1011 //void SetcTauK0Short( float x);
1012
ce7b3d98 1013 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
1014 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 1015 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 1016 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
1017 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
1018 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
1019 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
87ede832 1020
973a91f8 1021 tFemtoV0->SetidNeg(tAODv0->GetNegID());
1022 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
1023 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1024 tFemtoV0->SetidPos(tAODv0->GetPosID());
1025
1026 tFemtoV0->SetEtaV0(tAODv0->Eta());
1027 tFemtoV0->SetPhiV0(tAODv0->Phi());
1028 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
1029 //tFemtoV0->SetYV0(tAODv0->Y());
1030
1031 //void SetdedxNeg(float x);
1032 //void SeterrdedxNeg(float x);//Gael 04Fev2002
1033 //void SetlendedxNeg(float x);//Gael 04Fev2002
1034 //void SetdedxPos(float x);
1035 //void SeterrdedxPos(float x);//Gael 04Fev2002
1036 //void SetlendedxPos(float x);//Gael 04Fev2002
1037
ce7b3d98 1038 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
1039 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 1040
1041 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
1042 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
1043
1044 if(trackpos && trackneg)
87ede832 1045 {
1046 tFemtoV0->SetEtaPos(trackpos->Eta());
1047 tFemtoV0->SetEtaNeg(trackneg->Eta());
1048 tFemtoV0->SetptotPos(tAODv0->PProng(0));
1049 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
1050 tFemtoV0->SetptPos(trackpos->Pt());
1051 tFemtoV0->SetptNeg(trackneg->Pt());
1052
1053 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
1054 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
1055 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1056 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1057 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1058 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1059 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1060 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1061 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
1062 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
1063 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1064 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1065
1066 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1067 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1068 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1069 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1070 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1071 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1072
1073
1074 float bfield = 5*fMagFieldSign;
1075 float globalPositionsAtRadiiPos[9][3];
1076 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1077 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1078 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1079
1080 float globalPositionsAtRadiiNeg[9][3];
1081 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1082 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1083 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1084
1085 AliFemtoThreeVector tmpVec;
290ffcb9 1086 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetY(tpcEntrancePos[1]); tmpVec.SetZ(tpcEntrancePos[2]);
87ede832 1087 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1088
290ffcb9 1089 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetY(tpcExitPos[1]); tmpVec.SetZ(tpcExitPos[2]);
87ede832 1090 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1091
290ffcb9 1092 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetY(tpcEntranceNeg[1]); tmpVec.SetZ(tpcEntranceNeg[2]);
87ede832 1093 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1094
290ffcb9 1095 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
87ede832 1096 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1097
1098 AliFemtoThreeVector vecTpcPos[9];
1099 AliFemtoThreeVector vecTpcNeg[9];
1100 for(int i=0;i<9;i++)
973a91f8 1101 {
87ede832 1102 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1103 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1104 }
1105 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1106 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
ce7b3d98 1107
87ede832 1108 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
1109 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
ce7b3d98 1110
87ede832 1111 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1112 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
ce7b3d98 1113
87ede832 1114 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
1115 {
1116 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
973a91f8 1117 {
ce7b3d98 1118 tFemtoV0->SetPosNSigmaTOFK(-1000);
1119 tFemtoV0->SetNegNSigmaTOFK(-1000);
1120 tFemtoV0->SetPosNSigmaTOFP(-1000);
1121 tFemtoV0->SetNegNSigmaTOFP(-1000);
1122 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1123 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1124
1125 tFemtoV0->SetTOFProtonTimePos(-1000);
1126 tFemtoV0->SetTOFPionTimePos(-1000);
1127 tFemtoV0->SetTOFKaonTimePos(-1000);
1128 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1129 tFemtoV0->SetTOFPionTimeNeg(-1000);
1130 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1131 }
973a91f8 1132 }
87ede832 1133 else
973a91f8 1134 {
87ede832 1135 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1136 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1137 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1138 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1139 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1140 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1141
1142 double TOFSignalPos = trackpos->GetTOFsignal();
1143 double TOFSignalNeg = trackneg->GetTOFsignal();
1144 double pidPos[5];
1145 double pidNeg[5];
1146 trackpos->GetIntegratedTimes(pidPos);
1147 trackneg->GetIntegratedTimes(pidNeg);
1148
1149 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1150 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1151 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1152 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1153 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1154 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1155 }
87ede832 1156 }
1157 else
1158 {
1159 tFemtoV0->SetStatusPos(999);
1160 tFemtoV0->SetStatusNeg(999);
1161 }
973a91f8 1162 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1163}
1164
76ce4b5b 1165void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1166{
1167 fFilterBit = (1 << (ibit));
1168}
1169
64536eaf 1170
1171void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1172{
1173 fFilterMask = ibit;
1174}
1175
76ce4b5b 1176void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1177{
1178 fReadMC = a;
1179}
1180
973a91f8 1181
1182void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1183{
1184 fReadV0 = a;
1185}
1186
1777446b 1187void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1188{
1189 fEstEventMult = aType;
1190}
1191
76ce4b5b 1192AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1193{
1194 if (aLabel < 0) return 0;
1195 AliAODMCParticle *aodP;
1196 Int_t posstack = 0;
1197 if (aLabel > mcP->GetEntries())
1198 posstack = mcP->GetEntries();
1199 else
1200 posstack = aLabel;
1201
1202 aodP = (AliAODMCParticle *) mcP->At(posstack);
1203 if (aodP->GetLabel() > posstack) {
1204 do {
1205 aodP = (AliAODMCParticle *) mcP->At(posstack);
1206 if (aodP->GetLabel() == aLabel) return aodP;
1207 posstack--;
1208 }
1209 while (posstack > 0);
1210 }
1211 else {
1212 do {
1213 aodP = (AliAODMCParticle *) mcP->At(posstack);
1214 if (aodP->GetLabel() == aLabel) return aodP;
1215 posstack++;
1216 }
1217 while (posstack < mcP->GetEntries());
1218 }
87ede832 1219
76ce4b5b 1220 return 0;
1221}
1222
87ede832 1223void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1224 AliFemtoTrack *tFemtoTrack)
76ce4b5b 1225{
41cb6c1e 1226
50bb1be0 1227 if (fDCAglobalTrack) {
1228 // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
41cb6c1e 1229
50bb1be0 1230 double impact[2];
1231 double covimpact[3];
6691ea4f 1232
50bb1be0 1233 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1234 //cout << "sth went wrong with dca propagation" << endl;
1235 tFemtoTrack->SetImpactD(-1000.0);
1236 tFemtoTrack->SetImpactZ(-1000.0);
6691ea4f 1237
50bb1be0 1238 }
1239 else {
1240 tFemtoTrack->SetImpactD(impact[0]);
1241 tFemtoTrack->SetImpactZ(impact[1]);
1242 }
1243 }
41cb6c1e 1244
76ce4b5b 1245 double aodpid[10];
1246 tAodTrack->GetPID(aodpid);
1247 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1248 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1249 tFemtoTrack->SetPidProbPion(aodpid[2]);
1250 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1251 tFemtoTrack->SetPidProbProton(aodpid[4]);
1252
1253 aodpid[0] = -100000.0;
1254 aodpid[1] = -100000.0;
1255 aodpid[2] = -100000.0;
1256 aodpid[3] = -100000.0;
1257 aodpid[4] = -100000.0;
87ede832 1258
76ce4b5b 1259 double tTOF = 0.0;
87ede832 1260 Float_t probMis = 1.0;
76ce4b5b 1261
4eac0b05 1262 //what is that code? for what do we need it? nsigma values are not enough?
87ede832 1263 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1264 tTOF = tAodTrack->GetTOFsignal();
1265 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1266
87ede832 1267 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1268
1269 probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
1270 }
1271
1272
1273
1274 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1275
76ce4b5b 1276 ////// TPC ////////////////////////////////////////////
1277
87ede832 1278 float nsigmaTPCK=-1000.;
1279 float nsigmaTPCPi=-1000.;
1280 float nsigmaTPCP=-1000.;
1281
76ce4b5b 1282 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1283
1284 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1285 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1286 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1287 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1288 }
1289
1290 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1291 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1292 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
1293
87ede832 1294 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1295 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1296 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1297
1298 tFemtoTrack->SetTPCsignalN(1);
1299 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 1300 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
76ce4b5b 1301
87ede832 1302 ///////TOF//////////////////////
76ce4b5b 1303
87ede832 1304 float vp=-1000.;
1305 float nsigmaTOFPi=-1000.;
1306 float nsigmaTOFK=-1000.;
1307 float nsigmaTOFP=-1000.;
1308
1309 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1310 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
1311 (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
1312 && probMis < 0.01) // TOF mismatch probaility
1313 {
1314 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
76ce4b5b 1315 {
1316
1317 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1318 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1319 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1320
1321 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1322 Double_t tof=tAodTrack->GetTOFsignal();
1323 if(tof > 0.) vp=len/tof/0.03;
1324 }
87ede832 1325 }
1326 tFemtoTrack->SetVTOF(vp);
1327 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1328 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1329 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1330
76ce4b5b 1331
87ede832 1332 //////////////////////////////////////
76ce4b5b 1333
1334}
1335
1336void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1337{
1338 fCentRange[0] = min; fCentRange[1] = max;
87ede832 1339 fUsePreCent = 1;
1777446b 1340 fEstEventMult = kCentrality;
76ce4b5b 1341}
1342
973a91f8 1343void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1344{
1777446b 1345 if(anocent==false) {fEstEventMult=kCentrality;}
1346 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1347}
76ce4b5b 1348
1349void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1350{
1351 fAODpidUtil = aAODpidUtil;
1352 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1353}
1354
1777446b 1355void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1356{
1357 fAODheader = aAODheader;
1358}
1359
76ce4b5b 1360
ba3c23a4 1361void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1362{
1363 if(s>0)
1364 fMagFieldSign = 1;
1365 else if(s<0)
1366 fMagFieldSign = -1;
1367 else
1368 fMagFieldSign = 0;
1369}
1370
5e2038a9 1371void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1372{
87ede832 1373 fisEPVZ = iepvz;
5e2038a9 1374}
76ce4b5b 1375
ce7b3d98 1376void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1377{
1378 // Gets the global position of the track at nine different radii in the TPC
1379 // track is the track you want to propagate
1380 // bfield is the magnetic field of your event
1381 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1382
1383 // Initialize the array to something indicating there was no propagation
1384 for(Int_t i=0;i<9;i++){
1385 for(Int_t j=0;j<3;j++){
1386 globalPositionsAtRadii[i][j]=-9999.;
1387 }
1388 }
1389
1390 // Make a copy of the track to not change parameters of the track
1391 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1392 //printf("\nAfter CopyFromVTrack\n");
1393 //etp.Print();
1394
1395 // The global position of the the track
1396 Double_t xyz[3]={-9999.,-9999.,-9999.};
1397
1398 // Counter for which radius we want
1399 Int_t iR=0;
1400 // The radii at which we get the global positions
1401 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1402 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1403 // The global radius we are at
1404 Float_t globalRadius=0;
1405
1406 // Propagation is done in local x of the track
1407 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1408 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1409 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1410 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1411 // adds to the global radius
1412
1413 // Stop if the propagation was not succesful. This can happen for low pt tracks
1414 // that don't reach outer radii
1415 if(!etp.PropagateTo(x,bfield))break;
1416 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1417 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1418
1419 // Roughly reached the radius we want
1420 if(globalRadius > Rwanted[iR]){
1421
1422 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1423 while (globalRadius>Rwanted[iR]){
87ede832 1424 x-=.1;
1425 // printf("propagating to x %5.2f\n",x);
1426 if(!etp.PropagateTo(x,bfield))break;
1427 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1428 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1429 }
1430 //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]);
1431 globalPositionsAtRadii[iR][0]=xyz[0];
1432 globalPositionsAtRadii[iR][1]=xyz[1];
1433 globalPositionsAtRadii[iR][2]=xyz[2];
1434 // Indicate we want the next radius
1435 iR+=1;
1436 }
1437 if(iR>=8){
1438 // TPC edge reached
1439 return;
1440 }
1441 }
1442}
76ce4b5b 1443
734c121a 1444void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
1445{
1446 fpA2013 = pa2013;
1447}
50bb1be0 1448
1449void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
1450{
1451 fDCAglobalTrack = dcagt;
1452}
1453