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