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(kV0Centrality),
  60   fRotateToEventPlane(0),
  61   fESDpid(0),
  62   fIsPidOwner(0)
  63
  64 {
  65   //constructor with 0 parameters , look at default settings
  66 }
  67
  68 //__________________
  69 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  70   AliFemtoEventReader(aReader),
  71   fFileName(" "),
  72   fConstrained(true),
  73   fReadInner(false),
  74   fUseTPCOnly(false),
  75   fNumberofEvent(0),
  76   fCurEvent(0),
  77   fCurFile(0),
  78   fEvent(0x0),
  79   fStack(0x0),
  80   fGenHeader(0x0),
  81   fTrackType(kGlobal),
  82   fEstEventMult(kV0Centrality),
  83   fRotateToEventPlane(0),
  84   fESDpid(0),
  85   fIsPidOwner(0)
  86 {
  87   // Copy constructor
  88   fConstrained = aReader.fConstrained;
  89   fReadInner = aReader.fReadInner;
  90   fUseTPCOnly = aReader.fUseTPCOnly;
  91   fNumberofEvent = aReader.fNumberofEvent;
  92   fCurEvent = aReader.fCurEvent;
  93   fCurFile = aReader.fCurFile;
  94   fEvent = new AliESDEvent();
  95   fStack = aReader.fStack;
  96   fTrackType = aReader.fTrackType;
  97   fEstEventMult = aReader.fEstEventMult;
  98   fRotateToEventPlane = aReader.fRotateToEventPlane;
  99   fESDpid = aReader.fESDpid;
 100   fIsPidOwner = aReader.fIsPidOwner;
 101 }
 102 //__________________
 103 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
 104 {
 105   //Destructor
 106   delete fEvent;
 107 }
 108
 109 //__________________
 110 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
 111 {
 112   // Assignment operator
 113   if (this == &aReader)
 114     return *this;
 115
 116   fConstrained = aReader.fConstrained;
 117   fUseTPCOnly = aReader.fUseTPCOnly;
 118   fNumberofEvent = aReader.fNumberofEvent;
 119   fCurEvent = aReader.fCurEvent;
 120   fCurFile = aReader.fCurFile;
 121   if (fEvent) delete fEvent;
 122   fEvent = new AliESDEvent();
 123   fStack = aReader.fStack;
 124   fGenHeader = aReader.fGenHeader;
 125   fRotateToEventPlane = aReader.fRotateToEventPlane;
 126   fESDpid = aReader.fESDpid;
 127   fIsPidOwner = aReader.fIsPidOwner;
 128
 129   return *this;
 130 }
 131 //__________________
 132 // Simple report
 133 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
 134 {
 135   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
 136   return temp;
 137 }
 138
 139 //__________________
 140 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
 141 {
 142   // Select whether to read constrained or not constrained momentum
 143   fConstrained=constrained;
 144 }
 145 //__________________
 146 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
 147 {
 148   // Check whether we read constrained or not constrained momentum
 149   return fConstrained;
 150 }
 151 //__________________
 152 void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
 153 {
 154   fReadInner=readinner;
 155 }
 156
 157 bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
 158 {
 159   return fReadInner;
 160 }
 161
 162 //__________________
 163 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
 164 {
 165   fUseTPCOnly=usetpconly;
 166 }
 167
 168 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
 169 {
 170   return fUseTPCOnly;
 171 }
 172 void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
 173 {
 174   fEstEventMult = aType;
 175 }
 176 //__________________
 177 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 178 {
 179   // Get the event, read all the relevant information
 180   // and fill the AliFemtoEvent class
 181   // Returns a valid AliFemtoEvent
 182   AliFemtoEvent *hbtEvent = 0;
 183   string tFriendFileName;
 184
 185   // Get the friend information
 186   cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
 187   //  fEvent->SetESDfriend(fEventFriend);
 188
 189   hbtEvent = new AliFemtoEvent;
 190   //setting basic things
 191   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 192   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 193   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 194   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 195   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 196   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 197   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 198   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 199   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 200   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 201   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 202   //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 203
 204   if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
 205       (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
 206       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
 207       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
 208     hbtEvent->SetTriggerCluster(1);
 209   else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
 210            (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
 211     hbtEvent->SetTriggerCluster(2);
 212   else
 213     hbtEvent->SetTriggerCluster(0);
 214
 215   //Vertex
 216   double fV1[3];
 217   double fVCov[6];
 218   if (fUseTPCOnly) {
 219     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 220     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 221     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 222       fVCov[4] = -1001.0;
 223   }
 224   else {
 225     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 226     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 227     if (!fEvent->GetPrimaryVertex()->GetStatus())
 228       fVCov[4] = -1001.0;
 229   }
 230
 231   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 232   hbtEvent->SetPrimVertPos(vertex);
 233   hbtEvent->SetPrimVertCov(fVCov);
 234
 235   Double_t tReactionPlane = 0;
 236
 237   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 238   if (!hdh) {
 239     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 240     if (cdh) {
 241       TList *tGenHeaders = cdh->GetHeaders();
 242       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 243         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 244         if (hdh) break;
 245       }
 246     }
 247   }
 248
 249   if (hdh)
 250     {
 251       tReactionPlane = hdh->ReactionPlaneAngle();
 252       cout << "Got reaction plane " << tReactionPlane << endl;
 253     }
 254
 255   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 256
 257   Int_t spdetaonecount = 0;
 258
 259   for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
 260     if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
 261       spdetaonecount++;
 262
 263   //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
 264   hbtEvent->SetSPDMult(spdetaonecount);
 265
 266   //starting to reading tracks
 267   int nofTracks=0;  //number of reconstructed tracks in event
 268   nofTracks=fEvent->GetNumberOfTracks();
 269   int realnofTracks=0;//number of track which we use ina analysis
 270
 271   int tNormMult = 0;
 272   int tNormMultPos = 0;
 273   int tNormMultNeg = 0;
 274
 275   Float_t tTotalPt = 0.0;
 276
 277   Float_t b[2];
 278   Float_t bCov[3];
 279
 280   Int_t tTracklet=0, tITSTPC=0, tITSPure=0;
 281
 282   fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
 283
 284   hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
 285   hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
 286   //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
 287   hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 288
 289   Int_t *motherids;
 290   if (fStack) {
 291     motherids = new Int_t[fStack->GetNtrack()];
 292     for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 293
 294     // Read in mother ids
 295     TParticle *motherpart;
 296     for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 297       motherpart = fStack->Particle(ip);
 298       if (motherpart->GetDaughter(0) > 0)
 299         motherids[motherpart->GetDaughter(0)] = ip;
 300       if (motherpart->GetDaughter(1) > 0)
 301         motherids[motherpart->GetDaughter(1)] = ip;
 302
 303 //     if (motherpart->GetPdgCode() == 211) {
 304 //       cout << "Mother " << ip << " has daughters "
 305 //         << motherpart->GetDaughter(0) << " "
 306 //         << motherpart->GetDaughter(1) << " "
 307 //         << motherpart->Vx() << " "
 308 //         << motherpart->Vy() << " "
 309 //         << motherpart->Vz() << " "
 310 //         << endl;
 311
 312 //     }
 313     }
 314   }
 315   else {
 316     printf ("No Stack ???\n");
 317     delete hbtEvent;
 318     return 0;
 319   }
 320
 321   for (int i=0;i<nofTracks;i++)
 322     {
 323       //      cout << "Reading track " << i << endl;
 324       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 325
 326
 327       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 328       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 329
 330
 331       if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
 332           (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
 333         if (esdtrack->GetTPCNcls() > 70)
 334           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 335             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 336               esdtrack->GetImpactParameters(b,bCov);
 337               if ((b[0]<0.2) && (b[1] < 0.25)) {
 338                 tNormMult++;
 339                 tTotalPt += esdtrack->Pt();
 340               }
 341             }
 342           }
 343       }
 344       else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
 345         if (esdtrack->GetTPCNcls() > 100)
 346           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 347             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 348               esdtrack->GetImpactParameters(b,bCov);
 349               if ((b[0]<2.4) && (b[1] < 3.2)) {
 350                 tNormMult++;
 351                 tTotalPt += esdtrack->Pt();
 352               }
 353             }
 354           }
 355       }
 356
 357       hbtEvent->SetZDCEMEnergy(tTotalPt);
 358       //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
 359       //        if (esdtrack->GetTPCNcls() > 80)
 360       //          if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
 361       //            if (esdtrack->GetConstrainedParam())
 362       //              if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
 363       //                if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
 364       //                  if (esdtrack->GetSign() > 0)
 365       //                    tNormMultPos++;
 366       //                  else if (esdtrack->GetSign() < 0)
 367       //                    tNormMultNeg--;
 368       //                }
 369
 370       // If reading ITS-only tracks, reject all with TPC
 371       if (fTrackType == kITSOnly) {
 372         if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
 373         if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
 374         if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
 375         UChar_t iclm = esdtrack->GetITSClusterMap();
 376         Int_t incls = 0;
 377         for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
 378         if (incls<=3) {
 379           if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
 380           continue;
 381         }
 382       }
 383
 384
 385
 386       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 387       trackCopy->SetCharge((short)esdtrack->GetSign());
 388
 389       //in aliroot we have AliPID
 390       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 391       //we use only 5 first
 392       double esdpid[5];
 393       esdtrack->GetESDpid(esdpid);
 394       trackCopy->SetPidProbElectron(esdpid[0]);
 395       trackCopy->SetPidProbMuon(esdpid[1]);
 396       trackCopy->SetPidProbPion(esdpid[2]);
 397       trackCopy->SetPidProbKaon(esdpid[3]);
 398       trackCopy->SetPidProbProton(esdpid[4]);
 399
 400       esdpid[0] = -100000.0;
 401       esdpid[1] = -100000.0;
 402       esdpid[2] = -100000.0;
 403       esdpid[3] = -100000.0;
 404       esdpid[4] = -100000.0;
 405
 406       double tTOF = 0.0;
 407
 408       if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
 409         tTOF = esdtrack->GetTOFsignal();
 410         esdtrack->GetIntegratedTimes(esdpid);
 411       }
 412
 413       trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
 414
 415
 416       //////  TPC ////////////////////////////////////////////
 417       float nsigmaTPCK=-1000.;
 418       float nsigmaTPCPi=-1000.;
 419       float nsigmaTPCP=-1000.;
 420
 421      if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
 422         nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
 423         nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
 424         nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
 425
 426       }
 427       trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
 428       trackCopy->SetNSigmaTPCK(nsigmaTPCK);
 429       trackCopy->SetNSigmaTPCP(nsigmaTPCP);
 430
 431       ///// TOF ///////////////////////////////////////////////
 432
 433       float vp=-1000.;
 434       float nsigmaTOFPi=-1000.;
 435       float nsigmaTOFK=-1000.;
 436       float nsigmaTOFP=-1000.;
 437
 438         if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 439             (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 440             (esdtrack->GetStatus()&AliESDtrack::kTIME) &&
 441             !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch))
 442           {
 443
 444             //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 445             //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 446             //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
 447             // collect info from ESDpid class
 448
 449               if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
 450
 451
 452               double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
 453
 454               nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
 455               nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
 456               nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
 457
 458               Double_t len=esdtrack->GetIntegratedLength();
 459               Double_t tof=esdtrack->GetTOFsignal();
 460               if(tof > 0.) vp=len/tof/0.03;
 461             }
 462           }
 463
 464       trackCopy->SetVTOF(vp);
 465       trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
 466       trackCopy->SetNSigmaTOFK(nsigmaTOFK);
 467       trackCopy->SetNSigmaTOFP(nsigmaTOFP);
 468
 469
 470       double pxyz[3];
 471       double rxyz[3];
 472       double impact[2];
 473       double covimpact[3];
 474
 475       if (fUseTPCOnly) {
 476         if (!esdtrack->GetTPCInnerParam()) {
 477           cout << "No TPC inner param !" << endl;
 478           delete trackCopy;
 479           continue;
 480         }
 481
 482         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 483         param->GetXYZ(rxyz);
 484         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 485         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 486
 487           if (fReadInner == true) {
 488             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 489             tInfo->SetPDGPid(211);
 490             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 491             tInfo->SetMass(0.13957);
 492             //    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 493             //    tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 494             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 495             trackCopy->SetHiddenInfo(tInfo);
 496           }
 497
 498         if (fRotateToEventPlane) {
 499           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 500           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 501
 502           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 503           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 504         }
 505
 506         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 507         if (v.Mag() < 0.0001) {
 508           //      cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 509           delete trackCopy;
 510           continue;
 511         }
 512         trackCopy->SetP(v);//setting momentum
 513         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 514
 515         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 516         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 517         //setting helix I do not if it is ok
 518         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 519         trackCopy->SetHelix(helix);
 520
 521         //some stuff which could be useful
 522         trackCopy->SetImpactD(impact[0]);
 523         trackCopy->SetImpactZ(impact[1]);
 524         trackCopy->SetCdd(covimpact[0]);
 525         trackCopy->SetCdz(covimpact[1]);
 526         trackCopy->SetCzz(covimpact[2]);
 527         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 528
 529         delete param;
 530       }
 531       else {
 532           if (fReadInner == true) {
 533
 534             if (esdtrack->GetTPCInnerParam()) {
 535               AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 536               param->GetXYZ(rxyz);
 537               //            param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 538               param->GetPxPyPz(pxyz);//reading noconstarined momentum
 539               delete param;
 540
 541               AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 542               tInfo->SetPDGPid(211);
 543               tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 544               tInfo->SetMass(0.13957);
 545               //            tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 546               //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 547               tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 548               trackCopy->SetHiddenInfo(tInfo);
 549             }
 550           }
 551
 552
 553         if(fTrackType == kGlobal) {
 554         if (fConstrained==true)
 555           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 556         else
 557           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 558         }
 559           else if (fTrackType == kTPCOnly) {
 560             if (esdtrack->GetTPCInnerParam())
 561               esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
 562             else {
 563               delete trackCopy;
 564               continue;
 565             }
 566           }
 567           else if (fTrackType == kITSOnly) {
 568             if (fConstrained==true)
 569               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 570             else
 571               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 572           }
 573
 574         if (fRotateToEventPlane) {
 575           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 576           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 577
 578           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 579           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 580         }
 581
 582         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 583         if (v.Mag() < 0.0001) {
 584           //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 585           delete trackCopy;
 586           continue;
 587         }
 588
 589         trackCopy->SetP(v);//setting momentum
 590         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 591         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 592         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 593         //setting helix I do not if it is ok
 594         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 595         trackCopy->SetHelix(helix);
 596
 597         //some stuff which could be useful
 598         float imp[2];
 599         float cim[3];
 600         esdtrack->GetImpactParameters(imp,cim);
 601
 602         impact[0] = imp[0];
 603         impact[1] = imp[1];
 604         covimpact[0] = cim[0];
 605         covimpact[1] = cim[1];
 606         covimpact[2] = cim[2];
 607
 608         trackCopy->SetImpactD(impact[0]);
 609         trackCopy->SetImpactZ(impact[1]);
 610         trackCopy->SetCdd(covimpact[0]);
 611         trackCopy->SetCdz(covimpact[1]);
 612         trackCopy->SetCzz(covimpact[2]);
 613         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 614       }
 615
 616       trackCopy->SetTrackId(esdtrack->GetID());
 617       trackCopy->SetFlags(esdtrack->GetStatus());
 618       trackCopy->SetLabel(esdtrack->GetLabel());
 619
 620       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 621       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 622       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 623       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 624       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 625       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 626       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 627
 628
 629       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 630       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 631
 632       double xtpc[3];
 633       esdtrack->GetInnerXYZ(xtpc);
 634       xtpc[2] -= fV1[2];
 635       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 636
 637       esdtrack->GetOuterXYZ(xtpc);
 638       xtpc[2] -= fV1[2];
 639       trackCopy->SetNominalTPCExitPoint(xtpc);
 640
 641       int indexes[3];
 642       for (int ik=0; ik<3; ik++) {
 643         indexes[ik] = esdtrack->GetKinkIndex(ik);
 644       }
 645       trackCopy->SetKinkIndexes(indexes);
 646
 647
 648       // Fill the hidden information with the simulated data
 649       if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
 650         TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 651
 652         // Check the mother information
 653
 654         // Using the new way of storing the freeze-out information
 655         // Final state particle is stored twice on the stack
 656         // one copy (mother) is stored with original freeze-out information
 657         //   and is not tracked
 658         // the other one (daughter) is stored with primary vertex position
 659         //   and is tracked
 660
 661         // Freeze-out coordinates
 662         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 663         fpx = tPart->Vx() - fV1[0];
 664         fpy = tPart->Vy() - fV1[1];
 665         fpz = tPart->Vz() - fV1[2];
 666         fpt = tPart->T();
 667
 668         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 669         tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 670
 671         fpx *= 1e13;
 672         fpy *= 1e13;
 673         fpz *= 1e13;
 674         fpt *= 1e13;
 675
 676         //      cout << "Looking for mother ids " << endl;
 677         if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 678           //    cout << "Got mother id" << endl;
 679           TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 680           // Check if this is the same particle stored twice on the stack
 681           if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 682             // It is the same particle
 683             // Read in the original freeze-out information
 684             // and convert it from to [fm]
 685
 686             // EPOS style
 687             fpx = mother->Vx()*1e13*0.197327;
 688             fpy = mother->Vy()*1e13*0.197327;
 689             fpz = mother->Vz()*1e13*0.197327;
 690             fpt = mother->T() *1e13*0.197327;
 691
 692
 693             // Therminator style
 694 //          fpx = mother->Vx()*1e13;
 695 //          fpy = mother->Vy()*1e13;
 696 //          fpz = mother->Vz()*1e13;
 697 //          fpt = mother->T() *1e13*3e10;
 698
 699           }
 700         }
 701
 702         if (fRotateToEventPlane) {
 703           double tPhi = TMath::ATan2(fpy, fpx);
 704           double tRad = TMath::Hypot(fpx, fpy);
 705
 706           fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 707           fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 708         }
 709
 710         tInfo->SetPDGPid(tPart->GetPdgCode());
 711
 712         if (fRotateToEventPlane) {
 713           double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 714           double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 715
 716           tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 717                                  tRad*TMath::Sin(tPhi - tReactionPlane),
 718                                  tPart->Pz());
 719         }
 720         else
 721           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 722         Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 723                           tPart->Px()*tPart->Px() -
 724                           tPart->Py()*tPart->Py() -
 725                           tPart->Pz()*tPart->Pz());
 726         if (mass2>0.0)
 727           tInfo->SetMass(TMath::Sqrt(mass2));
 728         else
 729           tInfo->SetMass(0.0);
 730
 731         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 732         trackCopy->SetHiddenInfo(tInfo);
 733       }
 734       else {
 735         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 736         tInfo->SetMass(0.0);
 737         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 738         fpx = fV1[0]*1e13;
 739         fpy = fV1[1]*1e13;
 740         fpz = fV1[2]*1e13;
 741         fpt = 0.0;
 742         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 743
 744         tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
 745
 746         trackCopy->SetHiddenInfo(tInfo);
 747       }
 748       //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl;
 749
 750       //decision if we want this track
 751       //if we using diffrent labels we want that this label was use for first time
 752       //if we use hidden info we want to have match between sim data and ESD
 753       if (tGoodMomentum==true)
 754         {
 755           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 756           realnofTracks++;//real number of tracks
 757           //      delete trackCopy;
 758         }
 759       else
 760         {
 761           delete  trackCopy;
 762         }
 763
 764     }
 765
 766   delete [] motherids;
 767
 768   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 769
 770   AliCentrality *cent = fEvent->GetCentrality();
 771   if (cent) {
 772     hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 773     //    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 774     hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 775     //    hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 776
 777     if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
 778   }
 779
 780   if (fEstEventMult == kGlobalCount)
 781     hbtEvent->SetNormalizedMult(tNormMult);
 782   else if (fEstEventMult == kTracklet)
 783     hbtEvent->SetNormalizedMult(tTracklet);
 784   else if (fEstEventMult == kITSTPC)
 785     hbtEvent->SetNormalizedMult(tITSTPC);
 786   else if (fEstEventMult == kITSPure)
 787     hbtEvent->SetNormalizedMult(tITSPure);
 788   else if (fEstEventMult == kSPDLayer1)
 789     hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 790   else if (fEstEventMult == kV0Centrality) {
 791     // centrality between 0 (central) and 1 (very peripheral)
 792
 793     if (cent) {
 794       if (cent->GetCentralityPercentile("V0M") < 0.00001)
 795         hbtEvent->SetNormalizedMult(-1);
 796       else
 797         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
 798       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 799                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
 800     }
 801   }
 802
 803   if (tNormMultPos > tNormMultNeg)
 804     hbtEvent->SetZDCParticipants(tNormMultPos);
 805   else
 806     hbtEvent->SetZDCParticipants(tNormMultNeg);
 807
 808   AliEventplane* ep = fEvent->GetEventplane();
 809   if (ep) {
 810     hbtEvent->SetEP(ep);
 811     hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
 812   }
 813
 814
 815   fCurEvent++;
 816   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 817   return hbtEvent;
 818 }
 819 //___________________
 820 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
 821 {
 822   // The chain loads the ESD for us
 823   // You must provide the address where it can be found
 824   fEvent = aESD;
 825 }
 826 //___________________
 827 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
 828 {
 829   // The chain loads the stack for us
 830   // You must provide the address where it can be found
 831   fStack = aStack;
 832 }
 833 //___________________
 834 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 835 {
 836   // The chain loads the generator event header for us
 837   // You must provide the address where it can be found
 838   fGenHeader = aGenHeader;
 839 }
 840 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
 841 {
 842   fRotateToEventPlane=dorotate;
 843 }
 844 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
 845 {
 846   // Calculates the number of sigma to the vertex.
 847
 848   Float_t b[2];
 849   Float_t bRes[2];
 850   Float_t bCov[3];
 851
 852   b[0] = impact[0];
 853   b[1] = impact[1];
 854   bCov[0] = covar[0];
 855   bCov[1] = covar[1];
 856   bCov[2] = covar[2];
 857
 858   bRes[0] = TMath::Sqrt(bCov[0]);
 859   bRes[1] = TMath::Sqrt(bCov[2]);
 860
 861   // -----------------------------------
 862   // How to get to a n-sigma cut?
 863   //
 864   // The accumulated statistics from 0 to d is
 865   //
 866   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 867   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 868   //
 869   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 870   // Can this be expressed in a different way?
 871
 872   if (bRes[0] == 0 || bRes[1] ==0)
 873     return -1;
 874
 875   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 876
 877   // stupid rounding problem screws up everything:
 878   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 879   if (TMath::Exp(-d * d / 2) < 1e-10)
 880     return 1000;
 881
 882   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 883   return d;
 884 }
 885
 886 void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)
 887 {
 888   fTrackType = aType;
 889 }
 890