PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChain.cxx

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