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