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