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