PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx

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