PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESD.cxx

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