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=0;
 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=0;
 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