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