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