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