a6062ec71ace9503ba3fde0b32e09ad489552f73
[u/mrichter/AliRoot.git] / 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           fNumberofEvent=fTree->GetEntries();
184           cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
185           fCurEvent=0;
186           //sim data
187         }
188       else //no more data to read
189         {
190           cout<<"no more files "<<hbtEvent<<endl;
191           fReaderStatus=1;
192           return hbtEvent; 
193         }
194     }           
195   cout<<"starting to read event "<<fCurEvent<<endl;
196   fTree->GetEvent(fCurEvent);//getting next event
197   fEvent->SetESDfriend(fEventFriend);
198   vector<int> label_table;//to check labels
199         
200   hbtEvent = new AliFemtoEvent;
201   //setting basic things
202   //hbtEvent->SetEventNumber(fEvent->GetEventNumber());
203   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
204   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
205   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
206   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
207   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
208   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
209   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
210   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
211   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
212   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
213   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
214         
215   //Vertex
216   double fV1[3];
217   fEvent->GetVertex()->GetXYZ(fV1);
218
219   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
220   hbtEvent->SetPrimVertPos(vertex);
221         
222   //starting to reading tracks
223   int nofTracks=0;  //number of reconstructed tracks in event
224   nofTracks=fEvent->GetNumberOfTracks();
225   int realnofTracks=0;//number of track which we use ina analysis
226
227   // Clear the shared cluster list
228   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
229     fClusterPerPadrow[tPad]->clear();
230   }
231   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
232     fSharedList[tPad]->clear();
233   }
234
235
236   for (int i=0;i<nofTracks;i++) {
237     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
238
239     list<Int_t>::iterator tClustIter;
240
241     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
242     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
243     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
244       if (tTrackIndices[tNcl] >= 0) {
245         tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
246         if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
247           fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
248         }
249         else {
250           fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
251         }
252       }
253     }
254       
255   }
256
257   for (int i=0;i<nofTracks;i++)
258     {
259       bool  good_momentum=true; //flaga to chcek if we can read momentum of this track
260                 
261       AliFemtoTrack* trackCopy = new AliFemtoTrack();   
262       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
263       const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
264
265       trackCopy->SetCharge((short)esdtrack->GetSign());
266
267       //in aliroot we have AliPID 
268       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
269       //we use only 5 first
270       double esdpid[5];
271       esdtrack->GetESDpid(esdpid);
272       trackCopy->SetPidProbElectron(esdpid[0]);
273       trackCopy->SetPidProbMuon(esdpid[1]);
274       trackCopy->SetPidProbPion(esdpid[2]);
275       trackCopy->SetPidProbKaon(esdpid[3]);
276       trackCopy->SetPidProbProton(esdpid[4]);
277                                                 
278       double pxyz[3];
279       if (fConstrained==true)               
280         good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
281       else
282         good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
283       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
284       trackCopy->SetP(v);//setting momentum
285       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
286       const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
287       if (p.mag() == 0) {
288         delete trackCopy;
289         continue;
290       }
291       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
292       //setting helix I do not if it is ok
293       AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
294       trackCopy->SetHelix(helix);
295                 
296       trackCopy->SetTrackId(esdtrack->GetID());
297       trackCopy->SetFlags(esdtrack->GetStatus());
298       //trackCopy->SetLabel(esdtrack->GetLabel());
299                 
300       //some stuff which could be useful 
301       float impact[2];
302       float covimpact[3];
303       esdtrack->GetImpactParameters(impact,covimpact);
304       trackCopy->SetImpactD(impact[0]);
305       trackCopy->SetImpactZ(impact[1]);
306       trackCopy->SetCdd(covimpact[0]);
307       trackCopy->SetCdz(covimpact[1]);
308       trackCopy->SetCzz(covimpact[2]);
309       trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
310       trackCopy->SetITSncls(esdtrack->GetNcls(0));     
311       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
312       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
313       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
314       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
315       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
316
317       // Fill cluster per padrow information
318       Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
319       Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
320       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
321         if (tTrackIndices[tNcl] > 0)
322           trackCopy->SetTPCcluster(tNcl, 1);
323         else
324           trackCopy->SetTPCcluster(tNcl, 0);
325       }
326       
327       // Fill shared cluster information
328       list<Int_t>::iterator tClustIter;
329
330       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
331         if (tTrackIndices[tNcl] > 0) {
332           tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
333           if (tClustIter != fSharedList[tNcl]->end()) {
334             trackCopy->SetTPCshared(tNcl, 1);
335             cout << "Event next" <<  endl;
336             cout << "Track: " << i << endl;
337             cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
338           }
339           else {
340             trackCopy->SetTPCshared(tNcl, 0);
341           }
342         }
343       }
344
345       //decision if we want this track
346       //if we using diffrent labels we want that this label was use for first time 
347       //if we use hidden info we want to have match between sim data and ESD
348       if (good_momentum==true)
349         {
350           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
351           realnofTracks++;//real number of tracks
352         }
353       else
354         {
355           delete  trackCopy;
356         }
357                 
358     }
359
360   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event   
361   fCurEvent++;  
362   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
363   if (fCurEvent== fNumberofEvent)//if end of current file close all
364     {   
365       fTree->Reset(); 
366       delete fTree;
367       fEsdFile->Close();
368     }
369   return hbtEvent; 
370 }
371
372
373
374
375
376
377
378
379