1b446896 |
1 | |
2 | #include "AliHBTReaderITSv2.h" |
3 | |
a9bfdd7b |
4 | #include <iostream.h> |
5 | #include <fstream.h> |
6 | #include <TString.h> |
7 | #include <TObjString.h> |
8 | #include <TTree.h> |
9 | #include <TFile.h> |
10 | #include <TParticle.h> |
11 | |
12 | #include <AliRun.h> |
13 | #include <AliMagF.h> |
14 | #include <AliITStrackV2.h> |
15 | //#include <AliITSParam.h> |
16 | #include <AliITStrackerV2.h> |
17 | #include <AliITSgeom.h> |
18 | |
19 | #include "AliHBTRun.h" |
20 | #include "AliHBTEvent.h" |
21 | #include "AliHBTParticle.h" |
22 | #include "AliHBTParticleCut.h" |
23 | |
24 | |
1b446896 |
25 | ClassImp(AliHBTReaderITSv2) |
26 | |
a9bfdd7b |
27 | AliHBTReaderITSv2:: |
28 | AliHBTReaderITSv2(const Char_t* trackfilename, const Char_t* clusterfilename, |
29 | const Char_t* galicefilename) |
30 | :fTrackFileName(trackfilename),fClusterFileName(clusterfilename), |
31 | fGAliceFileName(galicefilename) |
32 | { |
33 | //constructor, |
34 | //Defaults: |
35 | // trackfilename = "AliITStracksV2.root" |
36 | // clusterfilename = "AliITSclustersV2.root" |
37 | // galicefilename = "galice.root" |
38 | fParticles = new AliHBTRun(); |
39 | fTracks = new AliHBTRun(); |
40 | fIsRead = kFALSE; |
41 | } |
42 | /********************************************************************/ |
43 | |
44 | AliHBTReaderITSv2:: |
45 | AliHBTReaderITSv2(TObjArray* dirs, const Char_t* trackfilename, |
46 | const Char_t* clusterfilename, const Char_t* galicefilename) |
47 | : AliHBTReader(dirs), |
48 | fTrackFileName(trackfilename),fClusterFileName(clusterfilename), |
49 | fGAliceFileName(galicefilename) |
50 | { |
51 | //constructor, |
52 | //Defaults: |
53 | // trackfilename = "AliITStracksV2.root" |
54 | // clusterfilename = "AliITSclustersV2.root" |
55 | // galicefilename = "galice.root" |
56 | |
57 | fParticles = new AliHBTRun(); |
58 | fTracks = new AliHBTRun(); |
59 | fIsRead = kFALSE; |
60 | } |
61 | /********************************************************************/ |
62 | |
63 | AliHBTReaderITSv2::~AliHBTReaderITSv2() |
64 | { |
65 | if (fParticles) delete fParticles; |
66 | if (fTracks) delete fTracks; |
67 | } |
68 | /********************************************************************/ |
69 | /********************************************************************/ |
70 | |
71 | AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n) |
72 | { |
73 | //returns Nth event with simulated particles |
74 | if (!fIsRead) |
75 | if(Read(fParticles,fTracks)) |
76 | { |
77 | Error("GetParticleEvent","Error in reading"); |
78 | return 0x0; |
79 | } |
80 | |
81 | return fParticles->GetEvent(n); |
82 | } |
83 | /********************************************************************/ |
84 | |
85 | AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n) |
86 | { |
87 | //returns Nth event with reconstructed tracks |
88 | if (!fIsRead) |
89 | if(Read(fParticles,fTracks)) |
90 | { |
91 | Error("GetTrackEvent","Error in reading"); |
92 | return 0x0; |
93 | } |
94 | return fTracks->GetEvent(n); |
95 | } |
96 | /********************************************************************/ |
97 | |
98 | Int_t AliHBTReaderITSv2::GetNumberOfPartEvents() |
99 | { |
100 | //returns number of events of particles |
101 | if (!fIsRead) |
102 | if(Read(fParticles,fTracks)) |
103 | { |
104 | Error("GetNumberOfPartEvents","Error in reading"); |
105 | return 0; |
106 | } |
107 | return fParticles->GetNumberOfEvents(); |
108 | } |
109 | |
110 | /********************************************************************/ |
111 | Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents() |
112 | { |
113 | //returns number of events of tracks |
114 | if (!fIsRead) |
115 | if(Read(fParticles,fTracks)) |
116 | { |
117 | Error("GetNumberOfTrackEvents","Error in reading"); |
118 | return 0; |
119 | } |
120 | return fTracks->GetNumberOfEvents(); |
121 | } |
122 | |
123 | |
124 | /********************************************************************/ |
125 | /********************************************************************/ |
126 | Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks) |
127 | { |
128 | Int_t Nevents = 0; //number of events found in given directory |
129 | Int_t Ndirs; //number of the directories to be read |
130 | Int_t Ntracks; //number of tracks in current event |
131 | Int_t currentdir = 0; //number of events in the current directory |
132 | Int_t totalNevents = 0; //total number of events read from all directories up to now |
133 | register Int_t i; //iterator |
134 | |
135 | TFile *aTracksFile;//file with tracks |
136 | TFile *aClustersFile;//file with clusters |
137 | TFile *aGAliceFile;//file name with galice |
138 | |
41515b05 |
139 | // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels |
a9bfdd7b |
140 | TTree *tracktree; // tree for tracks |
141 | |
142 | Double_t xk; |
143 | Double_t par[5]; //Kalman track parameters |
144 | Float_t phi, lam, pt;//angles and transverse momentum |
145 | Int_t label; //label of the current track |
146 | |
147 | char tname[100]; //buffer for tree name |
148 | AliITStrackV2 *iotrack= new AliITStrackV2(); //buffer track for reading data from tree |
149 | |
150 | if (!particles) //check if an object is instatiated |
151 | { |
152 | Error("Read"," particles object must instatiated before passing it to the reader"); |
153 | } |
154 | if (!tracks) //check if an object is instatiated |
155 | { |
156 | Error("Read"," tracks object must instatiated before passing it to the reader"); |
157 | } |
158 | particles->Reset();//clear runs == delete all old events |
159 | tracks->Reset(); |
160 | |
161 | if (fDirs) //if array with directories is supplied by user |
162 | { |
163 | Ndirs = fDirs->GetEntries(); //get the number if directories |
164 | } |
165 | else |
166 | { |
167 | Ndirs = 0; //if the array is not supplied read only from current directory |
168 | } |
169 | |
170 | // cout<<"Found "<<Ndirs<<" directory entries"<<endl; |
171 | |
172 | do //do while is good even if Ndirs==0 (than read from current directory) |
173 | { |
174 | if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) ) |
175 | { |
176 | Error("Read","Exiting due to problems with opening files. Errorcode %d",i); |
177 | delete iotrack; |
178 | return i; |
179 | } |
180 | |
181 | if (gAlice->TreeE())//check if tree E exists |
182 | { |
183 | Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice |
184 | cout<<"________________________________________________________\n"; |
185 | cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl; |
186 | cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl; |
187 | AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor()); |
188 | } |
189 | else |
190 | {//if not return an error |
191 | Error("Read","Can not find Header tree (TreeE) in gAlice"); |
192 | delete iotrack; |
193 | return 1; |
194 | } |
195 | |
196 | AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom"); |
197 | if (!geom) |
198 | { |
199 | Error("Read","Can't get the ITS geometry!"); |
200 | delete iotrack; |
201 | return 2; |
202 | } |
203 | |
204 | for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events |
205 | { |
206 | cout<<"Reading Event "<<currentEvent<<endl; |
207 | |
208 | aGAliceFile->cd(); |
209 | gAlice->GetEvent(currentEvent); |
9f69cb10 |
210 | |
211 | // TParticle * part = gAlice->Particle(0); |
a9bfdd7b |
212 | |
213 | aClustersFile->cd(); |
9f69cb10 |
214 | // Double_t orz=part->Vz(); |
41515b05 |
215 | // tracker = new AliITStrackerV2(geom,currentEvent,orz); //<---- this is for Massimo version |
216 | // tracker = new AliITStrackerV2(geom,currentEvent); |
a9bfdd7b |
217 | sprintf(tname,"TreeT_ITS_%d",currentEvent); |
218 | |
219 | tracktree=(TTree*)aTracksFile->Get(tname); |
220 | if (!tracktree) |
221 | { |
222 | Error("Read","Can't get a tree with ITS tracks"); |
223 | delete iotrack; |
41515b05 |
224 | // delete tracker; |
a9bfdd7b |
225 | return 4; |
226 | } |
227 | TBranch *tbranch=tracktree->GetBranch("tracks"); |
228 | |
229 | Ntracks=(Int_t)tracktree->GetEntries(); |
230 | |
231 | Int_t accepted = 0; |
232 | Int_t tpcfault = 0; |
233 | Int_t itsfault = 0; |
234 | for (i=0; i<Ntracks; i++) //loop over all tpc tracks |
235 | { |
9f69cb10 |
236 | if(i%100 == 0)cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<"\r"; |
a9bfdd7b |
237 | |
238 | tbranch->SetAddress(&iotrack); |
239 | tracktree->GetEvent(i); |
240 | |
241 | label=iotrack->GetLabel(); |
242 | if (label < 0) |
243 | { |
244 | tpcfault++; |
245 | continue; |
246 | } |
a9bfdd7b |
247 | |
248 | TParticle *p = (TParticle*)gAlice->Particle(label); |
249 | if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type |
250 | //if not take next partilce |
251 | |
252 | AliHBTParticle* part = new AliHBTParticle(*p); |
253 | if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts |
254 | //if it does not delete it and take next good track |
255 | |
256 | |
5f119c5b |
257 | iotrack->PropagateTo(3.,0.0028,65.19); |
a9bfdd7b |
258 | iotrack->PropagateToVertex(); |
259 | |
260 | iotrack->GetExternalParameters(xk,par); //get properties of the track |
261 | phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); |
262 | if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); |
263 | if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); |
264 | lam=par[3]; |
265 | pt=1.0/TMath::Abs(par[4]); |
266 | |
267 | Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum |
268 | Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum |
269 | Double_t tpz = pt * lam; //track z coordinate of momentum |
270 | |
271 | Double_t mass = p->GetMass(); |
272 | Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track |
273 | |
274 | AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.); |
275 | if(Pass(track))//check if meets all criteria of any of our cuts |
276 | //if it does not delete it and take next good track |
277 | { |
278 | delete track; |
279 | delete part; |
280 | continue; |
281 | } |
282 | particles->AddParticle(totalNevents,part);//put track and particle on the run |
283 | tracks->AddParticle(totalNevents,track); |
284 | accepted++; |
285 | }//end of loop over tracks in the event |
286 | |
287 | aTracksFile->Delete(tname); |
288 | aTracksFile->Delete("tracks"); |
41515b05 |
289 | // delete tracker; |
a9bfdd7b |
290 | |
291 | totalNevents++; |
a9bfdd7b |
292 | cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<endl; |
293 | |
294 | }//end of loop over events in current directory |
7836ee94 |
295 | CloseFiles(aTracksFile,aClustersFile,aGAliceFile); |
296 | currentdir++; |
a9bfdd7b |
297 | }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array |
298 | |
299 | delete iotrack; |
300 | fIsRead = kTRUE; |
301 | return 0; |
302 | } |
303 | |
304 | /********************************************************************/ |
305 | Int_t AliHBTReaderITSv2::OpenFiles |
306 | (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event) |
307 | { |
308 | //opens all the files |
309 | |
310 | |
311 | const TString& dirname = GetDirName(event); |
312 | if (dirname == "") |
313 | { |
314 | Error("OpenFiles","Can not get directory name"); |
315 | return 4; |
316 | } |
317 | |
318 | TString filename = dirname +"/"+ fTrackFileName; |
319 | aTracksFile = TFile::Open(filename.Data()); |
320 | if ( aTracksFile == 0x0 ) |
321 | { |
322 | Error("OpenFiles","Can't open file with tacks named %s",filename.Data()); |
323 | return 1; |
324 | } |
325 | if (!aTracksFile->IsOpen()) |
326 | { |
327 | Error("OpenFiles","Can't open file with tacks named %s",filename.Data()); |
328 | return 1; |
329 | } |
330 | |
331 | filename = dirname +"/"+ fClusterFileName; |
332 | aClustersFile = TFile::Open(filename.Data()); |
333 | if ( aClustersFile == 0x0 ) |
334 | { |
335 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); |
336 | return 2; |
337 | } |
338 | if (!aClustersFile->IsOpen()) |
339 | { |
340 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); |
341 | return 2; |
342 | } |
343 | |
344 | filename = dirname +"/"+ fGAliceFileName; |
345 | agAliceFile = TFile::Open(filename.Data()); |
346 | if ( agAliceFile== 0x0) |
347 | { |
348 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); |
349 | return 3; |
350 | } |
351 | if (!agAliceFile->IsOpen()) |
352 | { |
353 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); |
354 | return 3; |
355 | } |
356 | |
357 | if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) |
358 | { |
359 | Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data()); |
360 | return 5; |
361 | } |
362 | |
363 | return 0; |
364 | } |
365 | /********************************************************************/ |
366 | |
367 | /********************************************************************/ |
368 | |
369 | void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile) |
370 | { |
371 | //closes the files |
372 | tracksFile->Close(); |
373 | delete tracksFile; |
374 | tracksFile = 0x0; |
375 | |
376 | clustersFile->Close(); |
377 | delete clustersFile; |
378 | clustersFile = 0x0; |
379 | |
380 | delete gAlice; |
381 | gAlice = 0; |
382 | |
383 | if (gAliceFile) |
384 | { |
385 | gAliceFile->Close(); |
386 | delete gAliceFile; |
387 | gAliceFile = 0x0; |
388 | } |
389 | } |
390 | |
391 | /********************************************************************/ |
392 | |
393 | |