]>
Commit | Line | Data |
---|---|---|
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 | ||
257 | iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848); | |
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 |