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