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   fSharedList(0x0),
  33   fClusterPerPadrow(0x0)
  34 {
  35   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  36   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  37     fClusterPerPadrow[tPad] = new list<Int_t>();
  38   }
  39   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  40   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  41     fSharedList[tPad] = new list<Int_t>();
  42   }
  43 }
  44
  45 //__________________
  46 // Copy constructor
  47 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  48   fFileName(" "),
  49   fConstrained(true),
  50   fNumberofEvent(0),
  51   fCurEvent(0),
  52   fCurFile(0),
  53   fEvent(0x0),
  54   fEventFriend(0),
  55   fSharedList(0x0),
  56   fClusterPerPadrow(0x0)
  57 {
  58   fConstrained = aReader.fConstrained;
  59   fNumberofEvent = aReader.fNumberofEvent;
  60   fCurEvent = aReader.fCurEvent;
  61   fCurFile = aReader.fCurFile;
  62   fEvent = new AliESD(*aReader.fEvent);
  63   fEventFriend = aReader.fEventFriend;
  64   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  65   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  66     fClusterPerPadrow[tPad] = new list<Int_t>();
  67     list<Int_t>::iterator iter;
  68     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
  69       fClusterPerPadrow[tPad]->push_back(*iter);
  70     }
  71   }
  72   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  73   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  74     fSharedList[tPad] = new list<Int_t>();
  75     list<Int_t>::iterator iter;
  76     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
  77       fSharedList[tPad]->push_back(*iter);
  78     }
  79   }
  80 }
  81 //__________________
  82 //Destructor
  83 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
  84 {
  85   delete fEvent;
  86
  87   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  88     fClusterPerPadrow[tPad]->clear();
  89     delete fClusterPerPadrow[tPad];
  90   }
  91   delete [] fClusterPerPadrow;
  92   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  93     fSharedList[tPad]->clear();
  94     delete fSharedList[tPad];
  95   }
  96   delete [] fSharedList;
  97 }
  98
  99 //__________________
 100 // Assignment operator
 101 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
 102 {
 103   if (this == &aReader)
 104     return *this;
 105
 106   fConstrained = aReader.fConstrained;
 107   fNumberofEvent = aReader.fNumberofEvent;
 108   fCurEvent = aReader.fCurEvent;
 109   fCurFile = aReader.fCurFile;
 110   if (fEvent) delete fEvent;
 111   fEvent = new AliESD(*aReader.fEvent);
 112
 113   fEventFriend = aReader.fEventFriend;
 114
 115   if (fClusterPerPadrow) {
 116     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 117       fClusterPerPadrow[tPad]->clear();
 118       delete fClusterPerPadrow[tPad];
 119     }
 120     delete [] fClusterPerPadrow;
 121   }
 122
 123   if (fSharedList) {
 124     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 125       fSharedList[tPad]->clear();
 126       delete fSharedList[tPad];
 127     }
 128     delete [] fSharedList;
 129   }
 130
 131   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 132   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 133     fClusterPerPadrow[tPad] = new list<Int_t>();
 134     list<Int_t>::iterator iter;
 135     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 136       fClusterPerPadrow[tPad]->push_back(*iter);
 137     }
 138   }
 139   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 140   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 141     fSharedList[tPad] = new list<Int_t>();
 142     list<Int_t>::iterator iter;
 143     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 144       fSharedList[tPad]->push_back(*iter);
 145     }
 146   }
 147
 148   return *this;
 149 }
 150 //__________________
 151 // Simple report
 152 AliFemtoString AliFemtoEventReaderESDChain::Report()
 153 {
 154   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
 155   return temp;
 156 }
 157
 158 //__________________
 159 // Select whether to read constrained or not constrained momentum
 160 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
 161 {
 162   fConstrained=constrained;
 163 }
 164 //__________________
 165 // Check whether we read constrained or not constrained momentum
 166 bool AliFemtoEventReaderESDChain::GetConstrained() const
 167 {
 168   return fConstrained;
 169 }
 170 //__________________
 171 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
 172 {
 173   // Get the event, read all the relevant information
 174   // and fill the AliFemtoEvent class
 175   // Returns a valid AliFemtoEvent
 176   AliFemtoEvent *hbtEvent = 0;
 177   string tFriendFileName;
 178
 179   // Get the friend information
 180   cout<<"starting to read event "<<fCurEvent<<endl;
 181   fEvent->SetESDfriend(fEventFriend);
 182   vector<int> label_table;//to check labels
 183
 184   hbtEvent = new AliFemtoEvent;
 185   //setting basic things
 186   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 187   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 188   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 189   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 190   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 191   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 192   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 193   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 194   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 195   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 196   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 197   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 198
 199   //Vertex
 200   double fV1[3];
 201   fEvent->GetVertex()->GetXYZ(fV1);
 202
 203   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 204   hbtEvent->SetPrimVertPos(vertex);
 205
 206   //starting to reading tracks
 207   int nofTracks=0;  //number of reconstructed tracks in event
 208   nofTracks=fEvent->GetNumberOfTracks();
 209   int realnofTracks=0;//number of track which we use ina analysis
 210
 211   // Clear the shared cluster list
 212   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 213     fClusterPerPadrow[tPad]->clear();
 214   }
 215   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 216     fSharedList[tPad]->clear();
 217   }
 218
 219
 220   for (int i=0;i<nofTracks;i++) {
 221     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 222
 223     list<Int_t>::iterator tClustIter;
 224
 225     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 226     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 227     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 228       if (tTrackIndices[tNcl] >= 0) {
 229         tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 230         if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 231           fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 232         }
 233         else {
 234           fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 235         }
 236       }
 237     }
 238
 239   }
 240
 241   for (int i=0;i<nofTracks;i++)
 242     {
 243       bool  good_momentum=true; //flaga to chcek if we can read momentum of this track
 244
 245       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 246       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 247       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 248
 249       trackCopy->SetCharge((short)esdtrack->GetSign());
 250
 251       //in aliroot we have AliPID
 252       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 253       //we use only 5 first
 254       double esdpid[5];
 255       esdtrack->GetESDpid(esdpid);
 256       trackCopy->SetPidProbElectron(esdpid[0]);
 257       trackCopy->SetPidProbMuon(esdpid[1]);
 258       trackCopy->SetPidProbPion(esdpid[2]);
 259       trackCopy->SetPidProbKaon(esdpid[3]);
 260       trackCopy->SetPidProbProton(esdpid[4]);
 261
 262       double pxyz[3];
 263       if (fConstrained==true)
 264         good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 265       else
 266         good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 267
 268       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 269       if (v.mag() == 0) {
 270         //      cout << "Found 0 momentum ???? " <<endl;
 271         delete trackCopy;
 272         continue;
 273       }
 274       trackCopy->SetP(v);//setting momentum
 275       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 276       const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
 277       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
 278       //setting helix I do not if it is ok
 279       AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 280       trackCopy->SetHelix(helix);
 281
 282       trackCopy->SetTrackId(esdtrack->GetID());
 283       trackCopy->SetFlags(esdtrack->GetStatus());
 284       //trackCopy->SetLabel(esdtrack->GetLabel());
 285
 286       //some stuff which could be useful
 287       float impact[2];
 288       float covimpact[3];
 289       esdtrack->GetImpactParameters(impact,covimpact);
 290       trackCopy->SetImpactD(impact[0]);
 291       trackCopy->SetImpactZ(impact[1]);
 292       trackCopy->SetCdd(covimpact[0]);
 293       trackCopy->SetCdz(covimpact[1]);
 294       trackCopy->SetCzz(covimpact[2]);
 295       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 296       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 297       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 298       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 299       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 300       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 301       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 302
 303       // Fill cluster per padrow information
 304       Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 305       Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 306       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 307         if (tTrackIndices[tNcl] > 0)
 308           trackCopy->SetTPCcluster(tNcl, 1);
 309         else
 310           trackCopy->SetTPCcluster(tNcl, 0);
 311       }
 312
 313       // Fill shared cluster information
 314       list<Int_t>::iterator tClustIter;
 315
 316       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 317         if (tTrackIndices[tNcl] > 0) {
 318           tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
 319           if (tClustIter != fSharedList[tNcl]->end()) {
 320             trackCopy->SetTPCshared(tNcl, 1);
 321             cout << "Event next" <<  endl;
 322             cout << "Track: " << i << endl;
 323             cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
 324           }
 325           else {
 326             trackCopy->SetTPCshared(tNcl, 0);
 327           }
 328         }
 329       }
 330
 331       if (good_momentum==true)
 332         {
 333           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 334           realnofTracks++;//real number of tracks
 335           //      delete trackCopy;
 336         }
 337       else
 338         {
 339           delete  trackCopy;
 340         }
 341
 342     }
 343
 344   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 345   fCurEvent++;
 346   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 347   return hbtEvent;
 348 }
 349 //___________________
 350 // The chain loads the ESD for us
 351 // You must provide the address where it can be found
 352 void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
 353 {
 354   fEvent = aESD;
 355 }
 356 //___________________
 357 // We need the ESD tree to obtain
 358 // information about the friend file location
 359 void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 360 {
 361   fEventFriend = aFriend;
 362 }
 363
 364
 365
 366
 367
 368
 369
 370
 371