Make revisions needed for compilation from the main AliRoot tree
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Reader / AliFemtoEventReaderESD.cxx
CommitLineData
67427ff7 1/*
2 *$Id$
3 *$Log$
b2f60a91 4 *Revision 1.1.1.1 2007/04/25 15:38:41 panos
5 *Importing the HBT code dir
6 *
67427ff7 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
b2f60a91 33#include "Infrastructure/AliFmPhysicalHelixD.h"
34#include "Infrastructure/AliFmThreeVectorF.h"
67427ff7 35
b2f60a91 36#include "Base/SystemOfUnits.h"
67427ff7 37
b2f60a91 38#include "Infrastructure/AliFemtoEvent.h"
67427ff7 39
40ClassImp(AliFemtoEventReaderESD)
41
42#if !(ST_NO_NAMESPACES)
43 using namespace units;
44#endif
45
46using namespace std;
47//____________________________
48//constructor with 0 parameters , look at default settings
49AliFemtoEventReaderESD::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
73AliFemtoEventReaderESD::~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//__________________
93AliFemtoString 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
102void 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
133bool 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}
143void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
144{
145 fConstrained=constrained;
146}
147
148bool AliFemtoEventReaderESD::GetConstrained() const
149{
150 return fConstrained;
151}
152
153AliFemtoEvent* 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();
1e1d2996 167 // delete fTree;
67427ff7 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
1e1d2996 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
67427ff7 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
1e1d2996 210 hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7 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