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