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