PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESD.cxx

   1 /*
   2  *$Id$
   3  *$Log$
   4  *Revision 1.1.1.1  2007/04/25 15:38:41  panos
   5  *Importing the HBT code dir
   6  *
   7  *Revision 1.5  2007-04-03 16:00:08  mchojnacki
   8  *Changes to iprove memory managing
   9  *
  10  *Revision 1.4  2007/03/13 15:30:03  mchojnacki
  11  *adding reader for simulated data
  12  *
  13  *Revision 1.3  2007/03/08 14:58:03  mchojnacki
  14  *adding some alice stuff
  15  *
  16  *Revision 1.2  2007/03/07 13:36:16  mchojnacki
  17  *Add some comments
  18  *
  19  *Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  20  *First version on CVS
  21  *
  22  */
  23
  24 #include "AliFemtoEventReaderESD.h"
  25
  26 #include "TFile.h"
  27 #include "TTree.h"
  28 #include "AliESD.h"
  29 #include "AliESDtrack.h"
  30
  31 //#include "TSystem.h"
  32
  33 #include "Infrastructure/AliFmPhysicalHelixD.h"
  34 #include "Infrastructure/AliFmThreeVectorF.h"
  35
  36 #include "Base/SystemOfUnits.h"
  37
  38 #include "Infrastructure/AliFemtoEvent.h"
  39
  40 ClassImp(AliFemtoEventReaderESD)
  41
  42 #if !(ST_NO_NAMESPACES)
  43   using namespace units;
  44 #endif
  45
  46 using namespace std;
  47 //____________________________
  48 //constructor with 0 parameters , look at default settings
  49 AliFemtoEventReaderESD::AliFemtoEventReaderESD():
  50   fInputFile(" "),
  51   fFileName(" "),
  52   fConstrained(true),
  53   fNumberofEvent(0),
  54   fCurEvent(0),
  55   fCurFile(0),
  56   fTree(0x0),
  57   fEvent(0x0),
  58   fEsdFile(0x0),
  59   fEventFriend(0)
  60 {
  61   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  62   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  63     fClusterPerPadrow[tPad] = new list<Int_t>();
  64   }
  65   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  66   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  67     fSharedList[tPad] = new list<Int_t>();
  68   }
  69 }
  70
  71 //__________________
  72 //Destructor
  73 AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
  74 {
  75   //delete fListOfFiles;
  76   delete fTree;
  77   delete fEvent;
  78   delete fEsdFile;
  79
  80   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  81     fClusterPerPadrow[tPad]->clear();
  82     delete fClusterPerPadrow[tPad];
  83   }
  84   delete [] fClusterPerPadrow;
  85   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  86     fSharedList[tPad]->clear();
  87     delete fSharedList[tPad];
  88   }
  89   delete [] fSharedList;
  90 }
  91
  92 //__________________
  93 AliFemtoString AliFemtoEventReaderESD::Report()
  94 {
  95   AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
  96   return temp;
  97 }
  98
  99 //__________________
 100 //setting the name of file where names of ESD file are written
 101 //it takes only this files which have good trees
 102 void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
 103 {
 104   char buffer[256];
 105   fInputFile=string(inputFile);
 106   cout<<"Input File set on "<<fInputFile<<endl;
 107   ifstream infile(inputFile);
 108   if(infile.good()==true)
 109     {
 110       //checking if all give files have good tree inside
 111       while (infile.eof()==false)
 112         {
 113           infile.getline(buffer,256);
 114           //ifstream test_file(buffer);
 115           TFile *esdFile=TFile::Open(buffer,"READ");
 116           if (esdFile!=0x0)
 117             {
 118               TTree* tree = (TTree*) esdFile->Get("esdTree");
 119               if (tree!=0x0)
 120                 {
 121                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
 122                   fListOfFiles.push_back(string(buffer));
 123                   delete tree;
 124                 }
 125               esdFile->Close();
 126             }
 127           delete esdFile;
 128         }
 129     }
 130 }
 131
 132 //setting the next file to read
 133 bool AliFemtoEventReaderESD::GetNextFile()
 134 {
 135   if (fCurFile>=fListOfFiles.size())
 136     return false;
 137   fFileName=fListOfFiles.at(fCurFile);
 138   cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;
 139
 140   fCurFile++;
 141   return true;
 142 }
 143 void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
 144 {
 145   fConstrained=constrained;
 146 }
 147
 148 bool AliFemtoEventReaderESD::GetConstrained() const
 149 {
 150   return fConstrained;
 151 }
 152
 153 AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
 154 {
 155   AliFemtoEvent *hbtEvent = 0;
 156   string tFriendFileName;
 157
 158   if (fCurEvent==fNumberofEvent)//open next file
 159     {
 160       cout<<"next file"<<endl;
 161       if(GetNextFile())
 162         {
 163           delete fEventFriend;
 164           fEventFriend = 0;
 165           delete fEvent;//added 1.04.2007
 166           fEvent=new AliESD();
 167           delete fTree;
 168           //fTree=0;
 169           delete fEsdFile;
 170
 171           //ESD data
 172           fEsdFile=TFile::Open(fFileName.c_str(),"READ");
 173           fTree = (TTree*) fEsdFile->Get("esdTree");
 174           fTree->SetBranchAddress("ESD", &fEvent);
 175
 176           // Attach the friend tree with additional information
 177           tFriendFileName = fFileName;
 178           tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
 179           cout << "Reading friend " << tFriendFileName.c_str() << endl;;
 180           fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
 181           fTree->SetBranchAddress("ESDfriend",&fEventFriend);
 182
 183           fNumberofEvent=fTree->GetEntries();
 184           cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
 185           fCurEvent=0;
 186           //sim data
 187         }
 188       else //no more data to read
 189         {
 190           cout<<"no more files "<<hbtEvent<<endl;
 191           fReaderStatus=1;
 192           return hbtEvent;
 193         }
 194     }
 195   cout<<"starting to read event "<<fCurEvent<<endl;
 196   fTree->GetEvent(fCurEvent);//getting next event
 197   fEvent->SetESDfriend(fEventFriend);
 198   vector<int> label_table;//to check labels
 199
 200   hbtEvent = new AliFemtoEvent;
 201   //setting basic things
 202   //hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 203   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 204   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 205   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 206   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 207   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 208   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 209   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 210   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 211   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 212   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 213   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 214
 215   //Vertex
 216   double fV1[3];
 217   fEvent->GetVertex()->GetXYZ(fV1);
 218
 219   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 220   hbtEvent->SetPrimVertPos(vertex);
 221
 222   //starting to reading tracks
 223   int nofTracks=0;  //number of reconstructed tracks in event
 224   nofTracks=fEvent->GetNumberOfTracks();
 225   int realnofTracks=0;//number of track which we use ina analysis
 226
 227   // Clear the shared cluster list
 228   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 229     fClusterPerPadrow[tPad]->clear();
 230   }
 231   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 232     fSharedList[tPad]->clear();
 233   }
 234
 235
 236   for (int i=0;i<nofTracks;i++) {
 237     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 238
 239     list<Int_t>::iterator tClustIter;
 240
 241     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 242     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 243     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 244       if (tTrackIndices[tNcl] >= 0) {
 245         tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 246         if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 247           fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 248         }
 249         else {
 250           fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 251         }
 252       }
 253     }
 254
 255   }
 256
 257   for (int i=0;i<nofTracks;i++)
 258     {
 259       bool  good_momentum=true; //flaga to chcek if we can read momentum of this track
 260
 261       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 262       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 263       const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 264
 265       trackCopy->SetCharge((short)esdtrack->GetSign());
 266
 267       //in aliroot we have AliPID
 268       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 269       //we use only 5 first
 270       double esdpid[5];
 271       esdtrack->GetESDpid(esdpid);
 272       trackCopy->SetPidProbElectron(esdpid[0]);
 273       trackCopy->SetPidProbMuon(esdpid[1]);
 274       trackCopy->SetPidProbPion(esdpid[2]);
 275       trackCopy->SetPidProbKaon(esdpid[3]);
 276       trackCopy->SetPidProbProton(esdpid[4]);
 277
 278       double pxyz[3];
 279       if (fConstrained==true)
 280         good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 281       else
 282         good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 283       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 284       trackCopy->SetP(v);//setting momentum
 285       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 286       const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
 287       if (p.mag() == 0) {
 288         delete trackCopy;
 289         continue;
 290       }
 291       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
 292       //setting helix I do not if it is ok
 293       AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 294       trackCopy->SetHelix(helix);
 295
 296       trackCopy->SetTrackId(esdtrack->GetID());
 297       trackCopy->SetFlags(esdtrack->GetStatus());
 298       //trackCopy->SetLabel(esdtrack->GetLabel());
 299
 300       //some stuff which could be useful
 301       float impact[2];
 302       float covimpact[3];
 303       esdtrack->GetImpactParameters(impact,covimpact);
 304       trackCopy->SetImpactD(impact[0]);
 305       trackCopy->SetImpactZ(impact[1]);
 306       trackCopy->SetCdd(covimpact[0]);
 307       trackCopy->SetCdz(covimpact[1]);
 308       trackCopy->SetCzz(covimpact[2]);
 309       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 310       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 311       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 312       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 313       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 314       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 315       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 316
 317       // Fill cluster per padrow information
 318       Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 319       Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 320       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 321         if (tTrackIndices[tNcl] > 0)
 322           trackCopy->SetTPCcluster(tNcl, 1);
 323         else
 324           trackCopy->SetTPCcluster(tNcl, 0);
 325       }
 326
 327       // Fill shared cluster information
 328       list<Int_t>::iterator tClustIter;
 329
 330       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 331         if (tTrackIndices[tNcl] > 0) {
 332           tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
 333           if (tClustIter != fSharedList[tNcl]->end()) {
 334             trackCopy->SetTPCshared(tNcl, 1);
 335             cout << "Event next" <<  endl;
 336             cout << "Track: " << i << endl;
 337             cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
 338           }
 339           else {
 340             trackCopy->SetTPCshared(tNcl, 0);
 341           }
 342         }
 343       }
 344
 345       //decision if we want this track
 346       //if we using diffrent labels we want that this label was use for first time
 347       //if we use hidden info we want to have match between sim data and ESD
 348       if (good_momentum==true)
 349         {
 350           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 351           realnofTracks++;//real number of tracks
 352         }
 353       else
 354         {
 355           delete  trackCopy;
 356         }
 357
 358     }
 359
 360   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 361   fCurEvent++;
 362   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 363   if (fCurEvent== fNumberofEvent)//if end of current file close all
 364     {
 365       fTree->Reset();
 366       delete fTree;
 367       fEsdFile->Close();
 368     }
 369   return hbtEvent;
 370 }
 371
 372
 373
 374
 375
 376
 377
 378
 379