]>
Commit | Line | Data |
---|---|---|
1 | /* | |
2 | *$Id$ | |
3 | *$Log$ | |
4 | *Revision 1.1.1.1 2007/04/25 15:38:41 panos | |
5 | *Importing the HBT code dir | |
6 | * | |
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 | ||
33 | #include "Infrastructure/AliFmPhysicalHelixD.h" | |
34 | #include "Infrastructure/AliFmThreeVectorF.h" | |
35 | ||
36 | #include "Base/SystemOfUnits.h" | |
37 | ||
38 | #include "Infrastructure/AliFemtoEvent.h" | |
39 | ||
40 | ClassImp(AliFemtoEventReaderESD) | |
41 | ||
42 | #if !(ST_NO_NAMESPACES) | |
43 | using namespace units; | |
44 | #endif | |
45 | ||
46 | using namespace std; | |
47 | //____________________________ | |
48 | //constructor with 0 parameters , look at default settings | |
49 | AliFemtoEventReaderESD::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 | |
73 | AliFemtoEventReaderESD::~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 | //__________________ | |
93 | AliFemtoString 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 | |
102 | void 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 | |
133 | bool 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 | } | |
143 | void AliFemtoEventReaderESD::SetConstrained(const bool constrained) | |
144 | { | |
145 | fConstrained=constrained; | |
146 | } | |
147 | ||
148 | bool AliFemtoEventReaderESD::GetConstrained() const | |
149 | { | |
150 | return fConstrained; | |
151 | } | |
152 | ||
153 | AliFemtoEvent* 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(); | |
167 | delete fTree; | |
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 | ||
183 | fNumberofEvent=fTree->GetEntries(); | |
184 | cout<<"Number of Entries in file "<<fNumberofEvent<<endl; | |
185 | fCurEvent=0; | |
186 | //sim data | |
187 | } | |
188 | else //no more data to read | |
189 | { | |
190 | cout<<"no more files "<<hbtEvent<<endl; | |
191 | fReaderStatus=1; | |
192 | return hbtEvent; | |
193 | } | |
194 | } | |
195 | cout<<"starting to read event "<<fCurEvent<<endl; | |
196 | fTree->GetEvent(fCurEvent);//getting next event | |
197 | fEvent->SetESDfriend(fEventFriend); | |
198 | vector<int> label_table;//to check labels | |
199 | ||
200 | hbtEvent = new AliFemtoEvent; | |
201 | //setting basic things | |
202 | //hbtEvent->SetEventNumber(fEvent->GetEventNumber()); | |
203 | hbtEvent->SetRunNumber(fEvent->GetRunNumber()); | |
204 | //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks()); | |
205 | hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok | |
206 | hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy()); | |
207 | hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy()); | |
208 | hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy()); | |
209 | hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy()); | |
210 | hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy()); | |
211 | hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants()); | |
212 | hbtEvent->SetTriggerMask(fEvent->GetTriggerMask()); | |
213 | hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster()); | |
214 | ||
215 | //Vertex | |
216 | double fV1[3]; | |
217 | fEvent->GetVertex()->GetXYZ(fV1); | |
218 | ||
219 | AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]); | |
220 | hbtEvent->SetPrimVertPos(vertex); | |
221 | ||
222 | //starting to reading tracks | |
223 | int nofTracks=0; //number of reconstructed tracks in event | |
224 | nofTracks=fEvent->GetNumberOfTracks(); | |
225 | int realnofTracks=0;//number of track which we use ina analysis | |
226 | ||
227 | // Clear the shared cluster list | |
228 | for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) { | |
229 | fClusterPerPadrow[tPad]->clear(); | |
230 | } | |
231 | for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) { | |
232 | fSharedList[tPad]->clear(); | |
233 | } | |
234 | ||
235 | ||
236 | for (int i=0;i<nofTracks;i++) { | |
237 | const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track | |
238 | ||
239 | list<Int_t>::iterator tClustIter; | |
240 | ||
241 | Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster]; | |
242 | Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices); | |
243 | for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) { | |
244 | if (tTrackIndices[tNcl] >= 0) { | |
245 | tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]); | |
246 | if (tClustIter == fClusterPerPadrow[tNcl]->end()) { | |
247 | fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]); | |
248 | } | |
249 | else { | |
250 | fSharedList[tNcl]->push_back(tTrackIndices[tNcl]); | |
251 | } | |
252 | } | |
253 | } | |
254 | ||
255 | } | |
256 | ||
257 | for (int i=0;i<nofTracks;i++) | |
258 | { | |
259 | bool good_momentum=true; //flaga to chcek if we can read momentum of this track | |
260 | ||
261 | AliFemtoTrack* trackCopy = new AliFemtoTrack(); | |
262 | const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track | |
263 | const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack(); | |
264 | ||
265 | trackCopy->SetCharge((short)esdtrack->GetSign()); | |
266 | ||
267 | //in aliroot we have AliPID | |
268 | //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon | |
269 | //we use only 5 first | |
270 | double esdpid[5]; | |
271 | esdtrack->GetESDpid(esdpid); | |
272 | trackCopy->SetPidProbElectron(esdpid[0]); | |
273 | trackCopy->SetPidProbMuon(esdpid[1]); | |
274 | trackCopy->SetPidProbPion(esdpid[2]); | |
275 | trackCopy->SetPidProbKaon(esdpid[3]); | |
276 | trackCopy->SetPidProbProton(esdpid[4]); | |
277 | ||
278 | double pxyz[3]; | |
279 | if (fConstrained==true) | |
280 | good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum | |
281 | else | |
282 | good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum | |
283 | AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]); | |
284 | trackCopy->SetP(v);//setting momentum | |
285 | trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1])); | |
286 | const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]); | |
287 | if (p.mag() == 0) { | |
288 | delete trackCopy; | |
289 | continue; | |
290 | } | |
291 | const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]); | |
292 | //setting helix I do not if it is ok | |
293 | AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); | |
294 | trackCopy->SetHelix(helix); | |
295 | ||
296 | trackCopy->SetTrackId(esdtrack->GetID()); | |
297 | trackCopy->SetFlags(esdtrack->GetStatus()); | |
298 | //trackCopy->SetLabel(esdtrack->GetLabel()); | |
299 | ||
300 | //some stuff which could be useful | |
301 | float impact[2]; | |
302 | float covimpact[3]; | |
303 | esdtrack->GetImpactParameters(impact,covimpact); | |
304 | trackCopy->SetImpactD(impact[0]); | |
305 | trackCopy->SetImpactZ(impact[1]); | |
306 | trackCopy->SetCdd(covimpact[0]); | |
307 | trackCopy->SetCdz(covimpact[1]); | |
308 | trackCopy->SetCzz(covimpact[2]); | |
309 | trackCopy->SetITSchi2(esdtrack->GetITSchi2()); | |
310 | trackCopy->SetITSncls(esdtrack->GetNcls(0)); | |
311 | trackCopy->SetTPCchi2(esdtrack->GetTPCchi2()); | |
312 | trackCopy->SetTPCncls(esdtrack->GetTPCNcls()); | |
313 | trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF()); | |
314 | trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class | |
315 | trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); | |
316 | ||
317 | // Fill cluster per padrow information | |
318 | Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster]; | |
319 | Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices); | |
320 | for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) { | |
321 | if (tTrackIndices[tNcl] > 0) | |
322 | trackCopy->SetTPCcluster(tNcl, 1); | |
323 | else | |
324 | trackCopy->SetTPCcluster(tNcl, 0); | |
325 | } | |
326 | ||
327 | // Fill shared cluster information | |
328 | list<Int_t>::iterator tClustIter; | |
329 | ||
330 | for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) { | |
331 | if (tTrackIndices[tNcl] > 0) { | |
332 | tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]); | |
333 | if (tClustIter != fSharedList[tNcl]->end()) { | |
334 | trackCopy->SetTPCshared(tNcl, 1); | |
335 | cout << "Event next" << endl; | |
336 | cout << "Track: " << i << endl; | |
337 | cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl; | |
338 | } | |
339 | else { | |
340 | trackCopy->SetTPCshared(tNcl, 0); | |
341 | } | |
342 | } | |
343 | } | |
344 | ||
345 | //decision if we want this track | |
346 | //if we using diffrent labels we want that this label was use for first time | |
347 | //if we use hidden info we want to have match between sim data and ESD | |
348 | if (good_momentum==true) | |
349 | { | |
350 | hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis | |
351 | realnofTracks++;//real number of tracks | |
352 | } | |
353 | else | |
354 | { | |
355 | delete trackCopy; | |
356 | } | |
357 | ||
358 | } | |
359 | ||
360 | hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event | |
361 | fCurEvent++; | |
362 | cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl; | |
363 | if (fCurEvent== fNumberofEvent)//if end of current file close all | |
364 | { | |
365 | fTree->Reset(); | |
366 | delete fTree; | |
367 | fEsdFile->Close(); | |
368 | } | |
369 | return hbtEvent; | |
370 | } | |
371 | ||
372 | ||
373 | ||
374 | ||
375 | ||
376 | ||
377 | ||
378 | ||
379 |