PWGCF/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.2.2.2  2007/10/04 13:10:52  akisiel
  15  *Add Kink index storageAliFemtoEventReaderESD.cxx AliFemtoTrack.cxx AliFemtoTrack.h
  16  *
  17  *Revision 1.2.2.1  2007/09/30 11:38:59  akisiel
  18  *Adapt the readers to the new AliESDEvent structure
  19  *
  20  *Revision 1.2  2007/05/22 09:01:42  akisiel
  21  *Add the possibiloity to save cut settings in the ROOT file
  22  *
  23  *Revision 1.1  2007/05/16 10:22:11  akisiel
  24  *Making the directory structure of AliFemto flat. All files go into one common directory
  25  *
  26  *Revision 1.5  2007/05/03 09:45:20  akisiel
  27  *Fixing Effective C++ warnings
  28  *
  29  *Revision 1.4  2007/04/27 07:28:34  akisiel
  30  *Remove event number reading due to interface changes
  31  *
  32  *Revision 1.3  2007/04/27 07:25:16  akisiel
  33  *Make revisions needed for compilation from the main AliRoot tree
  34  *
  35  *Revision 1.1.1.1  2007/04/25 15:38:41  panos
  36  *Importing the HBT code dir
  37  *
  38  */
  39
  40 #include "AliFemtoEventReaderESD.h"
  41
  42 #include "TFile.h"
  43 #include "TTree.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDtrack.h"
  46 #include "AliESDVertex.h"
  47
  48 //#include "TSystem.h"
  49
  50 #include "AliFmPhysicalHelixD.h"
  51 #include "AliFmThreeVectorF.h"
  52
  53 #include "SystemOfUnits.h"
  54
  55 #include "AliFemtoEvent.h"
  56 #include "AliFemtoModelHiddenInfo.h"
  57
  58 ClassImp(AliFemtoEventReaderESD)
  59
  60 #if !(ST_NO_NAMESPACES)
  61   using namespace units;
  62 #endif
  63
  64 using namespace std;
  65 //____________________________
  66 //constructor with 0 parameters , look at default settings
  67 AliFemtoEventReaderESD::AliFemtoEventReaderESD():
  68   fInputFile(" "),
  69   fFileName(" "),
  70   fConstrained(true),
  71   fReadInner(false),
  72   fNumberofEvent(0),
  73   fCurEvent(0),
  74   fTree(0x0),
  75   fEsdFile(0x0),
  76   fEvent(0x0)
  77 {
  78   // default constructor
  79 }
  80
  81 AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
  82   AliFemtoEventReader(aReader),
  83   fInputFile(" "),
  84   fFileName(" "),
  85   fConstrained(true),
  86   fReadInner(false),
  87   fNumberofEvent(0),
  88   fCurEvent(0),
  89   fTree(0x0),
  90   fEsdFile(0x0),
  91   fEvent(0x0)
  92 {
  93   // copy constructor
  94   fInputFile = aReader.fInputFile;
  95   fFileName  = aReader.fFileName;
  96   fConstrained = aReader.fConstrained;
  97   fReadInner = aReader.fReadInner;
  98   fNumberofEvent = aReader.fNumberofEvent;
  99   fCurEvent = aReader.fCurEvent;
 100   //  fTree = aReader.fTree->CloneTree();
 101   //  fEvent = new AliESD(*aReader.fEvent);
 102   fEvent = new AliESDEvent();
 103   fEsdFile = new TFile(aReader.fEsdFile->GetName());
 104 }
 105 //__________________
 106 //Destructor
 107 AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
 108 {
 109   // destructor
 110   //delete fListOfFiles;
 111   delete fTree;
 112   delete fEvent;
 113   delete fEsdFile;
 114 }
 115
 116 //__________________
 117 AliFemtoEventReaderESD& AliFemtoEventReaderESD::operator=(const AliFemtoEventReaderESD& aReader)
 118 {
 119   // assignment operator
 120   if (this == &aReader)
 121     return *this;
 122
 123   fInputFile = aReader.fInputFile;
 124   fFileName  = aReader.fFileName;
 125   fConstrained = aReader.fConstrained;
 126   fReadInner = aReader.fReadInner;
 127   fNumberofEvent = aReader.fNumberofEvent;
 128   fCurEvent = aReader.fCurEvent;
 129   if (fTree) delete fTree;
 130   //  fTree = aReader.fTree->CloneTree();
 131   if (fEvent) delete fEvent;
 132   fEvent = new AliESDEvent();
 133   if (fEsdFile) delete fEsdFile;
 134   fEsdFile = new TFile(aReader.fEsdFile->GetName());
 135
 136   return *this;
 137 }
 138 //__________________
 139 AliFemtoString AliFemtoEventReaderESD::Report()
 140 {
 141   // create reader report
 142   AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
 143   return temp;
 144 }
 145
 146 //__________________
 147 void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
 148 {
 149   //setting the name of file where names of ESD file are written
 150   //it takes only this files which have good trees
 151   char buffer[256];
 152   fInputFile=string(inputFile);
 153   cout<<"Input File set on "<<fInputFile<<endl;
 154   ifstream infile(inputFile);
 155
 156   fTree = new TChain("esdTree");
 157
 158   if(infile.good()==true)
 159     {
 160       //checking if all give files have good tree inside
 161       while (infile.eof()==false)
 162         {
 163           infile.getline(buffer,256);
 164           //ifstream test_file(buffer);
 165           TFile *esdFile=TFile::Open(buffer,"READ");
 166           if (esdFile!=0x0)
 167             {
 168               TTree* tree = (TTree*) esdFile->Get("esdTree");
 169               if (tree!=0x0)
 170                 {
 171                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
 172                   fTree->AddFile(buffer);
 173                   delete tree;
 174                 }
 175               esdFile->Close();
 176             }
 177           delete esdFile;
 178         }
 179     }
 180 }
 181
 182 void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
 183 {
 184   fConstrained=constrained;
 185 }
 186
 187 bool AliFemtoEventReaderESD::GetConstrained() const
 188 {
 189   return fConstrained;
 190 }
 191
 192 void AliFemtoEventReaderESD::SetReadTPCInner(const bool readinner)
 193 {
 194   fReadInner=readinner;
 195 }
 196
 197 bool AliFemtoEventReaderESD::GetReadTPCInner() const
 198 {
 199   return fReadInner;
 200 }
 201
 202 AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
 203 {
 204   // read in a next hbt event from the chain
 205   // convert it to AliFemtoEvent and return
 206   // for further analysis
 207   AliFemtoEvent *hbtEvent = 0;
 208
 209   if (fCurEvent==fNumberofEvent)//open next file
 210     {
 211       if(fNumberofEvent==0)
 212         {
 213           //      delete fEvent;//added 1.04.2007
 214           fEvent=new AliESDEvent();
 215           //      delete fTree;
 216           //fTree=0;
 217           //      delete fEsdFile;
 218
 219           //ESD data
 220           //      fEsdFile=TFile::Open(fFileName.c_str(),"READ");
 221           //      fTree = (TTree*) fEsdFile->Get("esdTree");
 222           //      fTree->SetBranchAddress("ESD", &fEvent);
 223           fTree->SetBranchStatus("MuonTracks*",0);
 224           fTree->SetBranchStatus("PmdTracks*",0);
 225           fTree->SetBranchStatus("TrdTracks*",0);
 226           fTree->SetBranchStatus("V0s*",0);
 227           fTree->SetBranchStatus("Cascades*",0);
 228           fTree->SetBranchStatus("Kinks*",0);
 229           fTree->SetBranchStatus("CaloClusters*",0);
 230           fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
 231           fTree->SetBranchStatus("ESDfriend*",0);
 232           fEvent->ReadFromTree(fTree);
 233
 234           fNumberofEvent=fTree->GetEntries();
 235           cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
 236           fCurEvent=0;
 237           //sim data
 238         }
 239       else //no more data to read
 240         {
 241           cout<<"no more files "<<hbtEvent<<endl;
 242           fReaderStatus=1;
 243           return hbtEvent;
 244         }
 245     }
 246   cout<<"starting to read event "<<fCurEvent<<endl;
 247   fTree->GetEvent(fCurEvent);//getting next event
 248   cout << "Read event " << fEvent << " from file " << fTree << endl;
 249   //  vector<int> tLabelTable;//to check labels
 250
 251   hbtEvent = new AliFemtoEvent;
 252   //setting basic things
 253   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 254   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 255   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 256   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 257   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 258   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 259   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 260   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 261   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 262   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 263   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 264   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 265
 266   //Vertex
 267   double fV1[3];
 268   fEvent->GetVertex()->GetXYZ(fV1);
 269
 270   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 271   hbtEvent->SetPrimVertPos(vertex);
 272
 273   //starting to reading tracks
 274   int nofTracks=0;  //number of reconstructed tracks in event
 275   nofTracks=fEvent->GetNumberOfTracks();
 276   int realnofTracks=0;//number of track which we use ina analysis
 277   cout << "Event has " << nofTracks << " tracks " << endl;
 278
 279   for (int i=0;i<nofTracks;i++)
 280     {
 281       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 282
 283       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 284       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 285       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 286
 287       trackCopy->SetCharge((short)esdtrack->GetSign());
 288
 289       //in aliroot we have AliPID
 290       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 291       //we use only 5 first
 292       double esdpid[5];
 293       esdtrack->GetESDpid(esdpid);
 294       trackCopy->SetPidProbElectron(esdpid[0]);
 295       trackCopy->SetPidProbMuon(esdpid[1]);
 296       trackCopy->SetPidProbPion(esdpid[2]);
 297       trackCopy->SetPidProbKaon(esdpid[3]);
 298       trackCopy->SetPidProbProton(esdpid[4]);
 299
 300       double pxyz[3];
 301       if (fReadInner == true) {
 302
 303         if (esdtrack->GetTPCInnerParam()) {
 304           AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 305           param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 306           param->GetPxPyPz(pxyz);//reading noconstarined momentum
 307           delete param;
 308
 309           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 310           tInfo->SetPDGPid(211);
 311           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 312           tInfo->SetMass(0.13957);
 313           trackCopy->SetHiddenInfo(tInfo);
 314         }
 315       }
 316       if (fConstrained==true)
 317         tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 318       else
 319         tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 320       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 321       trackCopy->SetP(v);//setting momentum
 322       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 323       const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 324       if (ktP.Mag() == 0) {
 325         delete trackCopy;
 326         continue;
 327       }
 328       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
 329       //setting helix I do not if it is ok
 330       AliFmPhysicalHelixD helix(ktP,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 331       trackCopy->SetHelix(helix);
 332
 333       trackCopy->SetTrackId(esdtrack->GetID());
 334       trackCopy->SetFlags(esdtrack->GetStatus());
 335       trackCopy->SetLabel(esdtrack->GetLabel());
 336
 337       //some stuff which could be useful
 338       float impact[2];
 339       float covimpact[3];
 340       esdtrack->GetImpactParameters(impact,covimpact);
 341       trackCopy->SetImpactD(impact[0]);
 342       trackCopy->SetImpactZ(impact[1]);
 343       trackCopy->SetCdd(covimpact[0]);
 344       trackCopy->SetCdz(covimpact[1]);
 345       trackCopy->SetCzz(covimpact[2]);
 346       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 347       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 348       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 349       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 350       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 351       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 352       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 353
 354       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 355       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 356
 357       double pvrt[3];
 358       fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
 359
 360       double xtpc[3];
 361       esdtrack->GetInnerXYZ(xtpc);
 362       xtpc[2] -= pvrt[2];
 363       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 364
 365       esdtrack->GetOuterXYZ(xtpc);
 366       xtpc[2] -= pvrt[2];
 367       trackCopy->SetNominalTPCExitPoint(xtpc);
 368
 369       int indexes[3];
 370       for (int ik=0; ik<3; ik++) {
 371         indexes[ik] = esdtrack->GetKinkIndex(ik);
 372       }
 373       trackCopy->SetKinkIndexes(indexes);
 374       //decision if we want this track
 375       //if we using diffrent labels we want that this label was use for first time
 376       //if we use hidden info we want to have match between sim data and ESD
 377       if (tGoodMomentum==true)
 378         {
 379           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 380           realnofTracks++;//real number of tracks
 381         }
 382       else
 383         {
 384           delete  trackCopy;
 385         }
 386
 387     }
 388
 389   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 390   fCurEvent++;
 391   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 392 //   if (fCurEvent== fNumberofEvent)//if end of current file close all
 393 //     {
 394 //       fTree->Reset();
 395 //       delete fTree;
 396 //       fEsdFile->Close();
 397 //     }
 398   return hbtEvent;
 399 }
 400
 401
 402
 403
 404
 405
 406
 407
 408