Adaption to new fluka common blocks (E. Futo)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv2.cxx
CommitLineData
1b446896 1#include "AliHBTReaderITSv2.h"
2
d0c23b58 3#include <Riostream.h>
a9bfdd7b 4#include <TString.h>
5#include <TObjString.h>
6#include <TTree.h>
7#include <TFile.h>
8#include <TParticle.h>
9
10#include <AliRun.h>
11#include <AliMagF.h>
12#include <AliITStrackV2.h>
a9bfdd7b 13#include <AliITStrackerV2.h>
14#include <AliITSgeom.h>
15
16#include "AliHBTRun.h"
17#include "AliHBTEvent.h"
18#include "AliHBTParticle.h"
19#include "AliHBTParticleCut.h"
20
21
1b446896 22ClassImp(AliHBTReaderITSv2)
23
98295f4b 24AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* trackfilename,
25 const Char_t* clusterfilename,
26 const Char_t* galicefilename):
27 fTrackFileName(trackfilename),
28 fClusterFileName(clusterfilename),
29 fGAliceFileName(galicefilename)
a9bfdd7b 30{
31 //constructor,
32 //Defaults:
33 // trackfilename = "AliITStracksV2.root"
34 // clusterfilename = "AliITSclustersV2.root"
35 // galicefilename = "galice.root"
36 fParticles = new AliHBTRun();
37 fTracks = new AliHBTRun();
38 fIsRead = kFALSE;
39}
40/********************************************************************/
41
98295f4b 42AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs,
43 const Char_t* trackfilename,
44 const Char_t* clusterfilename,
45 const Char_t* galicefilename):
46 AliHBTReader(dirs),
47 fTrackFileName(trackfilename),
48 fClusterFileName(clusterfilename),
49 fGAliceFileName(galicefilename)
a9bfdd7b 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
63AliHBTReaderITSv2::~AliHBTReaderITSv2()
64 {
65 if (fParticles) delete fParticles;
66 if (fTracks) delete fTracks;
67 }
68/********************************************************************/
69/********************************************************************/
70
71AliHBTEvent* 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
85AliHBTEvent* 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
98Int_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/********************************************************************/
111Int_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/********************************************************************/
126Int_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
98295f4b 148 AliITStrackV2 *iotrack= 0x0; //buffer track for reading data from tree
a9bfdd7b 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);
5ee8b891 177 currentdir++;
178 continue;
a9bfdd7b 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
98295f4b 184 Info("Read","________________________________________________________");
185 Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
186 Float_t mf;
187 if (fUseMagFFromRun)
188 {
189 mf = gAlice->Field()->SolenoidField();
190 Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
191 }
192 else
193 {
194 Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
195 mf = fMagneticField*10.;
196 }
197 AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
198 if (iotrack == 0x0) iotrack = new AliITStrackV2();
a9bfdd7b 199 }
200 else
201 {//if not return an error
202 Error("Read","Can not find Header tree (TreeE) in gAlice");
5ee8b891 203 currentdir++;
204 continue;
a9bfdd7b 205 }
206
207 AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom");
208 if (!geom)
209 {
210 Error("Read","Can't get the ITS geometry!");
5ee8b891 211 currentdir++;
212 continue;
a9bfdd7b 213 }
214
215 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
216 {
217 cout<<"Reading Event "<<currentEvent<<endl;
218
219 aGAliceFile->cd();
220 gAlice->GetEvent(currentEvent);
9f69cb10 221
a9bfdd7b 222 aClustersFile->cd();
a9bfdd7b 223 sprintf(tname,"TreeT_ITS_%d",currentEvent);
224
225 tracktree=(TTree*)aTracksFile->Get(tname);
226 if (!tracktree)
227 {
228 Error("Read","Can't get a tree with ITS tracks");
5ee8b891 229 continue;
a9bfdd7b 230 }
231 TBranch *tbranch=tracktree->GetBranch("tracks");
232
233 Ntracks=(Int_t)tracktree->GetEntries();
234
235 Int_t accepted = 0;
236 Int_t tpcfault = 0;
237 Int_t itsfault = 0;
238 for (i=0; i<Ntracks; i++) //loop over all tpc tracks
239 {
9f69cb10 240 if(i%100 == 0)cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<"\r";
a9bfdd7b 241
242 tbranch->SetAddress(&iotrack);
243 tracktree->GetEvent(i);
244
245 label=iotrack->GetLabel();
246 if (label < 0)
247 {
248 tpcfault++;
249 continue;
250 }
a9bfdd7b 251
252 TParticle *p = (TParticle*)gAlice->Particle(label);
f1e2da22 253 if(p == 0x0) continue; //if returned pointer is NULL
254 if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
255
a9bfdd7b 256 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
257 //if not take next partilce
258
8089a083 259 AliHBTParticle* part = new AliHBTParticle(*p,i);
a9bfdd7b 260 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
261 //if it does not delete it and take next good track
262
263
5f119c5b 264 iotrack->PropagateTo(3.,0.0028,65.19);
a9bfdd7b 265 iotrack->PropagateToVertex();
266
267 iotrack->GetExternalParameters(xk,par); //get properties of the track
268 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
269 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
270 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
271 lam=par[3];
272 pt=1.0/TMath::Abs(par[4]);
273
274 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
275 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
276 Double_t tpz = pt * lam; //track z coordinate of momentum
277
278 Double_t mass = p->GetMass();
279 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
280
8089a083 281 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
a9bfdd7b 282 if(Pass(track))//check if meets all criteria of any of our cuts
283 //if it does not delete it and take next good track
284 {
285 delete track;
286 delete part;
287 continue;
288 }
289 particles->AddParticle(totalNevents,part);//put track and particle on the run
290 tracks->AddParticle(totalNevents,track);
291 accepted++;
292 }//end of loop over tracks in the event
293
294 aTracksFile->Delete(tname);
295 aTracksFile->Delete("tracks");
41515b05 296// delete tracker;
a9bfdd7b 297
298 totalNevents++;
a9bfdd7b 299 cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<endl;
300
301 }//end of loop over events in current directory
7836ee94 302 CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
303 currentdir++;
a9bfdd7b 304 }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
305
306 delete iotrack;
307 fIsRead = kTRUE;
308 return 0;
309}
310
311/********************************************************************/
312Int_t AliHBTReaderITSv2::OpenFiles
313(TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
314{
315 //opens all the files
316
317
318 const TString& dirname = GetDirName(event);
319 if (dirname == "")
320 {
321 Error("OpenFiles","Can not get directory name");
322 return 4;
323 }
324
325 TString filename = dirname +"/"+ fTrackFileName;
326 aTracksFile = TFile::Open(filename.Data());
327 if ( aTracksFile == 0x0 )
328 {
329 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
330 return 1;
331 }
332 if (!aTracksFile->IsOpen())
333 {
334 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
335 return 1;
336 }
337
338 filename = dirname +"/"+ fClusterFileName;
339 aClustersFile = TFile::Open(filename.Data());
340 if ( aClustersFile == 0x0 )
341 {
342 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
343 return 2;
344 }
345 if (!aClustersFile->IsOpen())
346 {
347 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
348 return 2;
349 }
350
351 filename = dirname +"/"+ fGAliceFileName;
352 agAliceFile = TFile::Open(filename.Data());
353 if ( agAliceFile== 0x0)
354 {
355 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
356 return 3;
357 }
358 if (!agAliceFile->IsOpen())
359 {
360 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
361 return 3;
362 }
363
364 if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
365 {
366 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
367 return 5;
368 }
369
370 return 0;
371}
372/********************************************************************/
373
374/********************************************************************/
375
376void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
377{
378 //closes the files
379 tracksFile->Close();
380 delete tracksFile;
381 tracksFile = 0x0;
382
383 clustersFile->Close();
384 delete clustersFile;
385 clustersFile = 0x0;
386
387 delete gAlice;
388 gAlice = 0;
389
390 if (gAliceFile)
391 {
392 gAliceFile->Close();
393 delete gAliceFile;
394 gAliceFile = 0x0;
395 }
396}
397
398/********************************************************************/
399
400