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