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