Add Reaction Plane aware analysis
[u/mrichter/AliRoot.git] / PWG2 / 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
19 #include "AliFmPhysicalHelixD.h"
20 #include "AliFmThreeVectorF.h"
21
22 #include "SystemOfUnits.h"
23
24 #include "AliFemtoEvent.h"
25 #include "AliFemtoModelHiddenInfo.h"
26 #include "AliFemtoModelGlobalHiddenInfo.h"
27
28 ClassImp(AliFemtoEventReaderAOD)
29
30 #if !(ST_NO_NAMESPACES)
31   using namespace units;
32 #endif
33
34 using namespace std;
35 //____________________________
36 //constructor with 0 parameters , look at default settings 
37 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
38   fNumberofEvent(0),
39   fCurEvent(0),
40   fEvent(0x0),
41   fAllTrue(160),
42   fAllFalse(160),
43   fFilterBit(0),
44   fPWG2AODTracks(0x0),
45   fReadMC(0),
46   fInputFile(" "),
47   fFileName(" "),
48   fTree(0x0),
49   fAodFile(0x0)
50 {
51   // default constructor
52   fAllTrue.ResetAllBits(kTRUE);
53   fAllFalse.ResetAllBits(kFALSE);
54 }
55
56 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
57   AliFemtoEventReader(),
58   fNumberofEvent(0),
59   fCurEvent(0),
60   fEvent(0x0),
61   fAllTrue(160),
62   fAllFalse(160),
63   fFilterBit(0),
64   fPWG2AODTracks(0x0),
65   fReadMC(0),
66   fInputFile(" "),
67   fFileName(" "),
68   fTree(0x0),
69   fAodFile(0x0)
70 {
71   // copy constructor
72   fInputFile = aReader.fInputFile;
73   fFileName  = aReader.fFileName;
74   fNumberofEvent = aReader.fNumberofEvent;
75   fCurEvent = aReader.fCurEvent;
76   fEvent = new AliAODEvent();
77   fAodFile = new TFile(aReader.fAodFile->GetName());
78   fAllTrue.ResetAllBits(kTRUE);
79   fAllFalse.ResetAllBits(kFALSE);
80   fFilterBit = aReader.fFilterBit;
81   fPWG2AODTracks = aReader.fPWG2AODTracks;
82 }
83 //__________________
84 //Destructor
85 AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
86 {
87   // destructor
88   delete fTree;
89   delete fEvent;
90   delete fAodFile;
91   if (fPWG2AODTracks) {
92     fPWG2AODTracks->Delete();
93     delete fPWG2AODTracks;
94   }
95 }
96
97 //__________________
98 AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
99 {
100   // assignment operator
101   if (this == &aReader)
102     return *this;
103
104   fInputFile = aReader.fInputFile;
105   fFileName  = aReader.fFileName;
106   fNumberofEvent = aReader.fNumberofEvent;
107   fCurEvent = aReader.fCurEvent;
108   if (fTree) delete fTree;
109   if (fEvent) delete fEvent;
110   fEvent = new AliAODEvent();
111   if (fAodFile) delete fAodFile;
112   fAodFile = new TFile(aReader.fAodFile->GetName());
113   fAllTrue.ResetAllBits(kTRUE);
114   fAllFalse.ResetAllBits(kFALSE);
115   fFilterBit = aReader.fFilterBit;
116   fPWG2AODTracks = aReader.fPWG2AODTracks;
117
118   return *this;
119 }
120 //__________________
121 AliFemtoString AliFemtoEventReaderAOD::Report()
122 {
123   // create reader report
124   AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
125   return temp;
126 }
127
128 //__________________
129 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
130 {
131   //setting the name of file where names of AOD file are written 
132   //it takes only this files which have good trees
133   char buffer[256];
134   fInputFile=string(inputFile);
135   ifstream infile(inputFile);
136
137   fTree = new TChain("aodTree");
138
139   if(infile.good()==true)
140     { 
141       //checking if all give files have good tree inside
142       while (infile.eof()==false)
143         {
144           infile.getline(buffer,256);
145           TFile *aodFile=TFile::Open(buffer,"READ");
146           if (aodFile!=0x0)
147             {   
148               TTree* tree = (TTree*) aodFile->Get("aodTree");
149               if (tree!=0x0)
150                 {
151 //                cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
152                   fTree->AddFile(buffer);
153                   delete tree;
154                 }
155               aodFile->Close(); 
156             }
157           delete aodFile;
158         }
159     }
160 }
161
162 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
163 {
164   // read in a next hbt event from the chain
165   // convert it to AliFemtoEvent and return
166   // for further analysis
167   AliFemtoEvent *hbtEvent = 0;
168
169   if (fCurEvent==fNumberofEvent)//open next file  
170     {
171       if(fNumberofEvent==0)     
172         {
173           fEvent=new AliAODEvent();
174           fEvent->ReadFromTree(fTree);
175
176           // Check for the existence of the additional information
177           fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
178
179           if (fPWG2AODTracks) {
180             cout << "Found additional PWG2 specific information in the AOD!" << endl;
181             cout << "Reading only tracks with the additional information" << endl;
182           }
183
184           fNumberofEvent=fTree->GetEntries();
185           //      cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
186           fCurEvent=0;
187         }
188       else //no more data to read
189         {
190           cout<<"no more files "<<hbtEvent<<endl;
191           fReaderStatus=1;
192           return hbtEvent; 
193         }
194     }           
195
196   cout<<"starting to read event "<<fCurEvent<<endl;
197   fTree->GetEvent(fCurEvent);//getting next event
198   //  cout << "Read event " << fEvent << " from file " << fTree << endl;
199         
200   hbtEvent = new AliFemtoEvent;
201
202   CopyAODtoFemtoEvent(hbtEvent);
203
204   fCurEvent++;  
205
206   return hbtEvent; 
207 }
208
209 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
210 {
211   // A function that reads in the AOD event
212   // and transfers the neccessary information into
213   // the internal AliFemtoEvent
214
215   // setting global event characteristics
216   tEvent->SetRunNumber(fEvent->GetRunNumber());
217   tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
218   tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
219   tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
220   tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
221   tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
222   tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
223   tEvent->SetZDCParticipants(0);
224   tEvent->SetTriggerMask(fEvent->GetTriggerMask());
225   tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
226
227   // Attempt to access MC header
228   AliAODMCHeader *mcH;
229   TClonesArray *mcP;
230   if (fReadMC) {
231     mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
232     if (!mcH) {
233       cout << "AOD MC information requested, but no header found!" << endl;
234     }
235
236     mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
237     if (!mcP) {
238       cout << "AOD MC information requested, but no particle array found!" << endl;
239     }
240   }
241
242   tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
243
244   Int_t *motherids;
245   if (mcP) {
246     motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
247     for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
248
249     // Read in mother ids
250     AliAODMCParticle *motherpart;
251     for (int ip=0; ip<mcP->GetEntries(); ip++) {
252       motherpart = (AliAODMCParticle *) mcP->At(ip);
253       if (motherpart->GetDaughter(0) > 0)
254         motherids[motherpart->GetDaughter(0)] = ip;
255       if (motherpart->GetDaughter(1) > 0)
256         motherids[motherpart->GetDaughter(1)] = ip;
257     }
258   }
259
260   // Primary Vertex position
261   double fV1[3];
262   fEvent->GetPrimaryVertex()->GetPosition(fV1);
263
264   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
265   tEvent->SetPrimVertPos(vertex);
266         
267   //starting to reading tracks
268   int nofTracks=0;  //number of reconstructed tracks in event
269
270   // Check to see whether the additional info exists
271   if (fPWG2AODTracks)
272     nofTracks=fPWG2AODTracks->GetEntries();
273   else
274     nofTracks=fEvent->GetNumberOfTracks();
275
276   int realnofTracks=0;   // number of track which we use in a analysis
277
278   for (int i=0;i<nofTracks;i++)
279     {
280       AliFemtoTrack* trackCopy = new AliFemtoTrack();   
281
282       if (fPWG2AODTracks) {
283         // Read tracks from the additional pwg2 specific AOD part
284         // if they exist
285         // Note that in that case all the AOD tracks without the 
286         // additional information will be ignored !
287         AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
288
289         // Getting the AOD track through the ref of the additional info
290         AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack(); 
291         if (!aodtrack->TestFilterBit(fFilterBit))
292           continue;
293         
294         CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
295         
296         if (mcP) {
297           // Fill the hidden information with the simulated data
298           Int_t pLabel = aodtrack->GetLabel();
299           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
300
301           // Check the mother information
302           
303           // Using the new way of storing the freeze-out information
304           // Final state particle is stored twice on the stack
305           // one copy (mother) is stored with original freeze-out information
306           //   and is not tracked
307           // the other one (daughter) is stored with primary vertex position
308           //   and is tracked
309           
310           // Freeze-out coordinates
311           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
312           fpx = tPart->Xv() - fV1[0];
313           fpy = tPart->Yv() - fV1[1];
314           fpz = tPart->Zv() - fV1[2];
315           fpt = tPart->T();
316
317           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
318           tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
319
320           fpx *= 1e13;
321           fpy *= 1e13;
322           fpz *= 1e13;
323           fpt *= 1e13;
324           
325           //      cout << "Looking for mother ids " << endl;
326           if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
327             //  cout << "Got mother id" << endl;
328             AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
329             // Check if this is the same particle stored twice on the stack
330             if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
331               // It is the same particle
332               // Read in the original freeze-out information
333               // and convert it from to [fm]
334               
335               // EPOS style 
336               //          fpx = mother->Xv()*1e13*0.197327;
337               //          fpy = mother->Yv()*1e13*0.197327;
338               //          fpz = mother->Zv()*1e13*0.197327;
339               //          fpt = mother->T() *1e13*0.197327*0.5;
340               
341               
342               // Therminator style 
343               fpx = mother->Xv()*1e13;
344               fpy = mother->Yv()*1e13;
345               fpz = mother->Zv()*1e13;
346               fpt = mother->T() *1e13*3e10;
347               
348             }
349           }
350           
351 //       if (fRotateToEventPlane) {
352 //      double tPhi = TMath::ATan2(fpy, fpx);
353 //      double tRad = TMath::Hypot(fpx, fpy);
354         
355 //      fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
356 //      fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
357 //       }
358
359           tInfo->SetPDGPid(tPart->GetPdgCode());
360
361 //        if (fRotateToEventPlane) {
362 //          double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
363 //          double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
364             
365 //          tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
366 //                                 tRad*TMath::Sin(tPhi - tReactionPlane),
367 //                                 tPart->Pz());
368 //        }
369 //       else
370           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
371           Double_t mass2 = (tPart->E() *tPart->E() -
372                             tPart->Px()*tPart->Px() -
373                             tPart->Py()*tPart->Py() -
374                             tPart->Pz()*tPart->Pz());
375           if (mass2>0.0)
376             tInfo->SetMass(TMath::Sqrt(mass2));
377           else 
378             tInfo->SetMass(0.0);
379           
380           tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
381           trackCopy->SetHiddenInfo(tInfo);
382
383         }
384
385         double pxyz[3];
386         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
387         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
388         // Check the sanity of the tracks - reject zero momentum tracks
389         if (ktP.mag() == 0) {
390           delete trackCopy;
391           continue;
392         }
393       }
394       else {
395         // No additional information exists
396         // Read in the normal AliAODTracks 
397         const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
398         
399 //      if (!aodtrack->TestFilterBit(fFilterBit))
400 //        continue;
401         
402         CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
403         
404         if (mcP) {
405           // Fill the hidden information with the simulated data
406           Int_t pLabel = aodtrack->GetLabel();
407           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
408           
409           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
410           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
411           if (!tPart) {
412             fpx = fV1[0];
413             fpy = fV1[1];
414             fpz = fV1[2];
415             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
416             tInfo->SetPDGPid(0);
417             tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
418             tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
419             tInfo->SetMass(0);
420           }
421           else {
422             // Check the mother information
423           
424             // Using the new way of storing the freeze-out information
425             // Final state particle is stored twice on the stack
426             // one copy (mother) is stored with original freeze-out information
427             //   and is not tracked
428             // the other one (daughter) is stored with primary vertex position
429             //   and is tracked
430             
431             // Freeze-out coordinates
432             fpx = tPart->Xv() - fV1[0];
433             fpy = tPart->Yv() - fV1[1];
434             fpz = tPart->Zv() - fV1[2];
435             //    fpt = tPart->T();
436             
437             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
438             
439             fpx *= 1e13;
440             fpy *= 1e13;
441             fpz *= 1e13;
442             //    fpt *= 1e13;
443             
444             //      cout << "Looking for mother ids " << endl;
445             if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
446               //        cout << "Got mother id" << endl;
447               AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
448               // Check if this is the same particle stored twice on the stack
449               if (mother) {
450                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
451                   // It is the same particle
452                   // Read in the original freeze-out information
453                   // and convert it from to [fm]
454                   
455                   // EPOS style 
456                   //      fpx = mother->Xv()*1e13*0.197327;
457                   //      fpy = mother->Yv()*1e13*0.197327;
458                   //      fpz = mother->Zv()*1e13*0.197327;
459                   //      fpt = mother->T() *1e13*0.197327*0.5;
460                   
461                   
462                   // Therminator style 
463                   fpx = mother->Xv()*1e13;
464                   fpy = mother->Yv()*1e13;
465                   fpz = mother->Zv()*1e13;
466                   //          fpt = mother->T() *1e13*3e10;
467                   
468                 }
469               }
470             }
471             
472             //       if (fRotateToEventPlane) {
473             //  double tPhi = TMath::ATan2(fpy, fpx);
474             //  double tRad = TMath::Hypot(fpx, fpy);
475             
476             //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
477             //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
478             //       }
479             
480             tInfo->SetPDGPid(tPart->GetPdgCode());
481             
482             //    if (fRotateToEventPlane) {
483             //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
484             //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
485             
486             //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
487             //                             tRad*TMath::Sin(tPhi - tReactionPlane),
488             //                             tPart->Pz());
489             //    }
490             //       else
491             tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
492             Double_t mass2 = (tPart->E() *tPart->E() -
493                               tPart->Px()*tPart->Px() -
494                               tPart->Py()*tPart->Py() -
495                               tPart->Pz()*tPart->Pz());
496             if (mass2>0.0)
497               tInfo->SetMass(TMath::Sqrt(mass2));
498             else 
499               tInfo->SetMass(0.0);
500             
501             tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
502           }
503           trackCopy->SetHiddenInfo(tInfo);
504         }
505
506         double pxyz[3];
507         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
508         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
509         // Check the sanity of the tracks - reject zero momentum tracks
510         if (ktP.mag() == 0) {
511           delete trackCopy;
512           continue;
513         }
514       }
515
516
517       tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
518       realnofTracks++;//real number of tracks           
519     }
520
521   tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event     
522
523   if (mcP) delete motherids;
524
525   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
526 }
527
528 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack, 
529                                                  AliFemtoTrack *tFemtoTrack, 
530                                                  AliPWG2AODTrack *tPWG2AODTrack)
531 {
532   // Copy the track information from the AOD into the internal AliFemtoTrack
533   // If it exists, use the additional information from the PWG2 AOD
534
535   // Primary Vertex position
536   double fV1[3];
537   fEvent->GetPrimaryVertex()->GetPosition(fV1);
538
539   tFemtoTrack->SetCharge(tAodTrack->Charge());
540   
541   //in aliroot we have AliPID 
542   //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
543   //we use only 5 first
544
545   // AOD pid has 10 components
546   double aodpid[10];
547   tAodTrack->GetPID(aodpid);
548   tFemtoTrack->SetPidProbElectron(aodpid[0]);
549   tFemtoTrack->SetPidProbMuon(aodpid[1]);
550   tFemtoTrack->SetPidProbPion(aodpid[2]);
551   tFemtoTrack->SetPidProbKaon(aodpid[3]);
552   tFemtoTrack->SetPidProbProton(aodpid[4]);
553                                                 
554   double pxyz[3];
555   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
556   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
557   tFemtoTrack->SetP(v);//setting momentum
558   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
559   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
560   //setting track helix 
561   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
562   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge())); 
563   tFemtoTrack->SetHelix(helix);
564                 
565   // Flags
566   tFemtoTrack->SetTrackId(tAodTrack->GetID());
567   tFemtoTrack->SetFlags(1);
568   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
569                 
570   // Track quality information 
571   float covmat[6];
572   tAodTrack->GetCovMatrix(covmat);
573   tFemtoTrack->SetImpactD(covmat[0]);
574   tFemtoTrack->SetImpactZ(covmat[2]);
575   tFemtoTrack->SetCdd(covmat[3]);
576   tFemtoTrack->SetCdz(covmat[4]);
577   tFemtoTrack->SetCzz(covmat[5]);
578   // This information is only available in the ESD
579   // We put in fake values or reasonable estimates
580   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());    
581   tFemtoTrack->SetITSncls(1);     
582   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());       
583   tFemtoTrack->SetTPCncls(1);       
584   tFemtoTrack->SetTPCnclsF(1);      
585   tFemtoTrack->SetTPCsignalN(1); 
586   tFemtoTrack->SetTPCsignalS(1); 
587
588   if (tPWG2AODTrack) {
589     // Copy the PWG2 specific information if it exists
590     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
591     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
592     
593     double xtpc[3] = {0,0,0};
594     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
595     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
596     tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
597     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
598   }
599   else {
600     // If not use dummy values
601     tFemtoTrack->SetTPCClusterMap(fAllTrue);
602     tFemtoTrack->SetTPCSharedMap(fAllFalse);
603     
604     double xtpc[3] = {0,0,0};
605     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
606     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
607   }
608
609   int indexes[3];
610   for (int ik=0; ik<3; ik++) {
611     indexes[ik] = 0;
612   }
613   tFemtoTrack->SetKinkIndexes(indexes);
614 }
615
616 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
617 {
618   fFilterBit = (1 << (ibit));
619 }
620
621 void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
622 {
623   fReadMC = a;
624 }
625
626 AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
627 {
628   if (aLabel < 0) return 0;
629   AliAODMCParticle *aodP;
630   Int_t posstack = 0;
631   if (aLabel > mcP->GetEntries())
632     posstack = mcP->GetEntries();
633   else
634     posstack = aLabel;
635
636   aodP = (AliAODMCParticle *) mcP->At(posstack);
637   if (aodP->GetLabel() > posstack) {
638     do {
639       aodP = (AliAODMCParticle *) mcP->At(posstack);
640       if (aodP->GetLabel() == aLabel) return aodP;
641       posstack--;
642     }
643     while (posstack > 0);
644   }
645   else {
646     do {
647       aodP = (AliAODMCParticle *) mcP->At(posstack);
648       if (aodP->GetLabel() == aLabel) return aodP;
649       posstack++;
650     }
651     while (posstack < mcP->GetEntries());
652   }
653   
654   return 0;
655 }
656
657
658
659
660