Dummy copy constructor
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv2.cxx
CommitLineData
88cb7938 1
1b446896 2#include "AliHBTReaderITSv2.h"
3
88cb7938 4#include <Riostream.h>
d0c23b58 5#include <Riostream.h>
a9bfdd7b 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>
88cb7938 13#include <AliRunLoader.h>
14#include <AliLoader.h>
a9bfdd7b 15#include <AliMagF.h>
88cb7938 16#include <AliITS.h>
a9bfdd7b 17#include <AliITStrackV2.h>
a9bfdd7b 18#include <AliITStrackerV2.h>
19#include <AliITSgeom.h>
20
21#include "AliHBTRun.h"
22#include "AliHBTEvent.h"
23#include "AliHBTParticle.h"
24#include "AliHBTParticleCut.h"
25
26
1b446896 27ClassImp(AliHBTReaderITSv2)
28
88cb7938 29AliHBTReaderITSv2::AliHBTReaderITSv2():fFileName("galice.root")
30{
31 //constructor,
32 //Defaults:
33 // galicefilename = "galice.root"
34 fParticles = 0x0;
35 fTracks = 0x0;
36 fIsRead = kFALSE;
37}
38/********************************************************************/
39
40AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):fFileName(galicefilename)
a9bfdd7b 41{
42 //constructor,
43 //Defaults:
a9bfdd7b 44 // galicefilename = "galice.root"
45 fParticles = new AliHBTRun();
46 fTracks = new AliHBTRun();
47 fIsRead = kFALSE;
48}
49/********************************************************************/
50
88cb7938 51AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename):
52 AliHBTReader(dirs), fFileName(galicefilename)
a9bfdd7b 53{
54 //constructor,
55 //Defaults:
a9bfdd7b 56 // galicefilename = "galice.root"
57
58 fParticles = new AliHBTRun();
59 fTracks = new AliHBTRun();
60 fIsRead = kFALSE;
61}
62/********************************************************************/
63
64AliHBTReaderITSv2::~AliHBTReaderITSv2()
65 {
66 if (fParticles) delete fParticles;
67 if (fTracks) delete fTracks;
68 }
69/********************************************************************/
70/********************************************************************/
71
72AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
73 {
74 //returns Nth event with simulated particles
88cb7938 75 if (!fIsRead)
76 {
77 if (fParticles == 0x0) fParticles = new AliHBTRun();
78 if (fTracks == 0x0) fTracks = new AliHBTRun();
a9bfdd7b 79 if(Read(fParticles,fTracks))
80 {
81 Error("GetParticleEvent","Error in reading");
82 return 0x0;
83 }
88cb7938 84 }
85 return (fParticles)?fParticles->GetEvent(n):0x0;
86}
a9bfdd7b 87/********************************************************************/
88
89AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
90 {
91 //returns Nth event with reconstructed tracks
88cb7938 92 if (!fIsRead)
93 {
94 if (fParticles == 0x0) fParticles = new AliHBTRun();
95 if (fTracks == 0x0) fTracks = new AliHBTRun();
a9bfdd7b 96 if(Read(fParticles,fTracks))
97 {
98 Error("GetTrackEvent","Error in reading");
99 return 0x0;
100 }
88cb7938 101 }
102 return (fTracks)?fTracks->GetEvent(n):0x0;
a9bfdd7b 103 }
104/********************************************************************/
105
106Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
107 {
108 //returns number of events of particles
88cb7938 109 if (!fIsRead)
110 {
111 if (fParticles == 0x0) fParticles = new AliHBTRun();
112 if (fTracks == 0x0) fTracks = new AliHBTRun();
a9bfdd7b 113 if(Read(fParticles,fTracks))
114 {
115 Error("GetNumberOfPartEvents","Error in reading");
116 return 0;
117 }
88cb7938 118 }
119 return (fParticles)?fParticles->GetNumberOfEvents():0;
a9bfdd7b 120 }
121
122/********************************************************************/
123Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
124 {
125 //returns number of events of tracks
126 if (!fIsRead)
88cb7938 127 {
128 if (fParticles == 0x0) fParticles = new AliHBTRun();
129 if (fTracks == 0x0) fTracks = new AliHBTRun();
a9bfdd7b 130 if(Read(fParticles,fTracks))
131 {
132 Error("GetNumberOfTrackEvents","Error in reading");
133 return 0;
134 }
88cb7938 135 }
136 return (fTracks)?fTracks->GetNumberOfEvents():0;
a9bfdd7b 137 }
a9bfdd7b 138/********************************************************************/
139/********************************************************************/
88cb7938 140
141
a9bfdd7b 142Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
143{
88cb7938 144//reads data
a9bfdd7b 145 Int_t Nevents = 0; //number of events found in given directory
146 Int_t Ndirs; //number of the directories to be read
147 Int_t Ntracks; //number of tracks in current event
148 Int_t currentdir = 0; //number of events in the current directory
149 Int_t totalNevents = 0; //total number of events read from all directories up to now
88cb7938 150 register Int_t i = 0; //iterator
a9bfdd7b 151
41515b05 152// AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
a9bfdd7b 153 TTree *tracktree; // tree for tracks
154
155 Double_t xk;
156 Double_t par[5]; //Kalman track parameters
157 Float_t phi, lam, pt;//angles and transverse momentum
158 Int_t label; //label of the current track
159
88cb7938 160 AliITStrackV2 *iotrack = 0x0; //buffer track for reading data from tree
a9bfdd7b 161
162 if (!particles) //check if an object is instatiated
163 {
164 Error("Read"," particles object must instatiated before passing it to the reader");
165 }
166 if (!tracks) //check if an object is instatiated
167 {
168 Error("Read"," tracks object must instatiated before passing it to the reader");
169 }
170 particles->Reset();//clear runs == delete all old events
171 tracks->Reset();
172
173 if (fDirs) //if array with directories is supplied by user
174 {
175 Ndirs = fDirs->GetEntries(); //get the number if directories
176 }
177 else
178 {
179 Ndirs = 0; //if the array is not supplied read only from current directory
180 }
181
182// cout<<"Found "<<Ndirs<<" directory entries"<<endl;
183
184 do //do while is good even if Ndirs==0 (than read from current directory)
185 {
88cb7938 186 TString filename = GetDirName(currentdir);
187 if (filename.IsNull())
a9bfdd7b 188 {
88cb7938 189 Error("Read","Can not get directory name");
5ee8b891 190 currentdir++;
191 continue;
a9bfdd7b 192 }
88cb7938 193 filename = filename +"/"+ fFileName;
194 AliRunLoader* rl = AliRunLoader::Open(filename);
195 if( rl == 0x0)
196 {
197 Error("Read","Exiting due to problems with opening files.");
198 currentdir++;
199 continue;
200 }
201
202 rl->LoadHeader();
203 rl->LoadKinematics();
204 rl->LoadgAlice();
205 gAlice = rl->GetAliRun();
206 AliITS* its = (AliITS*)gAlice->GetModule("ITS");
207
208 AliLoader* itsl = rl->GetLoader("ITSLoader");
a9bfdd7b 209
88cb7938 210 if ((its == 0x0) || ( itsl== 0x0))
211 {
212 Error("Read","Can not found ITS in this run");
213 delete rl;
214 rl = 0x0;
215 currentdir++;
216 continue;
217 }
218 Nevents = rl->GetNumberOfEvents();
219
220 if (Nevents > 0)//check if tree E exists
a9bfdd7b 221 {
98295f4b 222 Info("Read","________________________________________________________");
223 Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
224 Float_t mf;
225 if (fUseMagFFromRun)
226 {
227 mf = gAlice->Field()->SolenoidField();
228 Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
229 }
230 else
231 {
232 Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
233 mf = fMagneticField*10.;
234 }
235 AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
236 if (iotrack == 0x0) iotrack = new AliITStrackV2();
a9bfdd7b 237 }
238 else
239 {//if not return an error
88cb7938 240 Error("Read","No events in this run");
241 delete rl;
242 rl = 0x0;
5ee8b891 243 currentdir++;
244 continue;
a9bfdd7b 245 }
246
88cb7938 247 AliITSgeom *geom= its->GetITSgeom();
a9bfdd7b 248 if (!geom)
249 {
250 Error("Read","Can't get the ITS geometry!");
88cb7938 251 delete rl;
252 rl = 0x0;
5ee8b891 253 currentdir++;
254 continue;
a9bfdd7b 255 }
256
88cb7938 257 itsl->LoadTracks();
258
a9bfdd7b 259 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
260 {
261 cout<<"Reading Event "<<currentEvent<<endl;
88cb7938 262 rl->GetEvent(currentEvent);
263 tracktree=itsl->TreeT();
a9bfdd7b 264
a9bfdd7b 265 if (!tracktree)
266 {
267 Error("Read","Can't get a tree with ITS tracks");
5ee8b891 268 continue;
a9bfdd7b 269 }
270 TBranch *tbranch=tracktree->GetBranch("tracks");
a9bfdd7b 271 Ntracks=(Int_t)tracktree->GetEntries();
272
273 Int_t accepted = 0;
274 Int_t tpcfault = 0;
275 Int_t itsfault = 0;
276 for (i=0; i<Ntracks; i++) //loop over all tpc tracks
277 {
9f69cb10 278 if(i%100 == 0)cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<"\r";
a9bfdd7b 279
280 tbranch->SetAddress(&iotrack);
281 tracktree->GetEvent(i);
282
283 label=iotrack->GetLabel();
284 if (label < 0)
285 {
286 tpcfault++;
287 continue;
288 }
a9bfdd7b 289
290 TParticle *p = (TParticle*)gAlice->Particle(label);
f1e2da22 291 if(p == 0x0) continue; //if returned pointer is NULL
292 if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
293
a9bfdd7b 294 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
295 //if not take next partilce
296
8089a083 297 AliHBTParticle* part = new AliHBTParticle(*p,i);
a9bfdd7b 298 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
299 //if it does not delete it and take next good track
300
5f119c5b 301 iotrack->PropagateTo(3.,0.0028,65.19);
a9bfdd7b 302 iotrack->PropagateToVertex();
303
304 iotrack->GetExternalParameters(xk,par); //get properties of the track
305 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
306 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
307 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
308 lam=par[3];
309 pt=1.0/TMath::Abs(par[4]);
310
311 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
312 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
313 Double_t tpz = pt * lam; //track z coordinate of momentum
314
315 Double_t mass = p->GetMass();
316 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
317
8089a083 318 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
a9bfdd7b 319 if(Pass(track))//check if meets all criteria of any of our cuts
320 //if it does not delete it and take next good track
321 {
322 delete track;
323 delete part;
324 continue;
325 }
326 particles->AddParticle(totalNevents,part);//put track and particle on the run
327 tracks->AddParticle(totalNevents,track);
328 accepted++;
329 }//end of loop over tracks in the event
a9bfdd7b 330
331 totalNevents++;
a9bfdd7b 332 cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<endl;
333
334 }//end of loop over events in current directory
88cb7938 335 delete rl;
7836ee94 336 currentdir++;
a9bfdd7b 337 }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
338
339 delete iotrack;
340 fIsRead = kTRUE;
341 return 0;
342}
343
344/********************************************************************/
a9bfdd7b 345/********************************************************************/
346
347