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