PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESD.cxx

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