PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and   //
   4 // the model Kinematics information tailored for the Task framework and the   //
   5 // Reads in AliESDfriend to create shared hit/quality information             //
   6 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
   7 //                                                                            //
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoEventReaderESDChainKine.h"
  10
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "TList.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDtrack.h"
  16 #include "AliESDVertex.h"
  17
  18 #include "AliMultiplicity.h"
  19 #include "AliCentrality.h"
  20 #include "AliEventplane.h"
  21
  22 #include "AliFmPhysicalHelixD.h"
  23 #include "AliFmThreeVectorF.h"
  24
  25 #include "SystemOfUnits.h"
  26
  27 #include "AliFemtoEvent.h"
  28
  29 #include "TParticle.h"
  30 #include "AliFemtoModelHiddenInfo.h"
  31 #include "AliFemtoModelGlobalHiddenInfo.h"
  32 #include "AliGenHijingEventHeader.h"
  33 #include "AliGenCocktailEventHeader.h"
  34
  35 #include "AliVertexerTracks.h"
  36
  37 #include "AliPID.h"
  38
  39 ClassImp(AliFemtoEventReaderESDChainKine)
  40
  41 #if !(ST_NO_NAMESPACES)
  42 using namespace units;
  43 #endif
  44
  45 using namespace std;
  46 //____________________________
  47 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
  48   fFileName(" "),
  49   fConstrained(true),
  50   fReadInner(false),
  51   fUseTPCOnly(false),
  52   fNumberofEvent(0),
  53   fCurEvent(0),
  54   fCurFile(0),
  55   fEvent(0x0),
  56   fStack(0x0),
  57   fGenHeader(0x0),
  58   fTrackType(kGlobal),
  59   fEstEventMult(kReferenceITSTPC),
  60   fRotateToEventPlane(0),
  61   fESDpid(0),
  62   fIsPidOwner(0),
  63   fMagFieldSign(0),
  64   fReadV0(0),
  65   isKaonAnalysis(kFALSE),
  66   isProtonAnalysis(kFALSE),
  67   isPionAnalysis(kFALSE),
  68   fOnlyPrimaries(kFALSE)
  69
  70 {
  71   //constructor with 0 parameters , look at default settings
  72 }
  73
  74 //__________________
  75 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  76   AliFemtoEventReader(aReader),
  77   fFileName(" "),
  78   fConstrained(true),
  79   fReadInner(false),
  80   fUseTPCOnly(false),
  81   fNumberofEvent(0),
  82   fCurEvent(0),
  83   fCurFile(0),
  84   fEvent(0x0),
  85   fStack(0x0),
  86   fGenHeader(0x0),
  87   fTrackType(kGlobal),
  88   fEstEventMult(kReferenceITSTPC),
  89   fRotateToEventPlane(0),
  90   fESDpid(0),
  91   fIsPidOwner(0),
  92   fMagFieldSign(0),
  93   fReadV0(0),
  94   isKaonAnalysis(kFALSE),
  95   isProtonAnalysis(kFALSE),
  96   isPionAnalysis(kFALSE),
  97   fOnlyPrimaries(kFALSE)
  98
  99
 100 {
 101   // Copy constructor
 102   fConstrained = aReader.fConstrained;
 103   fReadInner = aReader.fReadInner;
 104   fUseTPCOnly = aReader.fUseTPCOnly;
 105   fNumberofEvent = aReader.fNumberofEvent;
 106   fCurEvent = aReader.fCurEvent;
 107   fCurFile = aReader.fCurFile;
 108   fEvent = new AliESDEvent();
 109   fStack = aReader.fStack;
 110   fTrackType = aReader.fTrackType;
 111   fEstEventMult = aReader.fEstEventMult;
 112   fRotateToEventPlane = aReader.fRotateToEventPlane;
 113   fESDpid = aReader.fESDpid;
 114   fIsPidOwner = aReader.fIsPidOwner;
 115   fMagFieldSign = aReader.fMagFieldSign;
 116   fReadV0 = aReader.fReadV0;
 117   isKaonAnalysis = aReader.isKaonAnalysis;
 118   isProtonAnalysis = aReader.isProtonAnalysis;
 119   isPionAnalysis = aReader.isPionAnalysis;
 120   fOnlyPrimaries = aReader.fOnlyPrimaries;
 121
 122 }
 123 //__________________
 124 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
 125 {
 126   //Destructor
 127   delete fEvent;
 128 }
 129
 130 //__________________
 131 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
 132 {
 133   // Assignment operator
 134   if (this == &aReader)
 135     return *this;
 136
 137   fConstrained = aReader.fConstrained;
 138   fUseTPCOnly = aReader.fUseTPCOnly;
 139   fNumberofEvent = aReader.fNumberofEvent;
 140   fCurEvent = aReader.fCurEvent;
 141   fCurFile = aReader.fCurFile;
 142   if (fEvent) delete fEvent;
 143   fEvent = new AliESDEvent();
 144   fStack = aReader.fStack;
 145   fGenHeader = aReader.fGenHeader;
 146   fRotateToEventPlane = aReader.fRotateToEventPlane;
 147   fESDpid = aReader.fESDpid;
 148   fIsPidOwner = aReader.fIsPidOwner;
 149   fMagFieldSign = aReader.fMagFieldSign;
 150   fReadV0 = aReader.fReadV0;
 151   isKaonAnalysis = aReader.isKaonAnalysis;
 152   isProtonAnalysis = aReader.isProtonAnalysis;
 153   isPionAnalysis = aReader.isPionAnalysis;
 154   fOnlyPrimaries = aReader.fOnlyPrimaries;
 155
 156   return *this;
 157 }
 158 //__________________
 159 // Simple report
 160 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
 161 {
 162   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
 163   return temp;
 164 }
 165
 166 //__________________
 167 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
 168 {
 169   // Select whether to read constrained or not constrained momentum
 170   fConstrained=constrained;
 171 }
 172 //__________________
 173 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
 174 {
 175   // Check whether we read constrained or not constrained momentum
 176   return fConstrained;
 177 }
 178 //__________________
 179 void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
 180 {
 181   fReadInner=readinner;
 182 }
 183
 184 bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
 185 {
 186   return fReadInner;
 187 }
 188
 189 //__________________
 190 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
 191 {
 192   fUseTPCOnly=usetpconly;
 193 }
 194
 195 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
 196 {
 197   return fUseTPCOnly;
 198 }
 199 void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
 200 {
 201   fEstEventMult = aType;
 202 }
 203 //__________________
 204 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 205 {
 206   // Get the event, read all the relevant information
 207   // and fill the AliFemtoEvent class
 208   // Returns a valid AliFemtoEvent
 209   AliFemtoEvent *hbtEvent = 0;
 210   string tFriendFileName;
 211
 212   // Get the friend information
 213   cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
 214   //  fEvent->SetESDfriend(fEventFriend);
 215
 216   hbtEvent = new AliFemtoEvent;
 217   //setting basic things
 218   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 219   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 220   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 221   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 222   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 223   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 224   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 225   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 226   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 227   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 228   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 229   //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 230
 231   if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
 232       (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
 233       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
 234       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
 235     hbtEvent->SetTriggerCluster(1);
 236   else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
 237            (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
 238     hbtEvent->SetTriggerCluster(2);
 239   else
 240     hbtEvent->SetTriggerCluster(0);
 241
 242   //Vertex
 243   double fV1[3];
 244   double fVCov[6];
 245   if (fUseTPCOnly) {
 246     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 247     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 248     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 249       fVCov[4] = -1001.0;
 250   }
 251   else {
 252     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 253     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 254     if (!fEvent->GetPrimaryVertex()->GetStatus())
 255       fVCov[4] = -1001.0;
 256   }
 257
 258   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 259   hbtEvent->SetPrimVertPos(vertex);
 260   hbtEvent->SetPrimVertCov(fVCov);
 261
 262   Double_t tReactionPlane = 0;
 263
 264   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 265   if (!hdh) {
 266     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 267     if (cdh) {
 268       TList *tGenHeaders = cdh->GetHeaders();
 269       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 270         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 271         if (hdh) break;
 272       }
 273     }
 274   }
 275
 276   if (hdh)
 277   {
 278     tReactionPlane = hdh->ReactionPlaneAngle();
 279     cout << "Got reaction plane " << tReactionPlane << endl;
 280   }
 281
 282   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 283
 284   Int_t spdetaonecount = 0;
 285
 286   for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
 287     if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
 288       spdetaonecount++;
 289
 290   //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
 291   hbtEvent->SetSPDMult(spdetaonecount);
 292
 293   //starting to reading tracks
 294   int nofTracks=0;  //number of reconstructed tracks in event
 295   nofTracks=fEvent->GetNumberOfTracks();
 296   int realnofTracks=0;//number of track which we use ina analysis
 297
 298   int tNormMult = 0;
 299   int tNormMultPos = 0;
 300   int tNormMultNeg = 0;
 301
 302   Float_t tTotalPt = 0.0;
 303
 304   Float_t b[2];
 305   Float_t bCov[3];
 306
 307   Int_t tTracklet=0, tITSTPC=0;
 308
 309   //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
 310
 311   hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
 312   hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
 313   //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
 314   hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 315
 316   Int_t *motherids;
 317   if (fStack) {
 318     motherids = new Int_t[fStack->GetNtrack()];
 319     for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 320
 321     // Read in mother ids
 322     TParticle *motherpart;
 323     for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 324       motherpart = fStack->Particle(ip);
 325       if (motherpart->GetDaughter(0) > 0)
 326         motherids[motherpart->GetDaughter(0)] = ip;
 327       if (motherpart->GetDaughter(1) > 0)
 328         motherids[motherpart->GetDaughter(1)] = ip;
 329
 330 //     if (motherpart->GetPdgCode() == 211) {
 331 //       cout << "Mother " << ip << " has daughters "
 332 //         << motherpart->GetDaughter(0) << " "
 333 //         << motherpart->GetDaughter(1) << " "
 334 //         << motherpart->Vx() << " "
 335 //         << motherpart->Vy() << " "
 336 //         << motherpart->Vz() << " "
 337 //         << endl;
 338
 339 //     }
 340     }
 341   }
 342   else {
 343     printf ("No Stack ???\n");
 344     delete hbtEvent;
 345     return 0;
 346   }
 347
 348   for (int i=0;i<nofTracks;i++)
 349   {
 350     //      cout << "Reading track " << i << endl;
 351     bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 352
 353
 354     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 355     //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 356
 357
 358     if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
 359         (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
 360       if (esdtrack->GetTPCNcls() > 70)
 361         if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 362           if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 363             esdtrack->GetImpactParameters(b,bCov);
 364             if ((b[0]<0.2) && (b[1] < 0.25)) {
 365               tNormMult++;
 366               tTotalPt += esdtrack->Pt();
 367             }
 368           }
 369         }
 370     }
 371     else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
 372       if (esdtrack->GetTPCNcls() > 100)
 373         if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 374           if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 375             esdtrack->GetImpactParameters(b,bCov);
 376             if ((b[0]<2.4) && (b[1] < 3.2)) {
 377               tNormMult++;
 378               tTotalPt += esdtrack->Pt();
 379             }
 380           }
 381         }
 382     }
 383
 384     hbtEvent->SetZDCEMEnergy(tTotalPt);
 385     //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
 386     //  if (esdtrack->GetTPCNcls() > 80)
 387     //    if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
 388     //      if (esdtrack->GetConstrainedParam())
 389     //        if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
 390     //          if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
 391     //            if (esdtrack->GetSign() > 0)
 392     //              tNormMultPos++;
 393     //            else if (esdtrack->GetSign() < 0)
 394     //              tNormMultNeg--;
 395     //          }
 396
 397     // If reading ITS-only tracks, reject all with TPC
 398     if (fTrackType == kITSOnly) {
 399       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
 400       if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
 401       if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
 402       UChar_t iclm = esdtrack->GetITSClusterMap();
 403       Int_t incls = 0;
 404       for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
 405       if (incls<=3) {
 406         if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
 407         continue;
 408       }
 409     }
 410
 411
 412
 413     AliFemtoTrack* trackCopy = new AliFemtoTrack();
 414     trackCopy->SetCharge((short)esdtrack->GetSign());
 415
 416     //in aliroot we have AliPID
 417     //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 418     //we use only 5 first
 419     double esdpid[5];
 420     esdtrack->GetESDpid(esdpid);
 421     trackCopy->SetPidProbElectron(esdpid[0]);
 422     trackCopy->SetPidProbMuon(esdpid[1]);
 423     trackCopy->SetPidProbPion(esdpid[2]);
 424     trackCopy->SetPidProbKaon(esdpid[3]);
 425     trackCopy->SetPidProbProton(esdpid[4]);
 426
 427     esdpid[0] = -100000.0;
 428     esdpid[1] = -100000.0;
 429     esdpid[2] = -100000.0;
 430     esdpid[3] = -100000.0;
 431     esdpid[4] = -100000.0;
 432
 433     double tTOF = 0.0;
 434
 435     if (esdtrack->GetStatus()&AliESDtrack::kTOFout&AliESDtrack::kTIME) {
 436       tTOF = esdtrack->GetTOFsignal();
 437       esdtrack->GetIntegratedTimes(esdpid);
 438     }
 439
 440     trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
 441
 442
 443     //////  TPC ////////////////////////////////////////////
 444     float nsigmaTPCK=-1000.;
 445     float nsigmaTPCPi=-1000.;
 446     float nsigmaTPCP=-1000.;
 447
 448     if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
 449       nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
 450       nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
 451       nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
 452
 453     }
 454     trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
 455     trackCopy->SetNSigmaTPCK(nsigmaTPCK);
 456     trackCopy->SetNSigmaTPCP(nsigmaTPCP);
 457
 458     ///// TOF ///////////////////////////////////////////////
 459
 460     float vp=-1000.;
 461     float nsigmaTOFPi=-1000.;
 462     float nsigmaTOFK=-1000.;
 463     float nsigmaTOFP=-1000.;
 464
 465     if (// (esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 466         (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 467         (esdtrack->GetStatus()&AliESDtrack::kTIME))
 468           {
 469
 470             //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 471             //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 472             //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
 473             // collect info from ESDpid class
 474
 475       if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFout&AliESDtrack::kTIME))) {
 476
 477
 478               double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
 479
 480               nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
 481               nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
 482               nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
 483
 484               Double_t len=esdtrack->GetIntegratedLength();
 485               Double_t tof=esdtrack->GetTOFsignal();
 486               if(tof > 0.) vp=len/tof/0.03;
 487             }
 488           }
 489
 490     trackCopy->SetVTOF(vp);
 491     trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
 492     trackCopy->SetNSigmaTOFK(nsigmaTOFK);
 493     trackCopy->SetNSigmaTOFP(nsigmaTOFP);
 494
 495
 496     double pxyz[3];
 497     double rxyz[3];
 498     double impact[2];
 499     double covimpact[3];
 500
 501     if (fUseTPCOnly) {
 502       if (!esdtrack->GetTPCInnerParam()) {
 503         cout << "No TPC inner param !" << endl;
 504         delete trackCopy;
 505         continue;
 506       }
 507
 508       AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 509       param->GetXYZ(rxyz);
 510       param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 511       param->GetPxPyPz(pxyz);//reading noconstarined momentum
 512
 513       if (fReadInner == true) {
 514         AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 515         tInfo->SetPDGPid(211);
 516         tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 517         tInfo->SetMass(0.13957);
 518         //        tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 519         //        tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 520         tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 521         trackCopy->SetHiddenInfo(tInfo);
 522       }
 523
 524       if (fRotateToEventPlane) {
 525         double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 526         double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 527
 528         pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 529         pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 530       }
 531
 532       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 533       if (v.Mag() < 0.0001) {
 534         //        cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 535         delete trackCopy;
 536         continue;
 537       }
 538       trackCopy->SetP(v);//setting momentum
 539       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 540
 541       const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 542       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 543       //setting helix I do not if it is ok
 544       AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 545       trackCopy->SetHelix(helix);
 546
 547       //some stuff which could be useful
 548       trackCopy->SetImpactD(impact[0]);
 549       trackCopy->SetImpactZ(impact[1]);
 550       trackCopy->SetCdd(covimpact[0]);
 551       trackCopy->SetCdz(covimpact[1]);
 552       trackCopy->SetCzz(covimpact[2]);
 553       trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 554
 555       delete param;
 556     }
 557     else {
 558       if (fReadInner == true) {
 559
 560         if (esdtrack->GetTPCInnerParam()) {
 561           AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 562           trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum());
 563           param->GetXYZ(rxyz);
 564           //        param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 565           param->GetPxPyPz(pxyz);//reading noconstarined momentum
 566           delete param;
 567
 568           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 569           tInfo->SetPDGPid(211);
 570           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 571           tInfo->SetMass(0.13957);
 572           //        tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 573           //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 574           tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 575           trackCopy->SetHiddenInfo(tInfo);
 576         }
 577       }
 578
 579
 580       if(fTrackType == kGlobal) {
 581         if (fConstrained==true)
 582           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 583         else
 584           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 585       }
 586       else if (fTrackType == kTPCOnly) {
 587         if (esdtrack->GetTPCInnerParam())
 588           esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
 589         else {
 590           delete trackCopy;
 591           continue;
 592         }
 593       }
 594       else if (fTrackType == kITSOnly) {
 595         if (fConstrained==true)
 596           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 597         else
 598           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 599       }
 600
 601       if (fRotateToEventPlane) {
 602         double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 603         double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 604
 605         pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 606         pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 607       }
 608
 609       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 610       if (v.Mag() < 0.0001) {
 611         //        cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 612         delete trackCopy;
 613         continue;
 614       }
 615
 616       trackCopy->SetP(v);//setting momentum
 617       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 618       const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 619       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 620       //setting helix I do not if it is ok
 621       AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 622       trackCopy->SetHelix(helix);
 623
 624       //some stuff which could be useful
 625       float imp[2];
 626       float cim[3];
 627       esdtrack->GetImpactParameters(imp,cim);
 628
 629       impact[0] = imp[0];
 630       impact[1] = imp[1];
 631       covimpact[0] = cim[0];
 632       covimpact[1] = cim[1];
 633       covimpact[2] = cim[2];
 634
 635       trackCopy->SetImpactD(impact[0]);
 636       trackCopy->SetImpactZ(impact[1]);
 637       trackCopy->SetCdd(covimpact[0]);
 638       trackCopy->SetCdz(covimpact[1]);
 639       trackCopy->SetCzz(covimpact[2]);
 640       trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 641     }
 642
 643     trackCopy->SetTrackId(esdtrack->GetID());
 644     trackCopy->SetFlags(esdtrack->GetStatus());
 645     trackCopy->SetLabel(esdtrack->GetLabel());
 646
 647     trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 648     trackCopy->SetITSncls(esdtrack->GetNcls(0));
 649     trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 650     trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 651     trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 652     trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 653     trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 654     trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
 655
 656     trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 657     trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 658
 659     double xtpc[3];
 660     esdtrack->GetInnerXYZ(xtpc);
 661     xtpc[2] -= fV1[2];
 662     trackCopy->SetNominalTPCEntrancePoint(xtpc);
 663
 664     esdtrack->GetOuterXYZ(xtpc);
 665     xtpc[2] -= fV1[2];
 666     trackCopy->SetNominalTPCExitPoint(xtpc);
 667
 668     int indexes[3];
 669     for (int ik=0; ik<3; ik++) {
 670       indexes[ik] = esdtrack->GetKinkIndex(ik);
 671     }
 672     trackCopy->SetKinkIndexes(indexes);
 673
 674     for (int ii=0; ii<6; ii++){
 675       trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii));
 676     }
 677
 678
 679     // Fill the hidden information with the simulated data
 680     if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
 681       TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 682
 683       // Check the mother information
 684
 685       // Using the new way of storing the freeze-out information
 686       // Final state particle is stored twice on the stack
 687       // one copy (mother) is stored with original freeze-out information
 688       //   and is not tracked
 689       // the other one (daughter) is stored with primary vertex position
 690       //   and is tracked
 691
 692       // Freeze-out coordinates
 693       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 694       fpx = tPart->Vx() - fV1[0];
 695       fpy = tPart->Vy() - fV1[1];
 696       fpz = tPart->Vz() - fV1[2];
 697       fpt = tPart->T();
 698
 699       AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 700       tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 701       trackCopy->SetGlobalEmissionPoint(fpx, fpy, fpz);
 702
 703       fpx *= 1e13;
 704       fpy *= 1e13;
 705       fpz *= 1e13;
 706       fpt *= 1e13;
 707
 708
 709       if (fOnlyPrimaries && !fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){
 710         delete trackCopy;
 711         continue;
 712       }
 713
 714       // fillDCA
 715
 716       if (TMath::Abs(impact[0]) > 0.001) {
 717         if (fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){
 718           trackCopy->SetImpactDprim(impact[0]);
 719           //cout << "prim" << endl;
 720
 721         }
 722         else if (fStack->IsSecondaryFromWeakDecay(TMath::Abs(esdtrack->GetLabel()))) {
 723           trackCopy->SetImpactDweak(impact[0]);
 724           //cout << "wea" << endl;
 725         }
 726         else if (fStack->IsSecondaryFromMaterial(TMath::Abs(esdtrack->GetLabel()))) {
 727           trackCopy->SetImpactDmat(impact[0]);
 728           //cout << "mat" << endl;
 729         }
 730       }
 731       //  end fillDCA
 732
 733       //      cout << "Looking for mother ids " << endl;
 734       if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 735         //      cout << "Got mother id" << endl;
 736         TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 737         // Check if this is the same particle stored twice on the stack
 738         if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 739           // It is the same particle
 740           // Read in the original freeze-out information
 741           // and convert it from to [fm]
 742
 743           // EPOS style
 744           fpx = mother->Vx()*1e13*0.197327;
 745           fpy = mother->Vy()*1e13*0.197327;
 746           fpz = mother->Vz()*1e13*0.197327;
 747           fpt = mother->T() *1e13*0.197327;
 748
 749
 750           // Therminator style
 751 //          fpx = mother->Vx()*1e13;
 752 //          fpy = mother->Vy()*1e13;
 753 //          fpz = mother->Vz()*1e13;
 754 //          fpt = mother->T() *1e13*3e10;
 755
 756         }
 757       }
 758
 759       if (fRotateToEventPlane) {
 760         double tPhi = TMath::ATan2(fpy, fpx);
 761         double tRad = TMath::Hypot(fpx, fpy);
 762
 763         fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 764         fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 765       }
 766
 767       tInfo->SetPDGPid(tPart->GetPdgCode());
 768       trackCopy->SetPDGPid(tPart->GetPdgCode());
 769
 770       if (fRotateToEventPlane) {
 771         double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 772         double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 773
 774         tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 775                                tRad*TMath::Sin(tPhi - tReactionPlane),
 776                                tPart->Pz());
 777         trackCopy->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 778                                    tRad*TMath::Sin(tPhi - tReactionPlane),
 779                                    tPart->Pz());
 780       }
 781       else
 782       {
 783         tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 784         trackCopy->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 785       }
 786       Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 787                         tPart->Px()*tPart->Px() -
 788                         tPart->Py()*tPart->Py() -
 789                         tPart->Pz()*tPart->Pz());
 790       if (mass2>0.0)
 791       {
 792         tInfo->SetMass(TMath::Sqrt(mass2));
 793         trackCopy->SetMass(TMath::Sqrt(mass2));
 794       }
 795       else
 796       {
 797         tInfo->SetMass(0.0);
 798         trackCopy->SetMass(0.0);
 799       }
 800
 801       tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 802       trackCopy->SetEmissionPoint(fpx, fpy, fpz, fpt);
 803       trackCopy->SetHiddenInfo(tInfo);
 804     }
 805     else {
 806       AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 807       tInfo->SetMass(0.0);
 808       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 809       fpx = fV1[0]*1e13;
 810       fpy = fV1[1]*1e13;
 811       fpz = fV1[2]*1e13;
 812       fpt = 0.0;
 813       tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 814
 815       tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
 816
 817       trackCopy->SetHiddenInfo(tInfo);
 818     }
 819     //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl;
 820
 821     //decision if we want this track
 822     //if we using diffrent labels we want that this label was use for first time
 823     //if we use hidden info we want to have match between sim data and ESD
 824
 825     bool trackAccept = true;
 826     if (isKaonAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 321) {
 827       trackAccept = false;
 828     }
 829     if (isProtonAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 2212) {
 830       trackAccept = false;
 831     }
 832     if (isPionAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 211) {
 833       trackAccept = false;
 834     }
 835
 836     if (tGoodMomentum==true  && trackAccept == true)
 837     {
 838
 839       hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 840       realnofTracks++;//real number of tracks
 841       //          delete trackCopy;
 842     }
 843     else
 844     {
 845       delete  trackCopy;
 846     }
 847
 848   }
 849
 850   delete [] motherids;
 851
 852   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 853
 854   AliCentrality *cent = fEvent->GetCentrality();
 855   if (cent) {
 856     hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 857     hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
 858     hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
 859     hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
 860     hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
 861     hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
 862     hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
 863     hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
 864     hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 865     hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
 866     //    hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 867     hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 868     hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
 869
 870     if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
 871   }
 872
 873   if (fEstEventMult == kGlobalCount)
 874     hbtEvent->SetNormalizedMult(tNormMult);
 875   else if (fEstEventMult == kReferenceITSTPC)
 876     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.2));
 877   else if(fEstEventMult == kReferenceITSSA)
 878     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.2));
 879   else if(fEstEventMult == kReferenceTracklets)
 880     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.2));
 881   else if (fEstEventMult == kSPDLayer1)
 882     hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 883   else if (fEstEventMult == kVZERO)
 884   {
 885     Float_t multV0 = 0;
 886     for (Int_t i=0; i<64; i++)
 887       multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
 888     hbtEvent->SetNormalizedMult(multV0);
 889   }
 890   else if (fEstEventMult == kCentrality) {
 891     // centrality between 0 (central) and 1 (very peripheral)
 892
 893     if (cent) {
 894       if (cent->GetCentralityPercentile("V0M") < 0.00001)
 895         hbtEvent->SetNormalizedMult(-1);
 896       else
 897         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
 898       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 899                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
 900     }
 901   }
 902   else if (fEstEventMult == kCentralityV0A) {
 903     // centrality between 0 (central) and 1 (very peripheral)
 904
 905     if (cent) {
 906       if (cent->GetCentralityPercentile("V0A") < 0.00001)
 907         hbtEvent->SetNormalizedMult(-1);
 908       else
 909         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A")));
 910       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 911                              10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A")));
 912     }
 913   }
 914   else if (fEstEventMult == kCentralityV0C) {
 915     // centrality between 0 (central) and 1 (very peripheral)
 916
 917     if (cent) {
 918       if (cent->GetCentralityPercentile("V0C") < 0.00001)
 919         hbtEvent->SetNormalizedMult(-1);
 920       else
 921         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C")));
 922       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 923                              10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C")));
 924     }
 925   }
 926   else if (fEstEventMult == kCentralityZNA) {
 927     // centrality between 0 (central) and 1 (very peripheral)
 928
 929     if (cent) {
 930       if (cent->GetCentralityPercentile("ZNA") < 0.00001)
 931         hbtEvent->SetNormalizedMult(-1);
 932       else
 933         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA")));
 934       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 935                              10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA")));
 936     }
 937   }
 938   else if (fEstEventMult == kCentralityZNC) {
 939     // centrality between 0 (central) and 1 (very peripheral)
 940
 941     if (cent) {
 942       if (cent->GetCentralityPercentile("ZNC") < 0.00001)
 943         hbtEvent->SetNormalizedMult(-1);
 944       else
 945         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC")));
 946       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 947                              10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC")));
 948     }
 949   }
 950   else if (fEstEventMult == kCentralityCL1) {
 951     // centrality between 0 (central) and 1 (very peripheral)
 952
 953     if (cent) {
 954       if (cent->GetCentralityPercentile("CL1") < 0.00001)
 955         hbtEvent->SetNormalizedMult(-1);
 956       else
 957         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1")));
 958       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 959                              10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1")));
 960     }
 961   }
 962   else if (fEstEventMult == kCentralityCL0) {
 963     // centrality between 0 (central) and 1 (very peripheral)
 964
 965     if (cent) {
 966       if (cent->GetCentralityPercentile("CL0") < 0.00001)
 967         hbtEvent->SetNormalizedMult(-1);
 968       else
 969         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0")));
 970       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 971                              10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0")));
 972     }
 973   }
 974   else if (fEstEventMult == kCentralityTRK) {
 975     // centrality between 0 (central) and 1 (very peripheral)
 976
 977     if (cent) {
 978       if (cent->GetCentralityPercentile("TRK") < 0.00001)
 979         hbtEvent->SetNormalizedMult(-1);
 980       else
 981         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK")));
 982       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 983                              10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK")));
 984     }
 985   }
 986   else if (fEstEventMult == kCentralityTKL) {
 987     // centrality between 0 (central) and 1 (very peripheral)
 988
 989     if (cent) {
 990       if (cent->GetCentralityPercentile("TKL") < 0.00001)
 991         hbtEvent->SetNormalizedMult(-1);
 992       else
 993         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL")));
 994       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 995                              10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL")));
 996     }
 997   }
 998   else if (fEstEventMult == kCentralityNPA) {
 999     // centrality between 0 (central) and 1 (very peripheral)
1000
1001     if (cent) {
1002       if (cent->GetCentralityPercentile("NPA") < 0.00001)
1003         hbtEvent->SetNormalizedMult(-1);
1004       else
1005         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA")));
1006       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
1007                              10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA")));
1008     }
1009   }
1010   else if (fEstEventMult == kCentralityCND) {
1011     // centrality between 0 (central) and 1 (very peripheral)
1012
1013     if (cent) {
1014       if (cent->GetCentralityPercentile("CND") < 0.00001)
1015         hbtEvent->SetNormalizedMult(-1);
1016       else
1017         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND")));
1018       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
1019                              10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND")));
1020     }
1021   }
1022   else if (fEstEventMult == kCentralityFMD) {
1023     // centrality between 0 (central) and 1 (very peripheral)
1024
1025     if (cent) {
1026       if (cent->GetCentralityPercentile("FMD") < 0.00001)
1027         hbtEvent->SetNormalizedMult(-1);
1028       else
1029         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD")));
1030       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
1031                              10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD")));
1032     }
1033   }
1034
1035
1036   if (tNormMultPos > tNormMultNeg)
1037     hbtEvent->SetZDCParticipants(tNormMultPos);
1038   else
1039     hbtEvent->SetZDCParticipants(tNormMultNeg);
1040
1041   AliEventplane* ep = fEvent->GetEventplane();
1042   if (ep) {
1043     hbtEvent->SetEP(ep);
1044     hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
1045   }
1046
1047   if(fReadV0)
1048   {
1049     for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
1050       AliESDv0* esdv0 = fEvent->GetV0(i);
1051       if (!esdv0) continue;
1052       //if(esdv0->GetNDaughters()>2) continue;
1053       //if(esdv0->GetNProngs()>2) continue;
1054       if(esdv0->Charge()!=0) continue;
1055       AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex());
1056       if(!trackPos) continue;
1057       AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex());
1058       if(!trackNeg) continue;
1059       if(trackPos->Charge()==trackNeg->Charge()) continue;
1060       AliFemtoV0* trackCopyV0 = new AliFemtoV0();
1061       CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent);
1062       hbtEvent->V0Collection()->push_back(trackCopyV0);
1063       //cout<<"Pushback v0 to v0collection"<<endl;
1064     }
1065   }
1066
1067
1068   fCurEvent++;
1069   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
1070   return hbtEvent;
1071 }
1072 //___________________
1073 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
1074 {
1075   // The chain loads the ESD for us
1076   // You must provide the address where it can be found
1077   fEvent = aESD;
1078 }
1079 //___________________
1080 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
1081 {
1082   // The chain loads the stack for us
1083   // You must provide the address where it can be found
1084   fStack = aStack;
1085 }
1086 //___________________
1087 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
1088 {
1089   // The chain loads the generator event header for us
1090   // You must provide the address where it can be found
1091   fGenHeader = aGenHeader;
1092 }
1093 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
1094 {
1095   fRotateToEventPlane=dorotate;
1096 }
1097 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
1098 {
1099   // Calculates the number of sigma to the vertex.
1100
1101   Float_t b[2];
1102   Float_t bRes[2];
1103   Float_t bCov[3];
1104
1105   b[0] = impact[0];
1106   b[1] = impact[1];
1107   bCov[0] = covar[0];
1108   bCov[1] = covar[1];
1109   bCov[2] = covar[2];
1110
1111   bRes[0] = TMath::Sqrt(bCov[0]);
1112   bRes[1] = TMath::Sqrt(bCov[2]);
1113
1114   // -----------------------------------
1115   // How to get to a n-sigma cut?
1116   //
1117   // The accumulated statistics from 0 to d is
1118   //
1119   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1120   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1121   //
1122   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
1123   // Can this be expressed in a different way?
1124
1125   if (bRes[0] == 0 || bRes[1] ==0)
1126     return -1;
1127
1128   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1129
1130   // stupid rounding problem screws up everything:
1131   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1132   if (TMath::Exp(-d * d / 2) < 1e-10)
1133     return 1000;
1134
1135   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1136   return d;
1137 }
1138
1139 void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)
1140 {
1141   fTrackType = aType;
1142 }
1143
1144 void AliFemtoEventReaderESDChainKine::SetReadV0(bool a)
1145 {
1146   fReadV0 = a;
1147 }
1148
1149 void AliFemtoEventReaderESDChainKine::SetKaonAnalysis(Bool_t a)
1150 {
1151   isKaonAnalysis = a;
1152 }
1153
1154 void AliFemtoEventReaderESDChainKine::SetProtonAnalysis(Bool_t a)
1155 {
1156   isProtonAnalysis = a;
1157 }
1158
1159 void AliFemtoEventReaderESDChainKine::SetPionAnalysis(Bool_t a)
1160 {
1161   isPionAnalysis = a;
1162 }
1163
1164 void AliFemtoEventReaderESDChainKine::SetOnlyPrimaries(Bool_t a)
1165 {
1166   fOnlyPrimaries = a;
1167 }
1168
1169 void AliFemtoEventReaderESDChainKine::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent)
1170 {
1171   double fPrimaryVtxPosition[3];
1172   tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition);
1173   tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2]));
1174
1175   //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego
1176   //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X());
1177   //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y());
1178   //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0
1179   //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z());
1180   //tFemtoV0->SetdecayVertexV0(decayvertex);
1181   tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters());
1182   tFemtoV0->SetmomV0X(tESDv0->Px());
1183   tFemtoV0->SetmomV0Y(tESDv0->Py());
1184   tFemtoV0->SetmomV0Z(tESDv0->Pz());
1185   AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz());
1186   tFemtoV0->SetmomV0(momv0);
1187   tFemtoV0->SetalphaV0(tESDv0->AlphaV0());
1188   tFemtoV0->SetptArmV0(tESDv0->PtArmV0());
1189   //tFemtoV0->SeteLambda(tESDv0->ELambda());
1190   //tFemtoV0->SeteK0Short(tESDv0->EK0Short());
1191   //tFemtoV0->SetePosProton(tESDv0->EPosProton());
1192   //tFemtoV0->SeteNegProton(tESDv0->ENegProton());
1193   tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2));
1194   tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4));
1195   tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2));
1196   //tFemtoV0->SetrapLambda(tESDv0->RapLambda());
1197   //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short());
1198   tFemtoV0->SetptV0(tESDv0->Pt());
1199   tFemtoV0->SetptotV0(tESDv0->P());
1200   tFemtoV0->SetEtaV0(tESDv0->Eta());
1201   tFemtoV0->SetPhiV0(tESDv0->Phi());
1202   tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2]));
1203   tFemtoV0->SetYV0(tESDv0->Y());
1204
1205   AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0);
1206   AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1);
1207
1208   if(trackpos && trackneg)
1209   {
1210     tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1211     tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1212     double MomPos[3];
1213     trackpos->PxPyPz(MomPos);
1214     tFemtoV0->SetmomPosX(MomPos[0]);
1215     tFemtoV0->SetmomPosY(MomPos[1]);
1216     tFemtoV0->SetmomPosZ(MomPos[2]);
1217     AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]);
1218     tFemtoV0->SetmomPos(mompos);
1219
1220     double MomNeg[3];
1221     trackneg->PxPyPz(MomNeg);
1222     tFemtoV0->SetmomNegX(MomNeg[0]);
1223     tFemtoV0->SetmomNegY(MomNeg[1]);
1224     tFemtoV0->SetmomNegZ(MomNeg[2]);
1225     AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]);
1226     tFemtoV0->SetmomNeg(momneg);
1227
1228     tFemtoV0->SetptPos(trackpos->Pt());
1229     tFemtoV0->SetptotPos(trackpos->P());
1230     tFemtoV0->SetptNeg(trackneg->Pt());
1231     tFemtoV0->SetptotNeg(trackneg->P());
1232
1233     tFemtoV0->SetidNeg(trackneg->GetID());
1234     //cout<<"tESDv0->GetNegID(): "<<tESDv0->GetNegID()<<endl;
1235     //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1236     tFemtoV0->SetidPos(trackpos->GetID());
1237
1238     tFemtoV0->SetEtaPos(trackpos->Eta());
1239     tFemtoV0->SetEtaNeg(trackneg->Eta());
1240
1241     tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos()
1242     tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg()
1243     tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1244     tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1245     tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1246     tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1247     tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1248     tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1249     tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls());
1250     tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls());
1251     tFemtoV0->SetStatusPos(trackpos->GetStatus());
1252     tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1253
1254     float bfield = 5*fMagFieldSign;
1255     float globalPositionsAtRadiiPos[9][3];
1256     GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1257     double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1258     double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1259
1260     float globalPositionsAtRadiiNeg[9][3];
1261     GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1262     double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1263     double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1264
1265     AliFemtoThreeVector tmpVec;
1266     tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
1267     tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1268
1269     tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
1270     tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1271
1272     tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
1273     tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1274
1275     tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
1276     tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1277
1278     AliFemtoThreeVector vecTpcPos[9];
1279     AliFemtoThreeVector vecTpcNeg[9];
1280     for(int i=0;i<9;i++)
1281     {
1282       vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1283       vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1284     }
1285     tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1286     tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
1287
1288     tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum();
1289     tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum();
1290
1291     tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1292     tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1293
1294
1295     if (fESDpid) {
1296       tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1297       tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1298       tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1299       tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1300       tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1301       tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1302     }
1303     else {
1304       tFemtoV0->SetPosNSigmaTPCK(-1000);
1305       tFemtoV0->SetNegNSigmaTPCK(-1000);
1306       tFemtoV0->SetPosNSigmaTPCP(-1000);
1307       tFemtoV0->SetNegNSigmaTPCP(-1000);
1308       tFemtoV0->SetPosNSigmaTPCPi(-1000);
1309       tFemtoV0->SetNegNSigmaTPCPi(-1000);
1310     }
1311
1312     if(// (tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 ||
1313        (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0)
1314     {
1315       if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
1316          (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0)
1317             {
1318               tFemtoV0->SetPosNSigmaTOFK(-1000);
1319               tFemtoV0->SetNegNSigmaTOFK(-1000);
1320               tFemtoV0->SetPosNSigmaTOFP(-1000);
1321               tFemtoV0->SetNegNSigmaTOFP(-1000);
1322               tFemtoV0->SetPosNSigmaTOFPi(-1000);
1323               tFemtoV0->SetNegNSigmaTOFPi(-1000);
1324             }
1325     }
1326     else
1327     {
1328       if (fESDpid) {
1329         tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1330         tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1331         tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1332         tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1333         tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1334         tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1335       }
1336       else {
1337         tFemtoV0->SetPosNSigmaTOFK(-1000);
1338         tFemtoV0->SetNegNSigmaTOFK(-1000);
1339         tFemtoV0->SetPosNSigmaTOFP(-1000);
1340         tFemtoV0->SetNegNSigmaTOFP(-1000);
1341         tFemtoV0->SetPosNSigmaTOFPi(-1000);
1342         tFemtoV0->SetNegNSigmaTOFPi(-1000);
1343       }
1344     }
1345
1346     if ((TMath::Abs(trackpos->GetLabel()) < fStack->GetNtrack()) && (TMath::Abs(trackneg->GetLabel()) < fStack->GetNtrack())) {
1347
1348       //cout<<"tESDv0->GetPdgCode(): "<<tESDv0->GetPdgCode()<<endl;
1349       // cout<<"Labels: "<<trackpos->GetLabel()<<" "<<trackneg->GetLabel()<<endl;
1350       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
1351       //TParticle *tPart = fStack->Particle(tESDv0->GetLabel()); //zle
1352
1353       int labelpos = TMath::Abs(trackpos->GetLabel());
1354       int labelneg = TMath::Abs(trackneg->GetLabel());
1355       TParticle *tPartPos = fStack->Particle(labelpos);
1356       TParticle *tPartNeg = fStack->Particle(labelneg);
1357
1358
1359       double impactpos[2];
1360       double impactneg[2];
1361       double covimpact[3];
1362
1363       AliExternalTrackParam *parampos = new AliExternalTrackParam(*trackpos->GetTPCInnerParam());
1364       parampos->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactpos, covimpact);
1365       AliExternalTrackParam *paramneg = new AliExternalTrackParam(*trackneg->GetTPCInnerParam());
1366       paramneg->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactneg, covimpact);
1367
1368
1369       // fillDCA
1370       if (TMath::Abs(impactpos[0]) > 0.001) {
1371         if (fStack->IsPhysicalPrimary(labelpos)){
1372           tFemtoV0->SetImpactDprimPos(impactpos[0]);
1373         }
1374         else if (fStack->IsSecondaryFromWeakDecay(labelpos)) {
1375           tFemtoV0->SetImpactDweakPos(impactpos[0]);
1376           //cout << "wea" << endl;
1377         }
1378         else if (fStack->IsSecondaryFromMaterial(labelpos)) {
1379           tFemtoV0->SetImpactDmatPos(impactpos[0]);
1380           //cout << "mat" << endl;
1381         }
1382       }
1383       if (TMath::Abs(impactneg[0]) > 0.001) {
1384         if (fStack->IsPhysicalPrimary(labelneg)){
1385           tFemtoV0->SetImpactDprimNeg(impactneg[0]);
1386           //cout << "prim" << endl;
1387         }
1388         else if (fStack->IsSecondaryFromWeakDecay(labelneg)) {
1389           tFemtoV0->SetImpactDweakNeg(impactneg[0]);
1390           //cout << "wea" << endl;
1391         }
1392         else if (fStack->IsSecondaryFromMaterial(labelneg)) {
1393           tFemtoV0->SetImpactDmatNeg(impactneg[0]);
1394           //cout << "mat" << endl;
1395         }
1396
1397       }
1398       //  end fillDCA
1399
1400
1401       //tInfo->SetPDGPid();
1402       //tInfo->SetMass();
1403       //tInfo->SetTrueMomentum();
1404       //tInfo->SetEmissionPoint();
1405
1406       // Freeze-out coordinates
1407       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
1408
1409       fpx = tPartPos->Vx() - fPrimaryVtxPosition[0];
1410       fpy = tPartPos->Vy() - fPrimaryVtxPosition[1];
1411       fpz = tPartPos->Vz() - fPrimaryVtxPosition[2];
1412       fpt = tPartPos->T();
1413
1414       fpx *= 1e13;
1415       fpy *= 1e13;
1416       fpz *= 1e13;
1417       fpt *= 1e13;
1418
1419       tInfo->SetPDGPidPos(tPartPos->GetPdgCode());
1420       tInfo->SetMassPos(tPartPos->GetMass());
1421       tInfo->SetTrueMomentumPos(tPartPos->Px(),tPartPos->Py(),tPartPos->Pz());
1422       tInfo->SetEmissionPointPos(fpx,fpy,fpz,fpt);
1423
1424       fpx = tPartNeg->Vx() - fPrimaryVtxPosition[0];
1425       fpy = tPartNeg->Vy() - fPrimaryVtxPosition[1];
1426       fpz = tPartNeg->Vz() - fPrimaryVtxPosition[2];
1427       fpt = tPartNeg->T();
1428
1429       fpx *= 1e13;
1430       fpy *= 1e13;
1431       fpz *= 1e13;
1432       fpt *= 1e13;
1433
1434       tInfo->SetPDGPidNeg(tPartNeg->GetPdgCode());
1435       tInfo->SetMassNeg(tPartNeg->GetMass());
1436       tInfo->SetTrueMomentumNeg(tPartNeg->Px(),tPartNeg->Py(),tPartNeg->Pz());
1437       tInfo->SetEmissionPointNeg(fpx,fpy,fpz,fpt);
1438
1439       tFemtoV0->SetHiddenInfo(tInfo);
1440     }
1441   }
1442   else
1443   {
1444     tFemtoV0->SetStatusPos(999);
1445     tFemtoV0->SetStatusNeg(999);
1446   }
1447   tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus());
1448 }
1449
1450 void AliFemtoEventReaderESDChainKine::SetMagneticFieldSign(int s)
1451 {
1452   if(s>0)
1453     fMagFieldSign = 1;
1454   else if(s<0)
1455     fMagFieldSign = -1;
1456   else
1457     fMagFieldSign = 0;
1458 }
1459
1460 void AliFemtoEventReaderESDChainKine::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1461 {
1462   // Gets the global position of the track at nine different radii in the TPC
1463   // track is the track you want to propagate
1464   // bfield is the magnetic field of your event
1465   // globalPositionsAtRadii is the array of global positions in the radii and xyz
1466
1467   // Initialize the array to something indicating there was no propagation
1468   for(Int_t i=0;i<9;i++){
1469     for(Int_t j=0;j<3;j++){
1470       globalPositionsAtRadii[i][j]=-9999.;
1471     }
1472   }
1473
1474   // Make a copy of the track to not change parameters of the track
1475   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1476   //printf("\nAfter CopyFromVTrack\n");
1477   //etp.Print();
1478
1479   // The global position of the the track
1480   Double_t xyz[3]={-9999.,-9999.,-9999.};
1481
1482   // Counter for which radius we want
1483   Int_t iR=0;
1484   // The radii at which we get the global positions
1485   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1486   Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1487   // The global radius we are at
1488   Float_t globalRadius=0;
1489
1490   // Propagation is done in local x of the track
1491   for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1492     // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1493     // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1494     // If the track's momentum is smaller than infinite, it will develop a y-component, which
1495     // adds to the global radius
1496
1497     // Stop if the propagation was not succesful. This can happen for low pt tracks
1498     // that don't reach outer radii
1499     if(!etp.PropagateTo(x,bfield))break;
1500     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1501     globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1502
1503     // Roughly reached the radius we want
1504     if(globalRadius > Rwanted[iR]){
1505
1506       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1507       while (globalRadius>Rwanted[iR]){
1508         x-=.1;
1509         //      printf("propagating to x %5.2f\n",x);
1510         if(!etp.PropagateTo(x,bfield))break;
1511         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1512         globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1513       }
1514       //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]);
1515       globalPositionsAtRadii[iR][0]=xyz[0];
1516       globalPositionsAtRadii[iR][1]=xyz[1];
1517       globalPositionsAtRadii[iR][2]=xyz[2];
1518       // Indicate we want the next radius
1519       iR+=1;
1520     }
1521     if(iR>=8){
1522       // TPC edge reached
1523       return;
1524     }
1525   }
1526 }