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 //constructor with 0 parameters , look at default settings
  41 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
  42   fNumberofEvent(0),
  43   fCurEvent(0),
  44   fEvent(0x0),
  45   fAllTrue(160),
  46   fAllFalse(160),
  47   fFilterBit(0),
  48   //  fPWG2AODTracks(0x0),
  49   fReadMC(0),
  50   fUsePreCent(0),
  51   fAODpidUtil(0),
  52   fInputFile(" "),
  53   fFileName(" "),
  54   fTree(0x0),
  55   fAodFile(0x0)
  56 {
  57   // default constructor
  58   fAllTrue.ResetAllBits(kTRUE);
  59   fAllFalse.ResetAllBits(kFALSE);
  60   fCentRange[0] = 0;
  61   fCentRange[1] = 1000;
  62 }
  63
  64 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
  65   AliFemtoEventReader(),
  66   fNumberofEvent(0),
  67   fCurEvent(0),
  68   fEvent(0x0),
  69   fAllTrue(160),
  70   fAllFalse(160),
  71   fFilterBit(0),
  72   //  fPWG2AODTracks(0x0),
  73   fReadMC(0),
  74   fUsePreCent(0),
  75   fAODpidUtil(0),
  76   fInputFile(" "),
  77   fFileName(" "),
  78   fTree(0x0),
  79   fAodFile(0x0)
  80 {
  81   // copy constructor
  82   fInputFile = aReader.fInputFile;
  83   fFileName  = aReader.fFileName;
  84   fNumberofEvent = aReader.fNumberofEvent;
  85   fCurEvent = aReader.fCurEvent;
  86   fEvent = new AliAODEvent();
  87   fAodFile = new TFile(aReader.fAodFile->GetName());
  88   fAllTrue.ResetAllBits(kTRUE);
  89   fAllFalse.ResetAllBits(kFALSE);
  90   fFilterBit = aReader.fFilterBit;
  91   //  fPWG2AODTracks = aReader.fPWG2AODTracks;
  92   fAODpidUtil = aReader.fAODpidUtil;
  93   fCentRange[0] = aReader.fCentRange[0];
  94   fCentRange[1] = aReader.fCentRange[1];
  95 }
  96 //__________________
  97 //Destructor
  98 AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
  99 {
 100   // destructor
 101   delete fTree;
 102   delete fEvent;
 103   delete fAodFile;
 104 //   if (fPWG2AODTracks) {
 105 //     fPWG2AODTracks->Delete();
 106 //     delete fPWG2AODTracks;
 107 //   }
 108 }
 109
 110 //__________________
 111 AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
 112 {
 113   // assignment operator
 114   if (this == &aReader)
 115     return *this;
 116
 117   fInputFile = aReader.fInputFile;
 118   fFileName  = aReader.fFileName;
 119   fNumberofEvent = aReader.fNumberofEvent;
 120   fCurEvent = aReader.fCurEvent;
 121   if (fTree) delete fTree;
 122   if (fEvent) delete fEvent;
 123   fEvent = new AliAODEvent();
 124   if (fAodFile) delete fAodFile;
 125   fAodFile = new TFile(aReader.fAodFile->GetName());
 126   fAllTrue.ResetAllBits(kTRUE);
 127   fAllFalse.ResetAllBits(kFALSE);
 128   fFilterBit = aReader.fFilterBit;
 129   //  fPWG2AODTracks = aReader.fPWG2AODTracks;
 130   fAODpidUtil = aReader.fAODpidUtil;
 131   fCentRange[0] = aReader.fCentRange[0];
 132   fCentRange[1] = aReader.fCentRange[1];
 133
 134   return *this;
 135 }
 136 //__________________
 137 AliFemtoString AliFemtoEventReaderAOD::Report()
 138 {
 139   // create reader report
 140   AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
 141   return temp;
 142 }
 143
 144 //__________________
 145 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
 146 {
 147   //setting the name of file where names of AOD file are written
 148   //it takes only this files which have good trees
 149   char buffer[256];
 150   fInputFile=string(inputFile);
 151   ifstream infile(inputFile);
 152
 153   fTree = new TChain("aodTree");
 154
 155   if(infile.good()==true)
 156     {
 157       //checking if all give files have good tree inside
 158       while (infile.eof()==false)
 159         {
 160           infile.getline(buffer,256);
 161           TFile *aodFile=TFile::Open(buffer,"READ");
 162           if (aodFile!=0x0)
 163             {
 164               TTree* tree = (TTree*) aodFile->Get("aodTree");
 165               if (tree!=0x0)
 166                 {
 167                   //              cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
 168                   fTree->AddFile(buffer);
 169                   delete tree;
 170                 }
 171               aodFile->Close();
 172             }
 173           delete aodFile;
 174         }
 175     }
 176 }
 177
 178 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
 179 {
 180   // read in a next hbt event from the chain
 181   // convert it to AliFemtoEvent and return
 182   // for further analysis
 183   AliFemtoEvent *hbtEvent = 0;
 184   cout<<"reader"<<endl;
 185   if (fCurEvent==fNumberofEvent)//open next file
 186     {
 187       if(fNumberofEvent==0)
 188         {
 189           fEvent=new AliAODEvent();
 190           fEvent->ReadFromTree(fTree);
 191
 192           // Check for the existence of the additional information
 193 //        fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
 194
 195 //        if (fPWG2AODTracks) {
 196 //          cout << "Found additional PWG2 specific information in the AOD!" << endl;
 197 //          cout << "Reading only tracks with the additional information" << endl;
 198 //        }
 199
 200           fNumberofEvent=fTree->GetEntries();
 201           //      cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
 202           fCurEvent=0;
 203         }
 204       else //no more data to read
 205         {
 206           cout<<"no more files "<<hbtEvent<<endl;
 207           fReaderStatus=1;
 208           return hbtEvent;
 209         }
 210     }
 211
 212   cout<<"starting to read event "<<fCurEvent<<endl;
 213   fTree->GetEvent(fCurEvent);//getting next event
 214   //  cout << "Read event " << fEvent << " from file " << fTree << endl;
 215
 216   hbtEvent = new AliFemtoEvent;
 217
 218   CopyAODtoFemtoEvent(hbtEvent);
 219   fCurEvent++;
 220
 221
 222   return hbtEvent;
 223 }
 224
 225 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 226 {
 227
 228   // A function that reads in the AOD event
 229   // and transfers the neccessary information into
 230   // the internal AliFemtoEvent
 231
 232   // setting global event characteristics
 233   tEvent->SetRunNumber(fEvent->GetRunNumber());
 234   tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 235   tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 236   tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 237   tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 238   tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 239   tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
 240   tEvent->SetZDCParticipants(0);
 241   tEvent->SetTriggerMask(fEvent->GetTriggerMask());
 242   tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 243
 244   // Attempt to access MC header
 245   AliAODMCHeader *mcH;
 246   TClonesArray *mcP=0;
 247   if (fReadMC) {
 248     mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
 249     if (!mcH) {
 250       cout << "AOD MC information requested, but no header found!" << endl;
 251     }
 252
 253     mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
 254     if (!mcP) {
 255       cout << "AOD MC information requested, but no particle array found!" << endl;
 256     }
 257   }
 258
 259   tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
 260
 261   Int_t *motherids=0;
 262   if (mcP) {
 263     motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
 264     for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
 265
 266     // Read in mother ids
 267     AliAODMCParticle *motherpart;
 268     for (int ip=0; ip<mcP->GetEntries(); ip++) {
 269       motherpart = (AliAODMCParticle *) mcP->At(ip);
 270       if (motherpart->GetDaughter(0) > 0)
 271         motherids[motherpart->GetDaughter(0)] = ip;
 272       if (motherpart->GetDaughter(1) > 0)
 273         motherids[motherpart->GetDaughter(1)] = ip;
 274     }
 275   }
 276
 277   // Primary Vertex position
 278   double fV1[3];
 279   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 280
 281   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 282   tEvent->SetPrimVertPos(vertex);
 283
 284   //starting to reading tracks
 285   int nofTracks=0;  //number of reconstructed tracks in event
 286
 287   // Check to see whether the additional info exists
 288 //   if (fPWG2AODTracks)
 289 //     nofTracks=fPWG2AODTracks->GetEntries();
 290 //   else
 291     nofTracks=fEvent->GetNumberOfTracks();
 292   cout<<"nofTracks: "<<nofTracks<<endl;
 293
 294   AliCentrality *cent = fEvent->GetCentrality();
 295
 296   if (cent && fUsePreCent) {
 297     if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
 298         (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
 299       {
 300         cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
 301
 302         return;
 303       }
 304   }
 305
 306   int realnofTracks=0;   // number of track which we use in a analysis
 307   int tracksPrim=0;
 308
 309   int labels[20000];
 310   for (int il=0; il<20000; il++) labels[il] = -1;
 311
 312   // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
 313   for (int i=0;i<nofTracks;i++) {
 314     const AliAODTrack *aodtrack=fEvent->GetTrack(i);
 315     if (!aodtrack->TestFilterBit(fFilterBit)) {
 316       labels[aodtrack->GetID()] = i;
 317     }
 318   }
 319
 320   //  int tNormMult = 0;
 321   for (int i=0;i<nofTracks;i++)
 322     {
 323       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 324
 325 //       if (fPWG2AODTracks) {
 326 //      // Read tracks from the additional pwg2 specific AOD part
 327 //      // if they exist
 328 //      // Note that in that case all the AOD tracks without the
 329 //      // additional information will be ignored !
 330 //      AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
 331
 332 //      // Getting the AOD track through the ref of the additional info
 333 //      AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
 334 //      if (!aodtrack->TestFilterBit(fFilterBit)) {
 335 //        delete trackCopy;
 336 //        continue;
 337 //      }
 338
 339
 340 //      if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
 341 //        if (aodtrack->Chi2perNDF() < 6.0)
 342 //          if (aodtrack->Eta() < 0.9)
 343 //            tNormMult++;
 344
 345
 346 //      CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
 347
 348 //      if (mcP) {
 349 //        // Fill the hidden information with the simulated data
 350 //        //      Int_t pLabel = aodtrack->GetLabel();
 351 //        AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 352
 353 //        // Check the mother information
 354
 355 //        // Using the new way of storing the freeze-out information
 356 //        // Final state particle is stored twice on the stack
 357 //        // one copy (mother) is stored with original freeze-out information
 358 //        //   and is not tracked
 359 //        // the other one (daughter) is stored with primary vertex position
 360 //        //   and is tracked
 361
 362 //        // Freeze-out coordinates
 363 //        double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 364 //        fpx = tPart->Xv() - fV1[0];
 365 //        fpy = tPart->Yv() - fV1[1];
 366 //        fpz = tPart->Zv() - fV1[2];
 367 //        fpt = tPart->T();
 368
 369 //        AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 370 //        tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 371
 372 //        fpx *= 1e13;
 373 //        fpy *= 1e13;
 374 //        fpz *= 1e13;
 375 //        fpt *= 1e13;
 376
 377 //        //      cout << "Looking for mother ids " << endl;
 378 //        if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 379 //          //  cout << "Got mother id" << endl;
 380 //          AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 381 //          // Check if this is the same particle stored twice on the stack
 382 //          if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 383 //            // It is the same particle
 384 //            // Read in the original freeze-out information
 385 //            // and convert it from to [fm]
 386
 387 //            // EPOS style
 388 //            //          fpx = mother->Xv()*1e13*0.197327;
 389 //            //          fpy = mother->Yv()*1e13*0.197327;
 390 //            //          fpz = mother->Zv()*1e13*0.197327;
 391 //            //          fpt = mother->T() *1e13*0.197327*0.5;
 392
 393
 394 //            // Therminator style
 395 //            fpx = mother->Xv()*1e13;
 396 //            fpy = mother->Yv()*1e13;
 397 //            fpz = mother->Zv()*1e13;
 398 //            fpt = mother->T() *1e13*3e10;
 399
 400 //          }
 401 //        }
 402
 403 //        //       if (fRotateToEventPlane) {
 404 //        //    double tPhi = TMath::ATan2(fpy, fpx);
 405 //        //    double tRad = TMath::Hypot(fpx, fpy);
 406
 407 //        //    fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 408 //        //    fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 409 //        //       }
 410
 411 //        tInfo->SetPDGPid(tPart->GetPdgCode());
 412
 413 //        //      if (fRotateToEventPlane) {
 414 //        //        double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 415 //        //        double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 416
 417 //        //        tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 418 //        //                               tRad*TMath::Sin(tPhi - tReactionPlane),
 419 //        //                               tPart->Pz());
 420 //        //      }
 421 //        //       else
 422 //        tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 423 //        Double_t mass2 = (tPart->E() *tPart->E() -
 424 //                          tPart->Px()*tPart->Px() -
 425 //                          tPart->Py()*tPart->Py() -
 426 //                          tPart->Pz()*tPart->Pz());
 427 //        if (mass2>0.0)
 428 //          tInfo->SetMass(TMath::Sqrt(mass2));
 429 //        else
 430 //          tInfo->SetMass(0.0);
 431
 432 //        tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 433 //        trackCopy->SetHiddenInfo(tInfo);
 434
 435 //      }
 436
 437 //      double pxyz[3];
 438 //      aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 439 //      const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 440 //      // Check the sanity of the tracks - reject zero momentum tracks
 441 //      if (ktP.Mag() == 0) {
 442 //        delete trackCopy;
 443 //        continue;
 444 //      }
 445 //       }
 446 //       else {
 447         // No additional information exists
 448         // Read in the normal AliAODTracks
 449
 450         //      const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 451         AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 452
 453         if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
 454
 455         if (!aodtrack->TestFilterBit(fFilterBit)) {
 456           delete trackCopy;
 457           continue;
 458         }
 459
 460         CopyAODtoFemtoTrack(aodtrack, trackCopy);
 461
 462         // copying PID information from the correspondent track
 463         //      const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
 464         AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
 465         CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
 466
 467         if (mcP) {
 468           // Fill the hidden information with the simulated data
 469           //      Int_t pLabel = aodtrack->GetLabel();
 470           AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 471
 472           AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 473           double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 474           if (!tPart) {
 475             fpx = fV1[0];
 476             fpy = fV1[1];
 477             fpz = fV1[2];
 478             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 479             tInfo->SetPDGPid(0);
 480             tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
 481             tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
 482             tInfo->SetMass(0);
 483           }
 484           else {
 485             // Check the mother information
 486
 487             // Using the new way of storing the freeze-out information
 488             // Final state particle is stored twice on the stack
 489             // one copy (mother) is stored with original freeze-out information
 490             //   and is not tracked
 491             // the other one (daughter) is stored with primary vertex position
 492             //   and is tracked
 493
 494             // Freeze-out coordinates
 495             fpx = tPart->Xv() - fV1[0];
 496             fpy = tPart->Yv() - fV1[1];
 497             fpz = tPart->Zv() - fV1[2];
 498             //    fpt = tPart->T();
 499
 500             tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 501
 502             fpx *= 1e13;
 503             fpy *= 1e13;
 504             fpz *= 1e13;
 505             //    fpt *= 1e13;
 506
 507             //      cout << "Looking for mother ids " << endl;
 508             if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 509               //        cout << "Got mother id" << endl;
 510               AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
 511               // Check if this is the same particle stored twice on the stack
 512               if (mother) {
 513                 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 514                   // It is the same particle
 515                   // Read in the original freeze-out information
 516                   // and convert it from to [fm]
 517
 518                   // EPOS style
 519                   //      fpx = mother->Xv()*1e13*0.197327;
 520                   //      fpy = mother->Yv()*1e13*0.197327;
 521                   //      fpz = mother->Zv()*1e13*0.197327;
 522                   //      fpt = mother->T() *1e13*0.197327*0.5;
 523
 524
 525                   // Therminator style
 526                   fpx = mother->Xv()*1e13;
 527                   fpy = mother->Yv()*1e13;
 528                   fpz = mother->Zv()*1e13;
 529                   //          fpt = mother->T() *1e13*3e10;
 530
 531                 }
 532               }
 533             }
 534
 535             //       if (fRotateToEventPlane) {
 536             //  double tPhi = TMath::ATan2(fpy, fpx);
 537             //  double tRad = TMath::Hypot(fpx, fpy);
 538
 539             //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 540             //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 541             //       }
 542
 543             tInfo->SetPDGPid(tPart->GetPdgCode());
 544
 545             //    if (fRotateToEventPlane) {
 546             //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 547             //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 548
 549             //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 550             //                             tRad*TMath::Sin(tPhi - tReactionPlane),
 551             //                             tPart->Pz());
 552             //    }
 553             //       else
 554             tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 555             Double_t mass2 = (tPart->E() *tPart->E() -
 556                               tPart->Px()*tPart->Px() -
 557                               tPart->Py()*tPart->Py() -
 558                               tPart->Pz()*tPart->Pz());
 559             if (mass2>0.0)
 560               tInfo->SetMass(TMath::Sqrt(mass2));
 561             else
 562               tInfo->SetMass(0.0);
 563
 564             tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 565           }
 566           trackCopy->SetHiddenInfo(tInfo);
 567         }
 568
 569         double pxyz[3];
 570         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 571         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 572         // Check the sanity of the tracks - reject zero momentum tracks
 573         if (ktP.Mag() == 0) {
 574           delete trackCopy;
 575           continue;
 576         }
 577         //    }
 578
 579
 580         tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 581         realnofTracks++;//real number of tracks
 582     }
 583
 584   tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 585   tEvent->SetNormalizedMult(tracksPrim);
 586
 587   //  AliCentrality *cent = fEvent->GetCentrality();
 588   if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
 589   //  if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
 590
 591   if (cent) {
 592     tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 593     //    tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 594     tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 595     //    tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 596   }
 597
 598
 599   if (mcP) delete [] motherids;
 600
 601   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 602 }
 603
 604 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
 605                                                  AliFemtoTrack *tFemtoTrack
 606                                                  //                                              AliPWG2AODTrack *tPWG2AODTrack
 607                                                  )
 608 {
 609   // Copy the track information from the AOD into the internal AliFemtoTrack
 610   // If it exists, use the additional information from the PWG2 AOD
 611
 612   // Primary Vertex position
 613   double fV1[3];
 614   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 615   //  fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 616
 617   tFemtoTrack->SetCharge(tAodTrack->Charge());
 618
 619   double pxyz[3];
 620   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
 621   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 622   tFemtoTrack->SetP(v);//setting momentum
 623   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 624   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 625   //setting track helix
 626   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 627   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
 628   tFemtoTrack->SetHelix(helix);
 629
 630   // Flags
 631   tFemtoTrack->SetTrackId(tAodTrack->GetID());
 632   tFemtoTrack->SetFlags(tAodTrack->GetFlags());
 633   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
 634
 635   // Track quality information
 636   float covmat[6];
 637   tAodTrack->GetCovMatrix(covmat);
 638
 639   double impact[2];
 640   double covimpact[3];
 641
 642   if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
 643     //cout << "sth went wrong with dca propagation" << endl;
 644     tFemtoTrack->SetImpactD(-1000.0);
 645     tFemtoTrack->SetImpactZ(-1000.0);
 646
 647   }
 648   else {
 649     tFemtoTrack->SetImpactD(impact[0]);
 650     tFemtoTrack->SetImpactZ(impact[1]);
 651   }
 652
 653   //   if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
 654   //     tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
 655   //   else
 656   //     tFemtoTrack->SetImpactD(0.0);
 657   //   tFemtoTrack->SetImpactD(tAodTrack->DCA());
 658
 659   //   tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
 660
 661
 662   //   tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
 663   //   tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
 664
 665
 666   //   cout
 667     //    << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
 668     //    << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
 669     //    << tAodTrack->Yv() - fV1[1]
 670 //     << "xv = " << tAodTrack->Xv() << endl
 671 //     << "fv1[0] = " << fV1[0]  << endl
 672 //     << "yv = " << tAodTrack->Yv()  << endl
 673 //     << "fv1[1] = " << fV1[1]  << endl
 674 //     << "zv = " << tAodTrack->Zv()  << endl
 675 //     << "fv1[2] = " << fV1[2]  << endl
 676 //     << "impact[0] = " << impact[0]  << endl
 677 //     << "impact[1] = " << impact[1]  << endl
 678 //     << endl << endl ;
 679
 680   tFemtoTrack->SetCdd(covmat[0]);
 681   tFemtoTrack->SetCdz(covmat[1]);
 682   tFemtoTrack->SetCzz(covmat[2]);
 683   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
 684   tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
 685   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 686   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
 687   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
 688   tFemtoTrack->SetTPCsignalN(1);
 689   tFemtoTrack->SetTPCsignalS(1);
 690   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
 691
 692 //   if (tPWG2AODTrack) {
 693 //     // Copy the PWG2 specific information if it exists
 694 //     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
 695 //     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
 696
 697 //     double xtpc[3] = {0,0,0};
 698 //     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
 699 //     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 700 //     tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
 701 //     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 702 //   }
 703 //   else {
 704     // If not use dummy values
 705   tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
 706   tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
 707
 708   double xtpc[3] = {0,0,0};
 709   tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 710   tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 711   //   }
 712
 713   //   //  cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << "  " << tAodTrack->Zv() << "  " << tAodTrack->GetTPCNcls() << endl;
 714
 715
 716   int indexes[3];
 717   for (int ik=0; ik<3; ik++) {
 718     indexes[ik] = 0;
 719   }
 720   tFemtoTrack->SetKinkIndexes(indexes);
 721 }
 722
 723 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
 724 {
 725   fFilterBit = (1 << (ibit));
 726 }
 727
 728 void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
 729 {
 730   fReadMC = a;
 731 }
 732
 733 AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
 734 {
 735   if (aLabel < 0) return 0;
 736   AliAODMCParticle *aodP;
 737   Int_t posstack = 0;
 738   if (aLabel > mcP->GetEntries())
 739     posstack = mcP->GetEntries();
 740   else
 741     posstack = aLabel;
 742
 743   aodP = (AliAODMCParticle *) mcP->At(posstack);
 744   if (aodP->GetLabel() > posstack) {
 745     do {
 746       aodP = (AliAODMCParticle *) mcP->At(posstack);
 747       if (aodP->GetLabel() == aLabel) return aodP;
 748       posstack--;
 749     }
 750     while (posstack > 0);
 751   }
 752   else {
 753     do {
 754       aodP = (AliAODMCParticle *) mcP->At(posstack);
 755       if (aodP->GetLabel() == aLabel) return aodP;
 756       posstack++;
 757     }
 758     while (posstack < mcP->GetEntries());
 759   }
 760
 761   return 0;
 762 }
 763
 764 void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
 765                                                  AliFemtoTrack *tFemtoTrack)
 766 {
 767   double aodpid[10];
 768   tAodTrack->GetPID(aodpid);
 769   tFemtoTrack->SetPidProbElectron(aodpid[0]);
 770   tFemtoTrack->SetPidProbMuon(aodpid[1]);
 771   tFemtoTrack->SetPidProbPion(aodpid[2]);
 772   tFemtoTrack->SetPidProbKaon(aodpid[3]);
 773   tFemtoTrack->SetPidProbProton(aodpid[4]);
 774
 775   aodpid[0] = -100000.0;
 776   aodpid[1] = -100000.0;
 777   aodpid[2] = -100000.0;
 778   aodpid[3] = -100000.0;
 779   aodpid[4] = -100000.0;
 780
 781   double tTOF = 0.0;
 782
 783   if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) {  //AliESDtrack::kTOFpid=0x8000
 784     tTOF = tAodTrack->GetTOFsignal();
 785     tAodTrack->GetIntegratedTimes(aodpid);
 786   }
 787
 788   tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
 789
 790   //////  TPC ////////////////////////////////////////////
 791
 792   float nsigmaTPCK=-1000.;
 793   float nsigmaTPCPi=-1000.;
 794   float nsigmaTPCP=-1000.;
 795
 796   //   cout<<"in reader fESDpid"<<fESDpid<<endl;
 797
 798   if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
 799     nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
 800     nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
 801     nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
 802   }
 803
 804   tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
 805   tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
 806   tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
 807
 808   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 809   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
 810   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
 811
 812   tFemtoTrack->SetTPCsignalN(1);
 813   tFemtoTrack->SetTPCsignalS(1);
 814   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
 815
 816   ///////TOF//////////////////////
 817
 818     float vp=-1000.;
 819     float nsigmaTOFPi=-1000.;
 820     float nsigmaTOFK=-1000.;
 821     float nsigmaTOFP=-1000.;
 822
 823     if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
 824         (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
 825         (tAodTrack->GetStatus() & AliESDtrack::kTIME) && //AliESDtrack::kTIME=0x80000000
 826         !(tAodTrack->GetStatus() & AliESDtrack::kTOFmismatch)) //AliESDtrack::kTOFmismatch=0x100000
 827       {
 828         if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
 829           {
 830
 831             nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
 832             nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
 833             nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
 834
 835             Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
 836             Double_t tof=tAodTrack->GetTOFsignal();
 837             if(tof > 0.) vp=len/tof/0.03;
 838           }
 839       }
 840     tFemtoTrack->SetVTOF(vp);
 841     tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
 842     tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
 843     tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
 844
 845
 846     //////////////////////////////////////
 847
 848 }
 849
 850 void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
 851 {
 852   fCentRange[0] = min; fCentRange[1] = max;
 853   fUsePreCent = 1;
 854 }
 855
 856
 857 void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
 858 {
 859   fAODpidUtil = aAODpidUtil;
 860   //  printf("fAODpidUtil: %x\n",fAODpidUtil);
 861 }
 862
 863
 864
 865
 866