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