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 "AliAODParticle.h"
  25 #include "TParticle.h"
  26 #include "AliFemtoModelHiddenInfo.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   fNumberofEvent(0),
  40   fCurEvent(0),
  41   fCurFile(0),
  42   fEvent(0x0),
  43   fStack(0x0)
  44 {
  45   //constructor with 0 parameters , look at default settings
  46 }
  47
  48 //__________________
  49 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  50   fFileName(" "),
  51   fConstrained(true),
  52   fNumberofEvent(0),
  53   fCurEvent(0),
  54   fCurFile(0),
  55   fEvent(0x0),
  56   fStack(0x0)
  57 {
  58   // Copy constructor
  59   fConstrained = aReader.fConstrained;
  60   fNumberofEvent = aReader.fNumberofEvent;
  61   fCurEvent = aReader.fCurEvent;
  62   fCurFile = aReader.fCurFile;
  63   fEvent = new AliESDEvent();
  64   fStack = aReader.fStack;
  65 }
  66 //__________________
  67 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
  68 {
  69   //Destructor
  70   delete fEvent;
  71 }
  72
  73 //__________________
  74 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
  75 {
  76   // Assignment operator
  77   if (this == &aReader)
  78     return *this;
  79
  80   fConstrained = aReader.fConstrained;
  81   fNumberofEvent = aReader.fNumberofEvent;
  82   fCurEvent = aReader.fCurEvent;
  83   fCurFile = aReader.fCurFile;
  84   if (fEvent) delete fEvent;
  85   fEvent = new AliESDEvent();
  86   fStack = aReader.fStack;
  87
  88   return *this;
  89 }
  90 //__________________
  91 // Simple report
  92 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
  93 {
  94   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
  95   return temp;
  96 }
  97
  98 //__________________
  99 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
 100 {
 101   // Select whether to read constrained or not constrained momentum
 102   fConstrained=constrained;
 103 }
 104 //__________________
 105 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
 106 {
 107   // Check whether we read constrained or not constrained momentum
 108   return fConstrained;
 109 }
 110 //__________________
 111 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 112 {
 113   // Get the event, read all the relevant information
 114   // and fill the AliFemtoEvent class
 115   // Returns a valid AliFemtoEvent
 116   AliFemtoEvent *hbtEvent = 0;
 117   string tFriendFileName;
 118
 119   // Get the friend information
 120   cout<<"starting to read event "<<fCurEvent<<endl;
 121   //  fEvent->SetESDfriend(fEventFriend);
 122
 123   hbtEvent = new AliFemtoEvent;
 124   //setting basic things
 125   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 126   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 127   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 128   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 129   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 130   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 131   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 132   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 133   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 134   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 135   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 136   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 137
 138   //Vertex
 139   double fV1[3];
 140   fEvent->GetVertex()->GetXYZ(fV1);
 141
 142   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 143   hbtEvent->SetPrimVertPos(vertex);
 144
 145   //starting to reading tracks
 146   int nofTracks=0;  //number of reconstructed tracks in event
 147   nofTracks=fEvent->GetNumberOfTracks();
 148   int realnofTracks=0;//number of track which we use ina analysis
 149
 150   for (int i=0;i<nofTracks;i++)
 151     {
 152       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 153
 154       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 155       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 156       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 157
 158       trackCopy->SetCharge((short)esdtrack->GetSign());
 159
 160       //in aliroot we have AliPID
 161       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 162       //we use only 5 first
 163       double esdpid[5];
 164       esdtrack->GetESDpid(esdpid);
 165       trackCopy->SetPidProbElectron(esdpid[0]);
 166       trackCopy->SetPidProbMuon(esdpid[1]);
 167       trackCopy->SetPidProbPion(esdpid[2]);
 168       trackCopy->SetPidProbKaon(esdpid[3]);
 169       trackCopy->SetPidProbProton(esdpid[4]);
 170
 171       double pxyz[3];
 172       if (fConstrained==true)
 173         tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 174       else
 175         tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 176
 177       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 178       if (v.mag() < 0.0001) {
 179         //      cout << "Found 0 momentum ???? " <<endl;
 180         delete trackCopy;
 181         continue;
 182       }
 183       trackCopy->SetP(v);//setting momentum
 184       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 185       const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 186       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 187       //setting helix I do not if it is ok
 188       AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 189       trackCopy->SetHelix(helix);
 190
 191       trackCopy->SetTrackId(esdtrack->GetID());
 192       trackCopy->SetFlags(esdtrack->GetStatus());
 193       trackCopy->SetLabel(esdtrack->GetLabel());
 194
 195       //some stuff which could be useful
 196       float impact[2];
 197       float covimpact[3];
 198       esdtrack->GetImpactParameters(impact,covimpact);
 199       trackCopy->SetImpactD(impact[0]);
 200       trackCopy->SetImpactZ(impact[1]);
 201       trackCopy->SetCdd(covimpact[0]);
 202       trackCopy->SetCdz(covimpact[1]);
 203       trackCopy->SetCzz(covimpact[2]);
 204       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 205       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 206       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 207       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 208       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 209       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 210       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 211
 212
 213       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 214       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 215
 216       double pvrt[3];
 217       fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
 218
 219       double xtpc[3];
 220       esdtrack->GetInnerXYZ(xtpc);
 221       xtpc[2] -= pvrt[2];
 222       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 223
 224       esdtrack->GetOuterXYZ(xtpc);
 225       xtpc[2] -= pvrt[2];
 226       trackCopy->SetNominalTPCExitPoint(xtpc);
 227
 228       int indexes[3];
 229       for (int ik=0; ik<3; ik++) {
 230         indexes[ik] = esdtrack->GetKinkIndex(ik);
 231       }
 232       trackCopy->SetKinkIndexes(indexes);
 233
 234       // Fill the hidden information with the simulated data
 235       TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 236       AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
 237       AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 238       tInfo->SetPDGPid(tParticle->GetMostProbable());
 239       tInfo->SetTrueMomentum(tParticle->Px(), tParticle->Py(), tParticle->Pz());
 240       Double_t mass2 = (tParticle->E()*tParticle->E() -
 241                         tParticle->Px()*tParticle->Px() -
 242                         tParticle->Py()*tParticle->Py() -
 243                         tParticle->Pz()*tParticle->Pz());
 244
 245
 246       if (mass2>0.0)
 247         tInfo->SetMass(TMath::Sqrt(mass2));
 248       else
 249         tInfo->SetMass(0.0);
 250       tInfo->SetEmissionPoint(tParticle->Vx(), tParticle->Vy(), tParticle->Vz(), tParticle->T());
 251       trackCopy->SetHiddenInfo(tInfo);
 252
 253       //decision if we want this track
 254       //if we using diffrent labels we want that this label was use for first time
 255       //if we use hidden info we want to have match between sim data and ESD
 256       if (tGoodMomentum==true)
 257         {
 258           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 259           realnofTracks++;//real number of tracks
 260           //      delete trackCopy;
 261         }
 262       else
 263         {
 264           delete  trackCopy;
 265         }
 266
 267     }
 268
 269   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 270   fCurEvent++;
 271   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 272   return hbtEvent;
 273 }
 274 //___________________
 275 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
 276 {
 277   // The chain loads the ESD for us
 278   // You must provide the address where it can be found
 279   fEvent = aESD;
 280 }
 281 //___________________
 282 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
 283 {
 284   // The chain loads the stack for us
 285   // You must provide the address where it can be found
 286   fStack = aStack;
 287 }
 288
 289
 290
 291
 292
 293
 294
 295
 296