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