PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESDChain.cxx

   1 #include "AliFemtoEventReaderESDChain.h"
   2
   3 #include "TFile.h"
   4 #include "TTree.h"
   5 #include "AliESD.h"
   6 #include "AliESDtrack.h"
   7
   8 #include "Infrastructure/AliFmPhysicalHelixD.h"
   9 #include "Infrastructure/AliFmThreeVectorF.h"
  10
  11 #include "Base/SystemOfUnits.h"
  12
  13 #include "Infrastructure/AliFemtoEvent.h"
  14
  15 ClassImp(AliFemtoEventReaderESDChain)
  16
  17 #if !(ST_NO_NAMESPACES)
  18   using namespace units;
  19 #endif
  20
  21 using namespace std;
  22 //____________________________
  23 //constructor with 0 parameters , look at default settings
  24 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  25   fFileName(" "),
  26   fConstrained(true),
  27   fNumberofEvent(0),
  28   fCurEvent(0),
  29   fCurFile(0),
  30   fEvent(0x0),
  31   fEventFriend(0)
  32 {
  33   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  34   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  35     fClusterPerPadrow[tPad] = new list<Int_t>();
  36   }
  37   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  38   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  39     fSharedList[tPad] = new list<Int_t>();
  40   }
  41 }
  42
  43 //__________________
  44 //Destructor
  45 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
  46 {
  47   delete fEvent;
  48
  49   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  50     fClusterPerPadrow[tPad]->clear();
  51     delete fClusterPerPadrow[tPad];
  52   }
  53   delete [] fClusterPerPadrow;
  54   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  55     fSharedList[tPad]->clear();
  56     delete fSharedList[tPad];
  57   }
  58   delete [] fSharedList;
  59 }
  60
  61 //__________________
  62 // Simple report
  63 AliFemtoString AliFemtoEventReaderESDChain::Report()
  64 {
  65   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
  66   return temp;
  67 }
  68
  69 //__________________
  70 // Select whether to read constrained or not constrained momentum
  71 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
  72 {
  73   fConstrained=constrained;
  74 }
  75 //__________________
  76 // Check whether we read constrained or not constrained momentum
  77 bool AliFemtoEventReaderESDChain::GetConstrained() const
  78 {
  79   return fConstrained;
  80 }
  81 //__________________
  82 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
  83 {
  84   // Get the event, read all the relevant information
  85   // and fill the AliFemtoEvent class
  86   // Returns a valid AliFemtoEvent
  87   AliFemtoEvent *hbtEvent = 0;
  88   string tFriendFileName;
  89
  90   // Get the friend information
  91   cout<<"starting to read event "<<fCurEvent<<endl;
  92   fEvent->SetESDfriend(fEventFriend);
  93   vector<int> label_table;//to check labels
  94
  95   hbtEvent = new AliFemtoEvent;
  96   //setting basic things
  97   //hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  98   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
  99   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 100   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 101   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 102   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 103   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 104   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 105   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 106   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 107   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 108   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 109
 110   //Vertex
 111   double fV1[3];
 112   fEvent->GetVertex()->GetXYZ(fV1);
 113
 114   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 115   hbtEvent->SetPrimVertPos(vertex);
 116
 117   //starting to reading tracks
 118   int nofTracks=0;  //number of reconstructed tracks in event
 119   nofTracks=fEvent->GetNumberOfTracks();
 120   int realnofTracks=0;//number of track which we use ina analysis
 121
 122   // Clear the shared cluster list
 123   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 124     fClusterPerPadrow[tPad]->clear();
 125   }
 126   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 127     fSharedList[tPad]->clear();
 128   }
 129
 130
 131   for (int i=0;i<nofTracks;i++) {
 132     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 133
 134     list<Int_t>::iterator tClustIter;
 135
 136     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 137     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 138     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 139       if (tTrackIndices[tNcl] >= 0) {
 140         tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 141         if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 142           fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 143         }
 144         else {
 145           fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 146         }
 147       }
 148     }
 149
 150   }
 151
 152   for (int i=0;i<nofTracks;i++)
 153     {
 154       bool  good_momentum=true; //flaga to chcek if we can read momentum of this track
 155
 156       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 157       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 158       const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 159
 160       trackCopy->SetCharge((short)esdtrack->GetSign());
 161
 162       //in aliroot we have AliPID
 163       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 164       //we use only 5 first
 165       double esdpid[5];
 166       esdtrack->GetESDpid(esdpid);
 167       trackCopy->SetPidProbElectron(esdpid[0]);
 168       trackCopy->SetPidProbMuon(esdpid[1]);
 169       trackCopy->SetPidProbPion(esdpid[2]);
 170       trackCopy->SetPidProbKaon(esdpid[3]);
 171       trackCopy->SetPidProbProton(esdpid[4]);
 172
 173       double pxyz[3];
 174       if (fConstrained==true)
 175         good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 176       else
 177         good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 178
 179       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 180       if (v.mag() == 0) {
 181         //      cout << "Found 0 momentum ???? " <<endl;
 182         delete trackCopy;
 183         continue;
 184       }
 185       trackCopy->SetP(v);//setting momentum
 186       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 187       const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
 188       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
 189       //setting helix I do not if it is ok
 190       AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 191       trackCopy->SetHelix(helix);
 192
 193       trackCopy->SetTrackId(esdtrack->GetID());
 194       trackCopy->SetFlags(esdtrack->GetStatus());
 195       //trackCopy->SetLabel(esdtrack->GetLabel());
 196
 197       //some stuff which could be useful
 198       float impact[2];
 199       float covimpact[3];
 200       esdtrack->GetImpactParameters(impact,covimpact);
 201       trackCopy->SetImpactD(impact[0]);
 202       trackCopy->SetImpactZ(impact[1]);
 203       trackCopy->SetCdd(covimpact[0]);
 204       trackCopy->SetCdz(covimpact[1]);
 205       trackCopy->SetCzz(covimpact[2]);
 206       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 207       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 208       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 209       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 210       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 211       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 212       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 213
 214       // Fill cluster per padrow information
 215       Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 216       Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 217       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 218         if (tTrackIndices[tNcl] > 0)
 219           trackCopy->SetTPCcluster(tNcl, 1);
 220         else
 221           trackCopy->SetTPCcluster(tNcl, 0);
 222       }
 223
 224       // Fill shared cluster information
 225       list<Int_t>::iterator tClustIter;
 226
 227       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 228         if (tTrackIndices[tNcl] > 0) {
 229           tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
 230           if (tClustIter != fSharedList[tNcl]->end()) {
 231             trackCopy->SetTPCshared(tNcl, 1);
 232             cout << "Event next" <<  endl;
 233             cout << "Track: " << i << endl;
 234             cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
 235           }
 236           else {
 237             trackCopy->SetTPCshared(tNcl, 0);
 238           }
 239         }
 240       }
 241
 242       if (good_momentum==true)
 243         {
 244           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 245           realnofTracks++;//real number of tracks
 246           //      delete trackCopy;
 247         }
 248       else
 249         {
 250           delete  trackCopy;
 251         }
 252
 253     }
 254
 255   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 256   fCurEvent++;
 257   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 258   return hbtEvent;
 259 }
 260 //___________________
 261 // The chain loads the ESD for us
 262 // You must provide the address where it can be found
 263 void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
 264 {
 265   fEvent = aESD;
 266 }
 267 //___________________
 268 // We need the ESD tree to obtain
 269 // information about the friend file location
 270 void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 271 {
 272   fEventFriend = aFriend;
 273 }
 274
 275
 276
 277
 278
 279
 280
 281
 282