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