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   //  fPWG2AODTracks(0x0),
  52   fReadMC(0),
  53   fReadV0(0),
  54   fUsePreCent(0),
  55   fNoCentrality(0),
  56   fAODpidUtil(0),
  57   fInputFile(" "),
  58   fFileName(" "),
  59   fTree(0x0),
  60   fAodFile(0x0)
  61 {
  62   // default constructor
  63   fAllTrue.ResetAllBits(kTRUE);
  64   fAllFalse.ResetAllBits(kFALSE);
  65   fCentRange[0] = 0;
  66   fCentRange[1] = 1000;
  67 }
  68
  69 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
  70   AliFemtoEventReader(),
  71   fNumberofEvent(0),
  72   fCurEvent(0),
  73   fEvent(0x0),
  74   fAllTrue(160),
  75   fAllFalse(160),
  76   fFilterBit(0),
  77   //  fPWG2AODTracks(0x0),
  78   fReadMC(0),
  79   fReadV0(0),
  80   fUsePreCent(0),
  81   fNoCentrality(0),
  82   fAODpidUtil(0),
  83   fInputFile(" "),
  84   fFileName(" "),
  85   fTree(0x0),
  86   fAodFile(0x0)
  87 {
  88   // copy constructor
  89   fReadMC = aReader.fReadMC;
  90   fReadV0 = aReader.fReadV0;
  91   fInputFile = aReader.fInputFile;
  92   fFileName  = aReader.fFileName;
  93   fNumberofEvent = aReader.fNumberofEvent;
  94   fCurEvent = aReader.fCurEvent;
  95   fEvent = new AliAODEvent();
  96   fAodFile = new TFile(aReader.fAodFile->GetName());
  97   fAllTrue.ResetAllBits(kTRUE);
  98   fAllFalse.ResetAllBits(kFALSE);
  99   fFilterBit = aReader.fFilterBit;
 100   //  fPWG2AODTracks = aReader.fPWG2AODTracks;
 101   fAODpidUtil = aReader.fAODpidUtil;
 102   fCentRange[0] = aReader.fCentRange[0];
 103   fCentRange[1] = aReader.fCentRange[1];
 104   fNoCentrality = aReader.fNoCentrality;
 105   fUsePreCent = aReader.fUsePreCent;
 106 }
 107 //__________________
 108 //Destructor
 109 AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
 110 {
 111   // destructor
 112   delete fTree;
 113   delete fEvent;
 114   delete fAodFile;
 115 //   if (fPWG2AODTracks) {
 116 //     fPWG2AODTracks->Delete();
 117 //     delete fPWG2AODTracks;
 118 //   }
 119 }
 120
 121 //__________________
 122 AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
 123 {
 124   // assignment operator
 125   if (this == &aReader)
 126     return *this;
 127
 128   fInputFile = aReader.fInputFile;
 129   fFileName  = aReader.fFileName;
 130   fNumberofEvent = aReader.fNumberofEvent;
 131   fCurEvent = aReader.fCurEvent;
 132   if (fTree) delete fTree;
 133   if (fEvent) delete fEvent;
 134   fEvent = new AliAODEvent();
 135   if (fAodFile) delete fAodFile;
 136   fAodFile = new TFile(aReader.fAodFile->GetName());
 137   fAllTrue.ResetAllBits(kTRUE);
 138   fAllFalse.ResetAllBits(kFALSE);
 139   fFilterBit = aReader.fFilterBit;
 140   //  fPWG2AODTracks = aReader.fPWG2AODTracks;
 141   fAODpidUtil = aReader.fAODpidUtil;
 142   fCentRange[0] = aReader.fCentRange[0];
 143   fCentRange[1] = aReader.fCentRange[1];
 144   fUsePreCent = aReader.fUsePreCent;
 145   fNoCentrality = aReader.fNoCentrality;
 146
 147   return *this;
 148 }
 149 //__________________
 150 AliFemtoString AliFemtoEventReaderAOD::Report()
 151 {
 152   // create reader report
 153   AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
 154   return temp;
 155 }
 156
 157 //__________________
 158 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
 159 {
 160   //setting the name of file where names of AOD file are written
 161   //it takes only this files which have good trees
 162   char buffer[256];
 163   fInputFile=string(inputFile);
 164   ifstream infile(inputFile);
 165
 166   fTree = new TChain("aodTree");
 167
 168   if(infile.good()==true)
 169     {
 170       //checking if all give files have good tree inside
 171       while (infile.eof()==false)
 172         {
 173           infile.getline(buffer,256);
 174           TFile *aodFile=TFile::Open(buffer,"READ");
 175           if (aodFile!=0x0)
 176             {
 177               TTree* tree = (TTree*) aodFile->Get("aodTree");
 178               if (tree!=0x0)
 179                 {
 180                   //              cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
 181                   fTree->AddFile(buffer);
 182                   delete tree;
 183                 }
 184               aodFile->Close();
 185             }
 186           delete aodFile;
 187         }
 188     }
 189 }
 190
 191 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
 192 {
 193   // read in a next hbt event from the chain
 194   // convert it to AliFemtoEvent and return
 195   // for further analysis
 196   AliFemtoEvent *hbtEvent = 0;
 197   cout<<"reader"<<endl;
 198   if (fCurEvent==fNumberofEvent)//open next file
 199     {
 200       if(fNumberofEvent==0)
 201         {
 202           fEvent=new AliAODEvent();
 203           fEvent->ReadFromTree(fTree);
 204
 205           // Check for the existence of the additional information
 206 //        fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
 207
 208 //        if (fPWG2AODTracks) {
 209 //          cout << "Found additional PWG2 specific information in the AOD!" << endl;
 210 //          cout << "Reading only tracks with the additional information" << endl;
 211 //        }
 212
 213           fNumberofEvent=fTree->GetEntries();
 214           //      cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
 215           fCurEvent=0;
 216         }
 217       else //no more data to read
 218         {
 219           cout<<"no more files "<<hbtEvent<<endl;
 220           fReaderStatus=1;
 221           return hbtEvent;
 222         }
 223     }
 224
 225   cout<<"starting to read event "<<fCurEvent<<endl;
 226   fTree->GetEvent(fCurEvent);//getting next event
 227   //  cout << "Read event " << fEvent << " from file " << fTree << endl;
 228
 229   hbtEvent = new AliFemtoEvent;
 230
 231   CopyAODtoFemtoEvent(hbtEvent);
 232   fCurEvent++;
 233
 234
 235   return hbtEvent;
 236 }
 237
 238 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 239 {
 240
 241   // A function that reads in the AOD event
 242   // and transfers the neccessary information into
 243   // the internal AliFemtoEvent
 244
 245   // setting global event characteristics
 246   tEvent->SetRunNumber(fEvent->GetRunNumber());
 247   tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 248   tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 249   tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 250   tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 251   tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 252   tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
 253   tEvent->SetZDCParticipants(0);
 254   tEvent->SetTriggerMask(fEvent->GetTriggerMask());
 255   tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 256
 257   // Attempt to access MC header
 258   AliAODMCHeader *mcH;
 259   TClonesArray *mcP=0;
 260   if (fReadMC) {
 261     mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
 262     if (!mcH) {
 263       cout << "AOD MC information requested, but no header found!" << endl;
 264     }
 265
 266     mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
 267     if (!mcP) {
 268       cout << "AOD MC information requested, but no particle array found!" << endl;
 269     }
 270   }
 271
 272   tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
 273
 274   Int_t *motherids=0;
 275   if (mcP) {
 276     motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
 277     for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
 278
 279     // Read in mother ids
 280     AliAODMCParticle *motherpart;
 281     for (int ip=0; ip<mcP->GetEntries(); ip++) {
 282       motherpart = (AliAODMCParticle *) mcP->At(ip);
 283       if (motherpart->GetDaughter(0) > 0)
 284         motherids[motherpart->GetDaughter(0)] = ip;
 285       if (motherpart->GetDaughter(1) > 0)
 286         motherids[motherpart->GetDaughter(1)] = ip;
 287     }
 288   }
 289
 290   // Primary Vertex position
 291   //  double fV1[3];
 292   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 293
 294   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 295   tEvent->SetPrimVertPos(vertex);
 296
 297   //starting to reading tracks
 298   int nofTracks=0;  //number of reconstructed tracks in event
 299
 300   // Check to see whether the additional info exists
 301 //   if (fPWG2AODTracks)
 302 //     nofTracks=fPWG2AODTracks->GetEntries();
 303 //   else
 304     nofTracks=fEvent->GetNumberOfTracks();
 305   cout<<"nofTracks: "<<nofTracks<<endl;
 306
 307   if (!fNoCentrality)
 308     AliCentrality *cent = fEvent->GetCentrality();
 309
 310   if (!fNoCentrality && cent && fUsePreCent) {
 311     if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
 312         (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
 313       {
 314         cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
 315
 316         return;
 317       }
 318   }
 319
 320   int realnofTracks=0;   // number of track which we use in a analysis
 321   int tracksPrim=0;
 322
 323   int labels[20000];
 324   for (int il=0; il<20000; il++) labels[il] = -1;
 325
 326   // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
 327   for (int i=0;i<nofTracks;i++) {
 328     const AliAODTrack *aodtrack=fEvent->GetTrack(i);
 329     if (!aodtrack->TestFilterBit(fFilterBit)) {
 330       if(aodtrack->GetID() < 0) continue;
 331       labels[aodtrack->GetID()] = i;
 332     }
 333   }
 334
 335   //  int tNormMult = 0;
 336   for (int i=0;i<nofTracks;i++)
 337     {
 338       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 339
 340 //       if (fPWG2AODTracks) {
 341 //      // Read tracks from the additional pwg2 specific AOD part
 342 //      // if they exist
 343 //      // Note that in that case all the AOD tracks without the
 344 //      // additional information will be ignored !
 345 //      AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
 346
 347 //      // Getting the AOD track through the ref of the additional info
 348 //      AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
 349 //      if (!aodtrack->TestFilterBit(fFilterBit)) {
 350 //        delete trackCopy;
 351 //        continue;
 352 //      }
 353
 354
 355 //      if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
 356 //        if (aodtrack->Chi2perNDF() < 6.0)
 357 //          if (aodtrack->Eta() < 0.9)
 358 //            tNormMult++;
 359
 360
 361 //      CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
 362
 363 //      if (mcP) {
 364 //        // Fill the hidden information with the simulated data
 365 //        //      Int_t pLabel = aodtrack->GetLabel();
 366 //        AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 367
 368 //        // Check the mother information
 369
 370 //        // Using the new way of storing the freeze-out information
 371 //        // Final state particle is stored twice on the stack
 372 //        // one copy (mother) is stored with original freeze-out information
 373 //        //   and is not tracked
 374 //        // the other one (daughter) is stored with primary vertex position
 375 //        //   and is tracked
 376
 377 //        // Freeze-out coordinates
 378 //        double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 379 //        fpx = tPart->Xv() - fV1[0];
 380 //        fpy = tPart->Yv() - fV1[1];
 381 //        fpz = tPart->Zv() - fV1[2];
 382 //        fpt = tPart->T();
 383
 384 //        AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 385 //        tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 386
 387 //        fpx *= 1e13;
 388 //        fpy *= 1e13;
 389 //        fpz *= 1e13;
 390 //        fpt *= 1e13;
 391
 392 //        //      cout << "Looking for mother ids " << endl;
 393 //        if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 394 //          //  cout << "Got mother id" << endl;
 395 //          AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 396 //          // Check if this is the same particle stored twice on the stack
 397 //          if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 398 //            // It is the same particle
 399 //            // Read in the original freeze-out information
 400 //            // and convert it from to [fm]
 401
 402 //            // EPOS style
 403 //            //          fpx = mother->Xv()*1e13*0.197327;
 404 //            //          fpy = mother->Yv()*1e13*0.197327;
 405 //            //          fpz = mother->Zv()*1e13*0.197327;
 406 //            //          fpt = mother->T() *1e13*0.197327*0.5;
 407
 408
 409 //            // Therminator style
 410 //            fpx = mother->Xv()*1e13;
 411 //            fpy = mother->Yv()*1e13;
 412 //            fpz = mother->Zv()*1e13;
 413 //            fpt = mother->T() *1e13*3e10;
 414
 415 //          }
 416 //        }
 417
 418 //        //       if (fRotateToEventPlane) {
 419 //        //    double tPhi = TMath::ATan2(fpy, fpx);
 420 //        //    double tRad = TMath::Hypot(fpx, fpy);
 421
 422 //        //    fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 423 //        //    fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 424 //        //       }
 425
 426 //        tInfo->SetPDGPid(tPart->GetPdgCode());
 427
 428 //        //      if (fRotateToEventPlane) {
 429 //        //        double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 430 //        //        double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 431
 432 //        //        tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 433 //        //                               tRad*TMath::Sin(tPhi - tReactionPlane),
 434 //        //                               tPart->Pz());
 435 //        //      }
 436 //        //       else
 437 //        tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 438 //        Double_t mass2 = (tPart->E() *tPart->E() -
 439 //                          tPart->Px()*tPart->Px() -
 440 //                          tPart->Py()*tPart->Py() -
 441 //                          tPart->Pz()*tPart->Pz());
 442 //        if (mass2>0.0)
 443 //          tInfo->SetMass(TMath::Sqrt(mass2));
 444 //        else
 445 //          tInfo->SetMass(0.0);
 446
 447 //        tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 448 //        trackCopy->SetHiddenInfo(tInfo);
 449
 450 //      }
 451
 452 //      double pxyz[3];
 453 //      aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 454 //      const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 455 //      // Check the sanity of the tracks - reject zero momentum tracks
 456 //      if (ktP.Mag() == 0) {
 457 //        delete trackCopy;
 458 //        continue;
 459 //      }
 460 //       }
 461 //       else {
 462         // No additional information exists
 463         // Read in the normal AliAODTracks
 464
 465         //      const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 466         AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 467
 468         if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
 469
 470         if (!aodtrack->TestFilterBit(fFilterBit)) {
 471           delete trackCopy;
 472           continue;
 473         }
 474
 475         CopyAODtoFemtoTrack(aodtrack, trackCopy);
 476
 477         // copying PID information from the correspondent track
 478         //      const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
 479         AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
 480         CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
 481
 482         if (mcP) {
 483           // Fill the hidden information with the simulated data
 484           //      Int_t pLabel = aodtrack->GetLabel();
 485           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 486
 487           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 488           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 489           if (!tPart) {
 490             fpx = fV1[0];
 491             fpy = fV1[1];
 492             fpz = fV1[2];
 493             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 494             tInfo->SetPDGPid(0);
 495             tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
 496             tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
 497             tInfo->SetMass(0);
 498           }
 499           else {
 500             // Check the mother information
 501
 502             // Using the new way of storing the freeze-out information
 503             // Final state particle is stored twice on the stack
 504             // one copy (mother) is stored with original freeze-out information
 505             //   and is not tracked
 506             // the other one (daughter) is stored with primary vertex position
 507             //   and is tracked
 508
 509             // Freeze-out coordinates
 510             fpx = tPart->Xv() - fV1[0];
 511             fpy = tPart->Yv() - fV1[1];
 512             fpz = tPart->Zv() - fV1[2];
 513             //    fpt = tPart->T();
 514
 515             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 516
 517             fpx *= 1e13;
 518             fpy *= 1e13;
 519             fpz *= 1e13;
 520             //    fpt *= 1e13;
 521
 522             //      cout << "Looking for mother ids " << endl;
 523             if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 524               //        cout << "Got mother id" << endl;
 525               AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 526               // Check if this is the same particle stored twice on the stack
 527               if (mother) {
 528                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 529                   // It is the same particle
 530                   // Read in the original freeze-out information
 531                   // and convert it from to [fm]
 532
 533                   // EPOS style
 534                   //      fpx = mother->Xv()*1e13*0.197327;
 535                   //      fpy = mother->Yv()*1e13*0.197327;
 536                   //      fpz = mother->Zv()*1e13*0.197327;
 537                   //      fpt = mother->T() *1e13*0.197327*0.5;
 538
 539
 540                   // Therminator style
 541                   fpx = mother->Xv()*1e13;
 542                   fpy = mother->Yv()*1e13;
 543                   fpz = mother->Zv()*1e13;
 544                   //          fpt = mother->T() *1e13*3e10;
 545
 546                 }
 547               }
 548             }
 549
 550             //       if (fRotateToEventPlane) {
 551             //  double tPhi = TMath::ATan2(fpy, fpx);
 552             //  double tRad = TMath::Hypot(fpx, fpy);
 553
 554             //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 555             //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 556             //       }
 557
 558             tInfo->SetPDGPid(tPart->GetPdgCode());
 559
 560             //    if (fRotateToEventPlane) {
 561             //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 562             //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 563
 564             //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 565             //                             tRad*TMath::Sin(tPhi - tReactionPlane),
 566             //                             tPart->Pz());
 567             //    }
 568             //       else
 569             tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 570             Double_t mass2 = (tPart->E() *tPart->E() -
 571                               tPart->Px()*tPart->Px() -
 572                               tPart->Py()*tPart->Py() -
 573                               tPart->Pz()*tPart->Pz());
 574             if (mass2>0.0)
 575               tInfo->SetMass(TMath::Sqrt(mass2));
 576             else
 577               tInfo->SetMass(0.0);
 578
 579             tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 580           }
 581           trackCopy->SetHiddenInfo(tInfo);
 582         }
 583
 584         double pxyz[3];
 585         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 586         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 587         // Check the sanity of the tracks - reject zero momentum tracks
 588         if (ktP.Mag() == 0) {
 589           delete trackCopy;
 590           continue;
 591         }
 592         //    }
 593
 594
 595         tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 596         realnofTracks++;//real number of tracks
 597     }
 598
 599   tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 600   tEvent->SetNormalizedMult(tracksPrim);
 601
 602   if (!fNoCentrality) {
 603     //  AliCentrality *cent = fEvent->GetCentrality();
 604     if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
 605     //  if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
 606
 607     if (cent) {
 608       tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 609       //    tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 610       tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 611       //    tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 612     }
 613   }
 614
 615   if (mcP) delete [] motherids;
 616
 617   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 618
 619   if(fReadV0)
 620     {
 621
 622       for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
 623         AliAODv0* aodv0 = fEvent->GetV0(i);
 624         if (!aodv0) continue;
 625         if(aodv0->GetNDaughters()>2) continue;
 626         if(aodv0->GetNProngs()>2) continue;
 627         if(aodv0->GetCharge()!=0) continue;
 628         if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
 629         AliFemtoV0* trackCopyV0 = new AliFemtoV0();
 630         CopyAODtoFemtoV0(aodv0, trackCopyV0);
 631         tEvent->V0Collection()->push_back(trackCopyV0);
 632         //cout<<"Pushback v0 to v0collection"<<endl;
 633       }
 634     }
 635
 636 }
 637
 638 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
 639                                                  AliFemtoTrack *tFemtoTrack
 640                                                  //                                              AliPWG2AODTrack *tPWG2AODTrack
 641                                                  )
 642 {
 643   // Copy the track information from the AOD into the internal AliFemtoTrack
 644   // If it exists, use the additional information from the PWG2 AOD
 645
 646   // Primary Vertex position
 647
 648   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 649   //  fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 650
 651   tFemtoTrack->SetCharge(tAodTrack->Charge());
 652
 653   double pxyz[3];
 654   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
 655   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 656   tFemtoTrack->SetP(v);//setting momentum
 657   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 658   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 659   //setting track helix
 660   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 661   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
 662   tFemtoTrack->SetHelix(helix);
 663
 664   // Flags
 665   tFemtoTrack->SetTrackId(tAodTrack->GetID());
 666   tFemtoTrack->SetFlags(tAodTrack->GetFlags());
 667   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
 668
 669   // Track quality information
 670   float covmat[6];
 671   tAodTrack->GetCovMatrix(covmat);
 672
 673   double impact[2];
 674   double covimpact[3];
 675
 676   if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
 677     //cout << "sth went wrong with dca propagation" << endl;
 678     tFemtoTrack->SetImpactD(-1000.0);
 679     tFemtoTrack->SetImpactZ(-1000.0);
 680
 681   }
 682   else {
 683     tFemtoTrack->SetImpactD(impact[0]);
 684     tFemtoTrack->SetImpactZ(impact[1]+fV1[2]);
 685   }
 686
 687   //   if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
 688   //     tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
 689   //   else
 690   //     tFemtoTrack->SetImpactD(0.0);
 691   //   tFemtoTrack->SetImpactD(tAodTrack->DCA());
 692
 693   //   tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
 694
 695
 696   //   tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
 697   //   tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
 698
 699
 700   //   cout
 701     //    << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
 702     //    << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
 703     //    << tAodTrack->Yv() - fV1[1]
 704 //     << "xv = " << tAodTrack->Xv() << endl
 705 //     << "fv1[0] = " << fV1[0]  << endl
 706 //     << "yv = " << tAodTrack->Yv()  << endl
 707 //     << "fv1[1] = " << fV1[1]  << endl
 708 //     << "zv = " << tAodTrack->Zv()  << endl
 709 //     << "fv1[2] = " << fV1[2]  << endl
 710 //     << "impact[0] = " << impact[0]  << endl
 711 //     << "impact[1] = " << impact[1]  << endl
 712 //     << endl << endl ;
 713
 714   tFemtoTrack->SetCdd(covmat[0]);
 715   tFemtoTrack->SetCdz(covmat[1]);
 716   tFemtoTrack->SetCzz(covmat[2]);
 717   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
 718   tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
 719   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 720   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
 721   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
 722   tFemtoTrack->SetTPCsignalN(1);
 723   tFemtoTrack->SetTPCsignalS(1);
 724   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
 725
 726 //   if (tPWG2AODTrack) {
 727 //     // Copy the PWG2 specific information if it exists
 728 //     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
 729 //     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
 730
 731 //     double xtpc[3] = {0,0,0};
 732 //     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
 733 //     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 734 //     tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
 735 //     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 736 //   }
 737 //   else {
 738     // If not use dummy values
 739   tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
 740   tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
 741
 742   double xtpc[3] = {0,0,0};
 743   tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 744   tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 745   //   }
 746
 747   //   //  cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << "  " << tAodTrack->Zv() << "  " << tAodTrack->GetTPCNcls() << endl;
 748
 749
 750   int indexes[3];
 751   for (int ik=0; ik<3; ik++) {
 752     indexes[ik] = 0;
 753   }
 754   tFemtoTrack->SetKinkIndexes(indexes);
 755 }
 756
 757 void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
 758 {
 759   //tFemtoV0->SetdecayLengthV0(tAODv0->DecayLengthV0()); //wrocic do tego
 760   tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
 761   tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
 762   tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
 763   AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
 764   tFemtoV0->SetdecayVertexV0(decayvertex);
 765   tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
 766   tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
 767   tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
 768   tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
 769   tFemtoV0->SetmomPosX(tAODv0->MomPosX());
 770   tFemtoV0->SetmomPosY(tAODv0->MomPosY());
 771   tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
 772   AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
 773   tFemtoV0->SetmomPos(mompos);
 774   tFemtoV0->SetmomNegX(tAODv0->MomNegX());
 775   tFemtoV0->SetmomNegY(tAODv0->MomNegY());
 776   tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
 777   AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
 778   tFemtoV0->SetmomNeg(momneg);
 779
 780   //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
 781   //void SettpcHitsPos(const int& i);
 782   //void SettpcHitsNeg(const int& i);
 783
 784   //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
 785   //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
 786
 787   tFemtoV0->SetmomV0X(tAODv0->MomV0X());
 788   tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
 789   tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
 790   AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
 791   tFemtoV0->SetmomV0(momv0);
 792   tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
 793   tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
 794   tFemtoV0->SeteLambda(tAODv0->ELambda());
 795   tFemtoV0->SeteK0Short(tAODv0->EK0Short());
 796   tFemtoV0->SetePosProton(tAODv0->EPosProton());
 797   tFemtoV0->SeteNegProton(tAODv0->ENegProton());
 798   tFemtoV0->SetmassLambda(tAODv0->MassLambda());
 799   tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
 800   tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
 801   tFemtoV0->SetrapLambda(tAODv0->RapLambda());
 802   tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
 803
 804   //void SetcTauLambda( float x);
 805   //void SetcTauK0Short( float x);
 806
 807   tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
 808   tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
 809   tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
 810   tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
 811   tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
 812   tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
 813
 814   tFemtoV0->SetidNeg(tAODv0->GetNegID());
 815   //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
 816   //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
 817   tFemtoV0->SetidPos(tAODv0->GetPosID());
 818
 819   tFemtoV0->SetEtaV0(tAODv0->Eta());
 820   tFemtoV0->SetPhiV0(tAODv0->Phi());
 821   tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
 822   //tFemtoV0->SetYV0(tAODv0->Y());
 823
 824   //void SetdedxNeg(float x);
 825   //void SeterrdedxNeg(float x);//Gael 04Fev2002
 826   //void SetlendedxNeg(float x);//Gael 04Fev2002
 827   //void SetdedxPos(float x);
 828   //void SeterrdedxPos(float x);//Gael 04Fev2002
 829   //void SetlendedxPos(float x);//Gael 04Fev2002
 830
 831   tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
 832   tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
 833
 834   AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
 835   AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
 836
 837   if(trackpos && trackneg)
 838     {
 839       //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
 840       //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
 841       tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
 842       tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
 843       tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
 844       tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
 845       tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
 846       tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
 847       tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
 848       tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
 849       tFemtoV0->SetStatusPos(trackpos->GetStatus());
 850       tFemtoV0->SetStatusNeg(trackneg->GetStatus());
 851
 852       tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
 853       tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
 854       tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
 855       tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
 856       tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
 857       tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
 858
 859
 860       if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFmismatch)>0)
 861         {
 862           if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFmismatch)>0)
 863             {
 864               tFemtoV0->SetPosNSigmaTOFK(-9999);
 865               tFemtoV0->SetNegNSigmaTOFK(-9999);
 866               tFemtoV0->SetPosNSigmaTOFP(-9999);
 867               tFemtoV0->SetNegNSigmaTOFP(-9999);
 868               tFemtoV0->SetPosNSigmaTOFPi(-9999);
 869               tFemtoV0->SetNegNSigmaTOFPi(-9999);
 870             }
 871         }
 872       else
 873         {
 874           tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
 875           tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
 876           tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
 877           tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
 878           tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
 879           tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
 880         }
 881     }
 882   else
 883     {
 884       tFemtoV0->SetStatusPos(999);
 885       tFemtoV0->SetStatusNeg(999);
 886     }
 887   tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
 888 }
 889
 890 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
 891 {
 892   fFilterBit = (1 << (ibit));
 893 }
 894
 895 void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
 896 {
 897   fReadMC = a;
 898 }
 899
 900
 901 void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
 902 {
 903   fReadV0 = a;
 904 }
 905
 906 AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
 907 {
 908   if (aLabel < 0) return 0;
 909   AliAODMCParticle *aodP;
 910   Int_t posstack = 0;
 911   if (aLabel > mcP->GetEntries())
 912     posstack = mcP->GetEntries();
 913   else
 914     posstack = aLabel;
 915
 916   aodP = (AliAODMCParticle *) mcP->At(posstack);
 917   if (aodP->GetLabel() > posstack) {
 918     do {
 919       aodP = (AliAODMCParticle *) mcP->At(posstack);
 920       if (aodP->GetLabel() == aLabel) return aodP;
 921       posstack--;
 922     }
 923     while (posstack > 0);
 924   }
 925   else {
 926     do {
 927       aodP = (AliAODMCParticle *) mcP->At(posstack);
 928       if (aodP->GetLabel() == aLabel) return aodP;
 929       posstack++;
 930     }
 931     while (posstack < mcP->GetEntries());
 932   }
 933
 934   return 0;
 935 }
 936
 937 void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
 938                                                  AliFemtoTrack *tFemtoTrack)
 939 {
 940   double aodpid[10];
 941   tAodTrack->GetPID(aodpid);
 942   tFemtoTrack->SetPidProbElectron(aodpid[0]);
 943   tFemtoTrack->SetPidProbMuon(aodpid[1]);
 944   tFemtoTrack->SetPidProbPion(aodpid[2]);
 945   tFemtoTrack->SetPidProbKaon(aodpid[3]);
 946   tFemtoTrack->SetPidProbProton(aodpid[4]);
 947
 948   aodpid[0] = -100000.0;
 949   aodpid[1] = -100000.0;
 950   aodpid[2] = -100000.0;
 951   aodpid[3] = -100000.0;
 952   aodpid[4] = -100000.0;
 953
 954   double tTOF = 0.0;
 955
 956   if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) {  //AliESDtrack::kTOFpid=0x8000
 957     tTOF = tAodTrack->GetTOFsignal();
 958     tAodTrack->GetIntegratedTimes(aodpid);
 959   }
 960
 961   tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
 962
 963   //////  TPC ////////////////////////////////////////////
 964
 965   float nsigmaTPCK=-1000.;
 966   float nsigmaTPCPi=-1000.;
 967   float nsigmaTPCP=-1000.;
 968
 969   //   cout<<"in reader fESDpid"<<fESDpid<<endl;
 970
 971   if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
 972     nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
 973     nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
 974     nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
 975   }
 976
 977   tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
 978   tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
 979   tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
 980
 981   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 982   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
 983   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
 984
 985   tFemtoTrack->SetTPCsignalN(1);
 986   tFemtoTrack->SetTPCsignalS(1);
 987   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
 988
 989   ///////TOF//////////////////////
 990
 991     float vp=-1000.;
 992     float nsigmaTOFPi=-1000.;
 993     float nsigmaTOFK=-1000.;
 994     float nsigmaTOFP=-1000.;
 995
 996     if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
 997         (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
 998         (tAodTrack->GetStatus() & AliESDtrack::kTIME) && //AliESDtrack::kTIME=0x80000000
 999         !(tAodTrack->GetStatus() & AliESDtrack::kTOFmismatch)) //AliESDtrack::kTOFmismatch=0x100000
1000       {
1001         if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
1002           {
1003
1004             nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1005             nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1006             nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1007
1008             Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1009             Double_t tof=tAodTrack->GetTOFsignal();
1010             if(tof > 0.) vp=len/tof/0.03;
1011           }
1012       }
1013     tFemtoTrack->SetVTOF(vp);
1014     tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1015     tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1016     tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1017
1018
1019     //////////////////////////////////////
1020
1021 }
1022
1023 void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1024 {
1025   fCentRange[0] = min; fCentRange[1] = max;
1026   fUsePreCent = 1;
1027   fNoCentrality = 0;
1028 }
1029
1030 void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1031 {
1032   fNoCentrality = anocent;
1033   if (fNoCentrality) fUsePreCent = 0;
1034 }
1035
1036 void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1037 {
1038   fAODpidUtil = aAODpidUtil;
1039   //  printf("fAODpidUtil: %x\n",fAODpidUtil);
1040 }
1041
1042
1043
1044
1045