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