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