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