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   int tracksPrim=0;
 278
 279   for (int i=0;i<nofTracks;i++)
 280     {
 281       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 282
 283       if (fPWG2AODTracks) {
 284         // Read tracks from the additional pwg2 specific AOD part
 285         // if they exist
 286         // Note that in that case all the AOD tracks without the
 287         // additional information will be ignored !
 288         AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
 289
 290         // Getting the AOD track through the ref of the additional info
 291         AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
 292         if (!aodtrack->TestFilterBit(fFilterBit)) {
 293           delete trackCopy;
 294           continue;
 295         }
 296
 297         CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
 298
 299         if (mcP) {
 300           // Fill the hidden information with the simulated data
 301           //      Int_t pLabel = aodtrack->GetLabel();
 302           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 303
 304           // Check the mother information
 305
 306           // Using the new way of storing the freeze-out information
 307           // Final state particle is stored twice on the stack
 308           // one copy (mother) is stored with original freeze-out information
 309           //   and is not tracked
 310           // the other one (daughter) is stored with primary vertex position
 311           //   and is tracked
 312
 313           // Freeze-out coordinates
 314           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 315           fpx = tPart->Xv() - fV1[0];
 316           fpy = tPart->Yv() - fV1[1];
 317           fpz = tPart->Zv() - fV1[2];
 318           fpt = tPart->T();
 319
 320           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 321           tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 322
 323           fpx *= 1e13;
 324           fpy *= 1e13;
 325           fpz *= 1e13;
 326           fpt *= 1e13;
 327
 328           //      cout << "Looking for mother ids " << endl;
 329           if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 330             //  cout << "Got mother id" << endl;
 331             AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 332             // Check if this is the same particle stored twice on the stack
 333             if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 334               // It is the same particle
 335               // Read in the original freeze-out information
 336               // and convert it from to [fm]
 337
 338               // EPOS style
 339               //          fpx = mother->Xv()*1e13*0.197327;
 340               //          fpy = mother->Yv()*1e13*0.197327;
 341               //          fpz = mother->Zv()*1e13*0.197327;
 342               //          fpt = mother->T() *1e13*0.197327*0.5;
 343
 344
 345               // Therminator style
 346               fpx = mother->Xv()*1e13;
 347               fpy = mother->Yv()*1e13;
 348               fpz = mother->Zv()*1e13;
 349               fpt = mother->T() *1e13*3e10;
 350
 351             }
 352           }
 353
 354 //       if (fRotateToEventPlane) {
 355 //      double tPhi = TMath::ATan2(fpy, fpx);
 356 //      double tRad = TMath::Hypot(fpx, fpy);
 357
 358 //      fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 359 //      fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 360 //       }
 361
 362           tInfo->SetPDGPid(tPart->GetPdgCode());
 363
 364 //        if (fRotateToEventPlane) {
 365 //          double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 366 //          double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 367
 368 //          tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 369 //                                 tRad*TMath::Sin(tPhi - tReactionPlane),
 370 //                                 tPart->Pz());
 371 //        }
 372 //       else
 373           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 374           Double_t mass2 = (tPart->E() *tPart->E() -
 375                             tPart->Px()*tPart->Px() -
 376                             tPart->Py()*tPart->Py() -
 377                             tPart->Pz()*tPart->Pz());
 378           if (mass2>0.0)
 379             tInfo->SetMass(TMath::Sqrt(mass2));
 380           else
 381             tInfo->SetMass(0.0);
 382
 383           tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 384           trackCopy->SetHiddenInfo(tInfo);
 385
 386         }
 387
 388         double pxyz[3];
 389         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 390         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 391         // Check the sanity of the tracks - reject zero momentum tracks
 392         if (ktP.Mag() == 0) {
 393           delete trackCopy;
 394           continue;
 395         }
 396       }
 397       else {
 398         // No additional information exists
 399         // Read in the normal AliAODTracks
 400         const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 401
 402         if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
 403
 404 //      if (!aodtrack->TestFilterBit(fFilterBit))
 405 //        continue;
 406
 407         CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
 408
 409         if (mcP) {
 410           // Fill the hidden information with the simulated data
 411           //      Int_t pLabel = aodtrack->GetLabel();
 412           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 413
 414           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 415           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 416           if (!tPart) {
 417             fpx = fV1[0];
 418             fpy = fV1[1];
 419             fpz = fV1[2];
 420             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 421             tInfo->SetPDGPid(0);
 422             tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
 423             tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
 424             tInfo->SetMass(0);
 425           }
 426           else {
 427             // Check the mother information
 428
 429             // Using the new way of storing the freeze-out information
 430             // Final state particle is stored twice on the stack
 431             // one copy (mother) is stored with original freeze-out information
 432             //   and is not tracked
 433             // the other one (daughter) is stored with primary vertex position
 434             //   and is tracked
 435
 436             // Freeze-out coordinates
 437             fpx = tPart->Xv() - fV1[0];
 438             fpy = tPart->Yv() - fV1[1];
 439             fpz = tPart->Zv() - fV1[2];
 440             //    fpt = tPart->T();
 441
 442             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 443
 444             fpx *= 1e13;
 445             fpy *= 1e13;
 446             fpz *= 1e13;
 447             //    fpt *= 1e13;
 448
 449             //      cout << "Looking for mother ids " << endl;
 450             if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 451               //        cout << "Got mother id" << endl;
 452               AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 453               // Check if this is the same particle stored twice on the stack
 454               if (mother) {
 455                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 456                   // It is the same particle
 457                   // Read in the original freeze-out information
 458                   // and convert it from to [fm]
 459
 460                   // EPOS style
 461                   //      fpx = mother->Xv()*1e13*0.197327;
 462                   //      fpy = mother->Yv()*1e13*0.197327;
 463                   //      fpz = mother->Zv()*1e13*0.197327;
 464                   //      fpt = mother->T() *1e13*0.197327*0.5;
 465
 466
 467                   // Therminator style
 468                   fpx = mother->Xv()*1e13;
 469                   fpy = mother->Yv()*1e13;
 470                   fpz = mother->Zv()*1e13;
 471                   //          fpt = mother->T() *1e13*3e10;
 472
 473                 }
 474               }
 475             }
 476
 477             //       if (fRotateToEventPlane) {
 478             //  double tPhi = TMath::ATan2(fpy, fpx);
 479             //  double tRad = TMath::Hypot(fpx, fpy);
 480
 481             //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 482             //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 483             //       }
 484
 485             tInfo->SetPDGPid(tPart->GetPdgCode());
 486
 487             //    if (fRotateToEventPlane) {
 488             //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 489             //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 490
 491             //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 492             //                             tRad*TMath::Sin(tPhi - tReactionPlane),
 493             //                             tPart->Pz());
 494             //    }
 495             //       else
 496             tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 497             Double_t mass2 = (tPart->E() *tPart->E() -
 498                               tPart->Px()*tPart->Px() -
 499                               tPart->Py()*tPart->Py() -
 500                               tPart->Pz()*tPart->Pz());
 501             if (mass2>0.0)
 502               tInfo->SetMass(TMath::Sqrt(mass2));
 503             else
 504               tInfo->SetMass(0.0);
 505
 506             tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 507           }
 508           trackCopy->SetHiddenInfo(tInfo);
 509         }
 510
 511         double pxyz[3];
 512         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 513         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 514         // Check the sanity of the tracks - reject zero momentum tracks
 515         if (ktP.Mag() == 0) {
 516           delete trackCopy;
 517           continue;
 518         }
 519       }
 520
 521
 522       tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 523       realnofTracks++;//real number of tracks
 524     }
 525
 526   tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 527   tEvent->SetNormalizedMult(tracksPrim);
 528
 529   AliCentrality *cent = fEvent->GetCentrality();
 530   if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
 531
 532   if (mcP) delete [] motherids;
 533
 534   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 535 }
 536
 537 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack,
 538                                                  AliFemtoTrack *tFemtoTrack,
 539                                                  AliPWG2AODTrack *tPWG2AODTrack)
 540 {
 541   // Copy the track information from the AOD into the internal AliFemtoTrack
 542   // If it exists, use the additional information from the PWG2 AOD
 543
 544   // Primary Vertex position
 545   double fV1[3];
 546   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 547
 548   tFemtoTrack->SetCharge(tAodTrack->Charge());
 549
 550   //in aliroot we have AliPID
 551   //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 552   //we use only 5 first
 553
 554   // AOD pid has 10 components
 555   double aodpid[10];
 556   tAodTrack->GetPID(aodpid);
 557   tFemtoTrack->SetPidProbElectron(aodpid[0]);
 558   tFemtoTrack->SetPidProbMuon(aodpid[1]);
 559   tFemtoTrack->SetPidProbPion(aodpid[2]);
 560   tFemtoTrack->SetPidProbKaon(aodpid[3]);
 561   tFemtoTrack->SetPidProbProton(aodpid[4]);
 562
 563   double pxyz[3];
 564   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
 565   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 566   tFemtoTrack->SetP(v);//setting momentum
 567   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 568   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 569   //setting track helix
 570   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 571   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
 572   tFemtoTrack->SetHelix(helix);
 573
 574   // Flags
 575   tFemtoTrack->SetTrackId(tAodTrack->GetID());
 576   tFemtoTrack->SetFlags(tAodTrack->GetFlags());
 577   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
 578
 579   // Track quality information
 580   float covmat[6];
 581   tAodTrack->GetCovMatrix(covmat);
 582
 583 //   if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
 584 //     tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
 585 //   else
 586 //     tFemtoTrack->SetImpactD(0.0);
 587 //   tFemtoTrack->SetImpactD(tAodTrack->DCA());
 588
 589 //   tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
 590   tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
 591   tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
 592
 593   //  cout << fV1[0] << " " << tAodTrack->XAtDCA() << " " << tAodTrack->Xv() << endl;
 594
 595   tFemtoTrack->SetCdd(covmat[0]);
 596   tFemtoTrack->SetCdz(covmat[1]);
 597   tFemtoTrack->SetCzz(covmat[2]);
 598   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
 599   tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
 600   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 601   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
 602   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
 603   tFemtoTrack->SetTPCsignalN(1);
 604   tFemtoTrack->SetTPCsignalS(1);
 605
 606   if (tPWG2AODTrack) {
 607     // Copy the PWG2 specific information if it exists
 608     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
 609     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
 610
 611     double xtpc[3] = {0,0,0};
 612     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
 613     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 614     tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
 615     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 616   }
 617   else {
 618     // If not use dummy values
 619     tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
 620     tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
 621
 622     double xtpc[3] = {0,0,0};
 623     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 624     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 625   }
 626
 627   //  cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << "  " << tAodTrack->Zv() << "  " << tAodTrack->GetTPCNcls() << endl;
 628
 629
 630   int indexes[3];
 631   for (int ik=0; ik<3; ik++) {
 632     indexes[ik] = 0;
 633   }
 634   tFemtoTrack->SetKinkIndexes(indexes);
 635 }
 636
 637 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
 638 {
 639   fFilterBit = (1 << (ibit));
 640 }
 641
 642 void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
 643 {
 644   fReadMC = a;
 645 }
 646
 647 AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
 648 {
 649   if (aLabel < 0) return 0;
 650   AliAODMCParticle *aodP;
 651   Int_t posstack = 0;
 652   if (aLabel > mcP->GetEntries())
 653     posstack = mcP->GetEntries();
 654   else
 655     posstack = aLabel;
 656
 657   aodP = (AliAODMCParticle *) mcP->At(posstack);
 658   if (aodP->GetLabel() > posstack) {
 659     do {
 660       aodP = (AliAODMCParticle *) mcP->At(posstack);
 661       if (aodP->GetLabel() == aLabel) return aodP;
 662       posstack--;
 663     }
 664     while (posstack > 0);
 665   }
 666   else {
 667     do {
 668       aodP = (AliAODMCParticle *) mcP->At(posstack);
 669       if (aodP->GetLabel() == aLabel) return aodP;
 670       posstack++;
 671     }
 672     while (posstack < mcP->GetEntries());
 673   }
 674
 675   return 0;
 676 }
 677
 678
 679
 680
 681