PWG2/FEMTOSCOPY/AliFemto/Reader/AliFemtoEventReaderESD.cxx

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