PWG2/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 "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
  24 #include "TParticle.h"
  25 #include "AliFemtoModelHiddenInfo.h"
  26 #include "AliGenHijingEventHeader.h"
  27
  28 ClassImp(AliFemtoEventReaderESDChainKine)
  29
  30 #if !(ST_NO_NAMESPACES)
  31   using namespace units;
  32 #endif
  33
  34 using namespace std;
  35 //____________________________
  36 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
  37   fFileName(" "),
  38   fConstrained(true),
  39   fUseTPCOnly(false),
  40   fNumberofEvent(0),
  41   fCurEvent(0),
  42   fCurFile(0),
  43   fEvent(0x0),
  44   fStack(0x0),
  45   fGenHeader(0x0)
  46 {
  47   //constructor with 0 parameters , look at default settings
  48 }
  49
  50 //__________________
  51 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  52   AliFemtoEventReader(aReader),
  53   fFileName(" "),
  54   fConstrained(true),
  55   fUseTPCOnly(false),
  56   fNumberofEvent(0),
  57   fCurEvent(0),
  58   fCurFile(0),
  59   fEvent(0x0),
  60   fStack(0x0),
  61   fGenHeader(0x0)
  62 {
  63   // Copy constructor
  64   fConstrained = aReader.fConstrained;
  65   fUseTPCOnly = aReader.fUseTPCOnly;
  66   fNumberofEvent = aReader.fNumberofEvent;
  67   fCurEvent = aReader.fCurEvent;
  68   fCurFile = aReader.fCurFile;
  69   fEvent = new AliESDEvent();
  70   fStack = aReader.fStack;
  71 }
  72 //__________________
  73 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
  74 {
  75   //Destructor
  76   delete fEvent;
  77 }
  78
  79 //__________________
  80 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
  81 {
  82   // Assignment operator
  83   if (this == &aReader)
  84     return *this;
  85
  86   fConstrained = aReader.fConstrained;
  87   fUseTPCOnly = aReader.fUseTPCOnly;
  88   fNumberofEvent = aReader.fNumberofEvent;
  89   fCurEvent = aReader.fCurEvent;
  90   fCurFile = aReader.fCurFile;
  91   if (fEvent) delete fEvent;
  92   fEvent = new AliESDEvent();
  93   fStack = aReader.fStack;
  94   fGenHeader = aReader.fGenHeader;
  95
  96   return *this;
  97 }
  98 //__________________
  99 // Simple report
 100 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
 101 {
 102   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
 103   return temp;
 104 }
 105
 106 //__________________
 107 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
 108 {
 109   // Select whether to read constrained or not constrained momentum
 110   fConstrained=constrained;
 111 }
 112 //__________________
 113 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
 114 {
 115   // Check whether we read constrained or not constrained momentum
 116   return fConstrained;
 117 }
 118 //__________________
 119 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 120 {
 121   // Get the event, read all the relevant information
 122   // and fill the AliFemtoEvent class
 123   // Returns a valid AliFemtoEvent
 124   AliFemtoEvent *hbtEvent = 0;
 125   string tFriendFileName;
 126
 127   // Get the friend information
 128   cout<<"starting to read event "<<fCurEvent<<endl;
 129   //  fEvent->SetESDfriend(fEventFriend);
 130
 131   hbtEvent = new AliFemtoEvent;
 132   //setting basic things
 133   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 134   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 135   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 136   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 137   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 138   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 139   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 140   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 141   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 142   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 143   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 144   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 145
 146   //Vertex
 147   double fV1[3];
 148   if (fUseTPCOnly) {
 149     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 150   }
 151   else {
 152     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 153   }
 154
 155   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 156   hbtEvent->SetPrimVertPos(vertex);
 157
 158   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 159
 160   Double_t tReactionPlane = 0;
 161   if (hdh)
 162     {
 163       tReactionPlane = hdh->ReactionPlaneAngle();
 164     }
 165   //starting to reading tracks
 166   int nofTracks=0;  //number of reconstructed tracks in event
 167   nofTracks=fEvent->GetNumberOfTracks();
 168   int realnofTracks=0;//number of track which we use ina analysis
 169
 170   Int_t *motherids;
 171   motherids = new Int_t[fStack->GetNtrack()];
 172   for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 173
 174   // Read in mother ids
 175   TParticle *motherpart;
 176   for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 177     motherpart = fStack->Particle(ip);
 178     if (motherpart->GetDaughter(0) > 0)
 179       motherids[motherpart->GetDaughter(0)] = ip;
 180     if (motherpart->GetDaughter(1) > 0)
 181       motherids[motherpart->GetDaughter(1)] = ip;
 182
 183 //     if (motherpart->GetPdgCode() == 211) {
 184 //       cout << "Mother " << ip << " has daughters "
 185 //         << motherpart->GetDaughter(0) << " "
 186 //         << motherpart->GetDaughter(1) << " "
 187 //         << motherpart->Vx() << " "
 188 //         << motherpart->Vy() << " "
 189 //         << motherpart->Vz() << " "
 190 //         << endl;
 191
 192 //     }
 193   }
 194
 195   for (int i=0;i<nofTracks;i++)
 196     {
 197       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 198
 199       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 200       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 201       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 202
 203       trackCopy->SetCharge((short)esdtrack->GetSign());
 204
 205       //in aliroot we have AliPID
 206       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 207       //we use only 5 first
 208       double esdpid[5];
 209       esdtrack->GetESDpid(esdpid);
 210       trackCopy->SetPidProbElectron(esdpid[0]);
 211       trackCopy->SetPidProbMuon(esdpid[1]);
 212       trackCopy->SetPidProbPion(esdpid[2]);
 213       trackCopy->SetPidProbKaon(esdpid[3]);
 214       trackCopy->SetPidProbProton(esdpid[4]);
 215
 216       double pxyz[3];
 217       double rxyz[3];
 218       double impact[2];
 219       double covimpact[3];
 220
 221       if (fUseTPCOnly) {
 222         if (!esdtrack->GetTPCInnerParam()) {
 223           delete trackCopy;
 224           continue;
 225         }
 226
 227         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 228         param->GetXYZ(rxyz);
 229         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 230         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 231
 232         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 233         if (v.mag() < 0.0001) {
 234           //    cout << "Found 0 momentum ???? " <<endl;
 235           delete trackCopy;
 236           continue;
 237         }
 238         trackCopy->SetP(v);//setting momentum
 239         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 240
 241         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 242         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 243         //setting helix I do not if it is ok
 244         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 245         trackCopy->SetHelix(helix);
 246
 247         //some stuff which could be useful
 248         trackCopy->SetImpactD(impact[0]);
 249         trackCopy->SetImpactZ(impact[1]);
 250         trackCopy->SetCdd(covimpact[0]);
 251         trackCopy->SetCdz(covimpact[1]);
 252         trackCopy->SetCzz(covimpact[2]);
 253         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 254
 255         delete param;
 256       }
 257       else {
 258         if (fConstrained==true)
 259           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 260         else
 261           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 262
 263         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 264         if (v.mag() < 0.0001) {
 265           //    cout << "Found 0 momentum ???? " <<endl;
 266           delete trackCopy;
 267           continue;
 268         }
 269         trackCopy->SetP(v);//setting momentum
 270         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 271         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 272         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 273         //setting helix I do not if it is ok
 274         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 275         trackCopy->SetHelix(helix);
 276
 277         //some stuff which could be useful
 278         float imp[2];
 279         float cim[3];
 280         esdtrack->GetImpactParameters(imp,cim);
 281
 282         impact[0] = imp[0];
 283         impact[1] = imp[1];
 284         covimpact[0] = cim[0];
 285         covimpact[1] = cim[1];
 286         covimpact[2] = cim[2];
 287
 288         trackCopy->SetImpactD(impact[0]);
 289         trackCopy->SetImpactZ(impact[1]);
 290         trackCopy->SetCdd(covimpact[0]);
 291         trackCopy->SetCdz(covimpact[1]);
 292         trackCopy->SetCzz(covimpact[2]);
 293         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 294       }
 295
 296       trackCopy->SetTrackId(esdtrack->GetID());
 297       trackCopy->SetFlags(esdtrack->GetStatus());
 298       trackCopy->SetLabel(esdtrack->GetLabel());
 299
 300       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 301       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 302       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 303       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 304       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 305       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 306       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 307
 308
 309       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 310       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 311
 312       double xtpc[3];
 313       esdtrack->GetInnerXYZ(xtpc);
 314       xtpc[2] -= fV1[2];
 315       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 316
 317       esdtrack->GetOuterXYZ(xtpc);
 318       xtpc[2] -= fV1[2];
 319       trackCopy->SetNominalTPCExitPoint(xtpc);
 320
 321       int indexes[3];
 322       for (int ik=0; ik<3; ik++) {
 323         indexes[ik] = esdtrack->GetKinkIndex(ik);
 324       }
 325       trackCopy->SetKinkIndexes(indexes);
 326
 327       // Fill the hidden information with the simulated data
 328       TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 329
 330       // Check the mother information
 331
 332       // Using the new way of storing the freeze-out information
 333       // Final state particle is stored twice on the stack
 334       // one copy (mother) is stored with original freeze-out information
 335       //   and is not tracked
 336       // the other one (daughter) is stored with primary vertex position
 337       //   and is tracked
 338
 339       // Freeze-out coordinates
 340       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 341       fpx = tPart->Vx();
 342       fpy = tPart->Vy();
 343       fpz = tPart->Vz();
 344       fpt = tPart->T();
 345
 346       if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 347         TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 348         // Check if this is the same particle stored twice on the stack
 349         if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 350           // It is the same particle
 351           // Read in the original freeze-out information
 352           // and convert it from to [fm]
 353           fpx = mother->Vx()*1e13;
 354           fpy = mother->Vy()*1e13;
 355           fpz = mother->Vz()*1e13;
 356           fpt = mother->T()*1e13*3e10;
 357
 358         }
 359       }
 360
 361       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 362       tInfo->SetPDGPid(tPart->GetPdgCode());
 363       tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 364       Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 365                         tPart->Px()*tPart->Px() -
 366                         tPart->Py()*tPart->Py() -
 367                         tPart->Pz()*tPart->Pz());
 368       if (mass2>0.0)
 369         tInfo->SetMass(TMath::Sqrt(mass2));
 370       else
 371         tInfo->SetMass(0.0);
 372
 373       tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 374       trackCopy->SetHiddenInfo(tInfo);
 375
 376       //decision if we want this track
 377       //if we using diffrent labels we want that this label was use for first time
 378       //if we use hidden info we want to have match between sim data and ESD
 379       if (tGoodMomentum==true)
 380         {
 381           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 382           realnofTracks++;//real number of tracks
 383           //      delete trackCopy;
 384         }
 385       else
 386         {
 387           delete  trackCopy;
 388         }
 389
 390     }
 391
 392   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 393   fCurEvent++;
 394   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 395   return hbtEvent;
 396 }
 397 //___________________
 398 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
 399 {
 400   // The chain loads the ESD for us
 401   // You must provide the address where it can be found
 402   fEvent = aESD;
 403 }
 404 //___________________
 405 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
 406 {
 407   // The chain loads the stack for us
 408   // You must provide the address where it can be found
 409   fStack = aStack;
 410 }
 411 //___________________
 412 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 413 {
 414   // The chain loads the generator event header for us
 415   // You must provide the address where it can be found
 416   fGenHeader = aGenHeader;
 417 }
 418
 419 //__________________
 420 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
 421 {
 422   fUseTPCOnly=usetpconly;
 423 }
 424
 425 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
 426 {
 427   return fUseTPCOnly;
 428 }
 429 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
 430 {
 431   // Calculates the number of sigma to the vertex.
 432
 433   Float_t b[2];
 434   Float_t bRes[2];
 435   Float_t bCov[3];
 436
 437   b[0] = impact[0];
 438   b[1] = impact[1];
 439   bCov[0] = covar[0];
 440   bCov[1] = covar[1];
 441   bCov[2] = covar[2];
 442
 443   bRes[0] = TMath::Sqrt(bCov[0]);
 444   bRes[1] = TMath::Sqrt(bCov[2]);
 445
 446   // -----------------------------------
 447   // How to get to a n-sigma cut?
 448   //
 449   // The accumulated statistics from 0 to d is
 450   //
 451   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 452   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 453   //
 454   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 455   // Can this be expressed in a different way?
 456
 457   if (bRes[0] == 0 || bRes[1] ==0)
 458     return -1;
 459
 460   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 461
 462   // stupid rounding problem screws up everything:
 463   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 464   if (TMath::Exp(-d * d / 2) < 1e-10)
 465     return 1000;
 466
 467   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 468   return d;
 469 }