PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChain.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoEventReaderESDChain - the reader class for the Alice ESD         ///
   4 /// tailored for the Task framework                                 ///
   5 /// Reads in AliESDfriend to create shared hit/quality information           ///
   6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoEventReaderESDChain.h"
  10
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDtrack.h"
  15 #include "AliESDVertex.h"
  16
  17 #include "AliFmPhysicalHelixD.h"
  18 #include "AliFmThreeVectorF.h"
  19
  20 #include "SystemOfUnits.h"
  21
  22 #include "AliFemtoEvent.h"
  23 #include "AliFemtoModelHiddenInfo.h"
  24
  25 ClassImp(AliFemtoEventReaderESDChain)
  26
  27 #if !(ST_NO_NAMESPACES)
  28   using namespace units;
  29 #endif
  30
  31 using namespace std;
  32 //____________________________
  33 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  34   fFileName(" "),
  35   fConstrained(true),
  36   fReadInner(false),
  37   fUseTPCOnly(false),
  38   fNumberofEvent(0),
  39   fCurEvent(0),
  40   fCurFile(0),
  41   fEvent(0x0),
  42   fUsePhysicsSel(kFALSE),
  43   fSelect(0x0),
  44   fTrackType(kGlobal)
  45 {
  46   //constructor with 0 parameters , look at default settings
  47 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  48 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  49 //     fClusterPerPadrow[tPad] = new list<Int_t>();
  50 //   }
  51 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  52 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  53 //     fSharedList[tPad] = new list<Int_t>();
  54 //   }
  55 }
  56
  57 //__________________
  58 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  59   AliFemtoEventReader(aReader),
  60   fFileName(" "),
  61   fConstrained(true),
  62   fReadInner(false),
  63   fUseTPCOnly(false),
  64   fNumberofEvent(0),
  65   fCurEvent(0),
  66   fCurFile(0),
  67   fEvent(0x0),
  68   fUsePhysicsSel(kFALSE),
  69   fSelect(0x0),
  70   fTrackType(kGlobal)
  71 {
  72   // Copy constructor
  73   fConstrained = aReader.fConstrained;
  74   fReadInner = aReader.fReadInner;
  75   fUseTPCOnly = aReader.fUseTPCOnly;
  76   fNumberofEvent = aReader.fNumberofEvent;
  77   fCurEvent = aReader.fCurEvent;
  78   fCurFile = aReader.fCurFile;
  79   //  fEvent = new AliESD(*aReader.fEvent);
  80   fEvent = new AliESDEvent();
  81   fUsePhysicsSel = aReader.fUsePhysicsSel;
  82   if (aReader.fUsePhysicsSel)
  83     fSelect = new AliPhysicsSelection();
  84   fTrackType = aReader.fTrackType;
  85 //   fEventFriend = aReader.fEventFriend;
  86 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  87 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  88 //     fClusterPerPadrow[tPad] = new list<Int_t>();
  89 //     list<Int_t>::iterator iter;
  90 //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
  91 //       fClusterPerPadrow[tPad]->push_back(*iter);
  92 //     }
  93 //   }
  94 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  95 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  96 //     fSharedList[tPad] = new list<Int_t>();
  97 //     list<Int_t>::iterator iter;
  98 //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
  99 //       fSharedList[tPad]->push_back(*iter);
 100 //     }
 101 //   }
 102 }
 103 //__________________
 104 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
 105 {
 106   //Destructor
 107   delete fEvent;
 108
 109 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 110 //     fClusterPerPadrow[tPad]->clear();
 111 //     delete fClusterPerPadrow[tPad];
 112 //   }
 113 //   delete [] fClusterPerPadrow;
 114 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 115 //     fSharedList[tPad]->clear();
 116 //     delete fSharedList[tPad];
 117 //   }
 118 //   delete [] fSharedList;
 119   if (fSelect) delete fSelect;
 120 }
 121
 122 //__________________
 123 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
 124 {
 125   // Assignment operator
 126   if (this == &aReader)
 127     return *this;
 128
 129   fConstrained = aReader.fConstrained;
 130   fReadInner = aReader.fReadInner;
 131   fUseTPCOnly = aReader.fUseTPCOnly;
 132   fNumberofEvent = aReader.fNumberofEvent;
 133   fCurEvent = aReader.fCurEvent;
 134   fCurFile = aReader.fCurFile;
 135   if (fEvent) delete fEvent;
 136   fEvent = new AliESDEvent();
 137   fTrackType = aReader.fTrackType;
 138
 139   fUsePhysicsSel = aReader.fUsePhysicsSel;
 140   if (aReader.fUsePhysicsSel)
 141     fSelect = new AliPhysicsSelection();
 142   //  fEventFriend = aReader.fEventFriend;
 143
 144 //   if (fClusterPerPadrow) {
 145 //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 146 //       fClusterPerPadrow[tPad]->clear();
 147 //       delete fClusterPerPadrow[tPad];
 148 //     }
 149 //     delete [] fClusterPerPadrow;
 150 //   }
 151
 152 //   if (fSharedList) {
 153 //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 154 //       fSharedList[tPad]->clear();
 155 //       delete fSharedList[tPad];
 156 //     }
 157 //     delete [] fSharedList;
 158 //   }
 159
 160 //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 161 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 162 //     fClusterPerPadrow[tPad] = new list<Int_t>();
 163 //     list<Int_t>::iterator iter;
 164 //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 165 //       fClusterPerPadrow[tPad]->push_back(*iter);
 166 //     }
 167 //   }
 168 //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 169 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 170 //     fSharedList[tPad] = new list<Int_t>();
 171 //     list<Int_t>::iterator iter;
 172 //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 173 //       fSharedList[tPad]->push_back(*iter);
 174 //     }
 175 //   }
 176
 177   return *this;
 178 }
 179 //__________________
 180 // Simple report
 181 AliFemtoString AliFemtoEventReaderESDChain::Report()
 182 {
 183   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
 184   return temp;
 185 }
 186
 187 //__________________
 188 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
 189 {
 190   // Select whether to read constrained or not constrained momentum
 191   fConstrained=constrained;
 192 }
 193 //__________________
 194 bool AliFemtoEventReaderESDChain::GetConstrained() const
 195 {
 196   // Check whether we read constrained or not constrained momentum
 197   return fConstrained;
 198 }
 199 //__________________
 200 void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
 201 {
 202   fReadInner=readinner;
 203 }
 204
 205 bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
 206 {
 207   return fReadInner;
 208 }
 209
 210 //__________________
 211 void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
 212 {
 213   fUseTPCOnly=usetpconly;
 214 }
 215
 216 bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
 217 {
 218   return fUseTPCOnly;
 219 }
 220
 221 void AliFemtoEventReaderESDChain::SetUsePhysicsSelection(const bool usephysics)
 222 {
 223   fUsePhysicsSel = usephysics;
 224   if (!fSelect) fSelect = new AliPhysicsSelection();
 225 }
 226
 227 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
 228 {
 229   // Get the event, read all the relevant information
 230   // and fill the AliFemtoEvent class
 231   // Returns a valid AliFemtoEvent
 232   AliFemtoEvent *hbtEvent = 0;
 233   string tFriendFileName;
 234
 235   // Get the friend information
 236   cout<<"starting to read event "<<fCurEvent<<endl;
 237   //  fEvent->SetESDfriend(fEventFriend);
 238   if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
 239
 240   hbtEvent = new AliFemtoEvent;
 241
 242   if (fUsePhysicsSel) {
 243     hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fEvent));
 244     if (!(fSelect->IsCollisionCandidate(fEvent)))
 245       printf("Event not a collision candidate\n");
 246   }
 247   else
 248     hbtEvent->SetIsCollisionCandidate(kTRUE);
 249
 250   //setting basic things
 251   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 252   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 253   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 254   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 255   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 256   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 257   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 258   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 259   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 260   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 261   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 262   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 263
 264   //Vertex
 265   double fV1[3];
 266   double fVCov[6];
 267   if (fUseTPCOnly) {
 268     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 269     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 270     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 271       fVCov[4] = -1001.0;
 272   }
 273   else {
 274     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 275     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 276     if (!fEvent->GetPrimaryVertex()->GetStatus())
 277       fVCov[4] = -1001.0;
 278   }
 279
 280   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 281   hbtEvent->SetPrimVertPos(vertex);
 282   hbtEvent->SetPrimVertCov(fVCov);
 283
 284   //starting to reading tracks
 285   int nofTracks=0;  //number of reconstructed tracks in event
 286   nofTracks=fEvent->GetNumberOfTracks();
 287   int realnofTracks=0;//number of track which we use ina analysis
 288
 289 //   // Clear the shared cluster list
 290 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 291 //     fClusterPerPadrow[tPad]->clear();
 292 //   }
 293 //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 294 //     fSharedList[tPad]->clear();
 295 //   }
 296
 297
 298 //   for (int i=0;i<nofTracks;i++) {
 299 //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 300
 301 //     list<Int_t>::iterator tClustIter;
 302
 303 //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 304 //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 305 //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 306 //       if (tTrackIndices[tNcl] >= 0) {
 307 //      tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 308 //      if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 309 //        fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 310 //      }
 311 //      else {
 312 //        fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 313 //      }
 314 //       }
 315 //     }
 316
 317 //   }
 318
 319   int tNormMult = 0;
 320   for (int i=0;i<nofTracks;i++)
 321     {
 322       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 323
 324       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 325       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 326       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
 327         if (esdtrack->GetTPCNcls() > 80)
 328           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
 329             if (esdtrack->GetConstrainedParam())
 330               if (esdtrack->GetConstrainedParam()->Eta() < 0.9)
 331                 tNormMult++;
 332
 333       // If reading ITS-only tracks, reject all with TPC
 334       if (fTrackType == kITSOnly) {
 335         if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
 336         if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
 337         if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
 338         UChar_t iclm = esdtrack->GetITSClusterMap();
 339         Int_t incls = 0;
 340         for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
 341         if (incls<=3) {
 342           cout << "Rejecting track with " << incls << " clusters" << endl;
 343           continue;
 344         }
 345       }
 346
 347       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 348       trackCopy->SetCharge((short)esdtrack->GetSign());
 349
 350       //in aliroot we have AliPID
 351       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 352       //we use only 5 first
 353       double esdpid[5];
 354       esdtrack->GetESDpid(esdpid);
 355       trackCopy->SetPidProbElectron(esdpid[0]);
 356       trackCopy->SetPidProbMuon(esdpid[1]);
 357       trackCopy->SetPidProbPion(esdpid[2]);
 358       trackCopy->SetPidProbKaon(esdpid[3]);
 359       trackCopy->SetPidProbProton(esdpid[4]);
 360
 361       double pxyz[3];
 362       double rxyz[3];
 363       double impact[2];
 364       double covimpact[3];
 365
 366       if (fUseTPCOnly) {
 367         if (!esdtrack->GetTPCInnerParam()) {
 368           delete trackCopy;
 369           continue;
 370         }
 371
 372
 373         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 374         param->GetXYZ(rxyz);
 375         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 376         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 377
 378         if (fReadInner == true) {
 379           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 380           tInfo->SetPDGPid(211);
 381           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 382           tInfo->SetMass(0.13957);
 383           //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 384           //      tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 385           tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 386           trackCopy->SetHiddenInfo(tInfo);
 387         }
 388
 389         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 390         if (v.Mag() < 0.0001) {
 391           //    cout << "Found 0 momentum ???? " <<endl;
 392           delete trackCopy;
 393           continue;
 394         }
 395         trackCopy->SetP(v);//setting momentum
 396         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 397
 398         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 399         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 400         //setting helix I do not if it is ok
 401         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 402         trackCopy->SetHelix(helix);
 403
 404         //some stuff which could be useful
 405         trackCopy->SetImpactD(impact[0]);
 406         trackCopy->SetImpactZ(impact[1]);
 407         trackCopy->SetCdd(covimpact[0]);
 408         trackCopy->SetCdz(covimpact[1]);
 409         trackCopy->SetCzz(covimpact[2]);
 410         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 411
 412         delete param;
 413       }
 414       else {
 415         if (fReadInner == true) {
 416
 417           if (esdtrack->GetTPCInnerParam()) {
 418             AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 419             param->GetXYZ(rxyz);
 420             param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 421             param->GetPxPyPz(pxyz);//reading noconstarined momentum
 422             delete param;
 423
 424             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 425             tInfo->SetPDGPid(211);
 426             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 427             tInfo->SetMass(0.13957);
 428             //      tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 429             //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 430             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 431             trackCopy->SetHiddenInfo(tInfo);
 432           }
 433         }
 434         if (fConstrained==true)
 435           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 436         else
 437           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 438
 439         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 440         if (v.Mag() < 0.0001) {
 441           //    cout << "Found 0 momentum ???? " <<endl;
 442           delete trackCopy;
 443           continue;
 444         }
 445         trackCopy->SetP(v);//setting momentum
 446         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 447         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 448         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 449         //setting helix I do not if it is ok
 450         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 451         trackCopy->SetHelix(helix);
 452
 453         //some stuff which could be useful
 454         float imp[2];
 455         float cim[3];
 456         esdtrack->GetImpactParameters(imp,cim);
 457
 458         impact[0] = imp[0];
 459         impact[1] = imp[1];
 460         covimpact[0] = cim[0];
 461         covimpact[1] = cim[1];
 462         covimpact[2] = cim[2];
 463
 464         trackCopy->SetImpactD(impact[0]);
 465         trackCopy->SetImpactZ(impact[1]);
 466         trackCopy->SetCdd(covimpact[0]);
 467         trackCopy->SetCdz(covimpact[1]);
 468         trackCopy->SetCzz(covimpact[2]);
 469         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 470       }
 471
 472       trackCopy->SetTrackId(esdtrack->GetID());
 473       trackCopy->SetFlags(esdtrack->GetStatus());
 474       trackCopy->SetLabel(esdtrack->GetLabel());
 475
 476       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 477       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 478       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 479       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 480       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 481       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 482       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 483
 484       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 485       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 486
 487       double xtpc[3];
 488       esdtrack->GetInnerXYZ(xtpc);
 489       xtpc[2] -= fV1[2];
 490       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 491
 492       esdtrack->GetOuterXYZ(xtpc);
 493       xtpc[2] -= fV1[2];
 494       trackCopy->SetNominalTPCExitPoint(xtpc);
 495
 496       int indexes[3];
 497       for (int ik=0; ik<3; ik++) {
 498         indexes[ik] = esdtrack->GetKinkIndex(ik);
 499       }
 500       trackCopy->SetKinkIndexes(indexes);
 501       //decision if we want this track
 502       //if we using diffrent labels we want that this label was use for first time
 503       //if we use hidden info we want to have match between sim data and ESD
 504       if (tGoodMomentum==true)
 505         {
 506           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 507           realnofTracks++;//real number of tracks
 508           //      delete trackCopy;
 509         }
 510       else
 511         {
 512           delete  trackCopy;
 513         }
 514
 515     }
 516
 517   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 518   hbtEvent->SetNormalizedMult(tNormMult);
 519   fCurEvent++;
 520   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 521   return hbtEvent;
 522 }
 523 //___________________
 524 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
 525 {
 526   // The chain loads the ESD for us
 527   // You must provide the address where it can be found
 528   fEvent = aESD;
 529 }
 530 //___________________
 531 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 532 // {
 533 //   // We need the ESD tree to obtain
 534 //   // information about the friend file location
 535 //   fEventFriend = aFriend;
 536 // }
 537
 538 //____________________________________________________________________
 539 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
 540 {
 541   // Calculates the number of sigma to the vertex.
 542
 543   Float_t b[2];
 544   Float_t bRes[2];
 545   Float_t bCov[3];
 546
 547   b[0] = impact[0];
 548   b[1] = impact[1];
 549   bCov[0] = covar[0];
 550   bCov[1] = covar[1];
 551   bCov[2] = covar[2];
 552
 553   bRes[0] = TMath::Sqrt(bCov[0]);
 554   bRes[1] = TMath::Sqrt(bCov[2]);
 555
 556   // -----------------------------------
 557   // How to get to a n-sigma cut?
 558   //
 559   // The accumulated statistics from 0 to d is
 560   //
 561   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 562   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 563   //
 564   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 565   // Can this be expressed in a different way?
 566
 567   if (bRes[0] == 0 || bRes[1] ==0)
 568     return -1;
 569
 570   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 571
 572   // stupid rounding problem screws up everything:
 573   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 574   if (TMath::Exp(-d * d / 2) < 1e-10)
 575     return 1000;
 576
 577   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 578   return d;
 579 }
 580
 581 void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
 582 {
 583   fTrackType = aType;
 584 }
 585
 586
 587
 588
 589
 590
 591