39245b93e814c2205320a35bbe7cb8ed1856e83d
[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 //        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