PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESD.cxx

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