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