Bug fix
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv1.cxx
CommitLineData
1b446896 1
64e553ff 2#include "AliHBTReaderITSv1.h"
3#include "AliHBTEvent.h"
4#include "AliHBTRun.h"
5#include "AliHBTParticle.h"
6#include "AliHBTParticleCut.h"
1b446896 7
64e553ff 8#include <iostream.h>
9
10#include <TROOT.h>
11#include <TFile.h>
12#include <TTree.h>
13#include <TBranch.h>
14#include <TObjArray.h>
15#include <TParticle.h>
16#include <TString.h>
17#include <TObjString.h>
18
19#include <AliRun.h>
20#include <AliMagF.h>
21#include <AliKalmanTrack.h>
22#include <AliITSIOTrack.h>
1b446896 23
64e553ff 24ClassImp(AliHBTReaderITSv1)
25/********************************************************************/
1b446896 26
64e553ff 27AliHBTReaderITSv1::
28AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
29 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
1b446896 30 {
31 fParticles = new AliHBTRun();
32 fTracks = new AliHBTRun();
33 fIsRead = kFALSE;
34 }
64e553ff 35/********************************************************************/
36
37AliHBTReaderITSv1::
38AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
39 AliHBTReader(dirs),
40 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
41 {
42 fParticles = new AliHBTRun();
43 fTracks = new AliHBTRun();
44 fIsRead = kFALSE;
45 }
46/********************************************************************/
47
48AliHBTReaderITSv1::~AliHBTReaderITSv1()
1b446896 49{
50 delete fParticles;
51 delete fTracks;
52}
64e553ff 53/********************************************************************/
54
1b446896 55AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
56 {
57 //returns Nth event with simulated particles
64e553ff 58 if (!fIsRead)
59 if(Read(fParticles,fTracks))
60 {
61 Error("GetParticleEvent","Error in reading");
62 return 0x0;
63 }
64
1b446896 65 return fParticles->GetEvent(n);
66 }
67/********************************************************************/
64e553ff 68
1b446896 69AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
70 {
71 //returns Nth event with reconstructed tracks
64e553ff 72 if (!fIsRead)
73 if(Read(fParticles,fTracks))
74 {
75 Error("GetTrackEvent","Error in reading");
76 return 0x0;
77 }
1b446896 78 return fTracks->GetEvent(n);
79 }
80/********************************************************************/
81
82Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
83 {
84 //returns number of events of particles
64e553ff 85 if (!fIsRead)
86 if(Read(fParticles,fTracks))
87 {
88 Error("GetNumberOfPartEvents","Error in reading");
89 return 0;
90 }
1b446896 91 return fParticles->GetNumberOfEvents();
92 }
93
94/********************************************************************/
95Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
96 {
97 //returns number of events of tracks
64e553ff 98 if (!fIsRead)
99 if(Read(fParticles,fTracks))
100 {
101 Error("GetNumberOfTrackEvents","Error in reading");
102 return 0;
103 }
1b446896 104 return fTracks->GetNumberOfEvents();
105 }
106/********************************************************************/
107
108Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
64e553ff 109{
429f368c 110 cout<<"AliHBTReaderITSv1::Read()"<<endl;
64e553ff 111 Int_t Nevents = 0;
112 AliITSIOTrack *iotrack=new AliITSIOTrack;
113 Int_t currentdir = 0;
114 Int_t Ndirs;
b0d56759 115 Int_t totalNevents = 0;
116
64e553ff 117 if (fDirs)
118 {
119 Ndirs = fDirs->GetEntries();
120 }
121 else
122 {
123 Ndirs = 0;
124 }
125
126 do //do while is good even if
127 {
128 TFile* gAliceFile = OpenGAliceFile(currentdir);
129 if(gAliceFile == 0x0)
1b446896 130 {
64e553ff 131 Error("Read","Can not open the file with gAlice");
132 delete iotrack;
133 return 1;
1b446896 134 }
64e553ff 135 TFile *file = OpenTrackFile(currentdir);
136 if(file == 0x0)
1b446896 137 {
64e553ff 138 Error("Read","Can not open the file with ITS tracks V1");
139 delete iotrack;
140 return 2;
1b446896 141 }
64e553ff 142
143 if (gAlice->TreeE())//check if tree E exists
144 {
145 Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
146 cout<<"________________________________________________________\n";
147 cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
148 cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
149 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
150 }
151 else
152 {//if not return an error
153 Error("Read","Can not find Header tree (TreeE) in gAlice");
154 delete iotrack;
155 return 4;
156 }
b0d56759 157
158 Int_t naccepted = 0;
64e553ff 159 char tname[30];
b0d56759 160
64e553ff 161 for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
162 {
b0d56759 163 cout<<"Reading Event "<<currentEvent;
64e553ff 164
165 sprintf(tname,"TreeT%d",currentEvent);
166 file->cd();
167 TTree *tracktree=(TTree*)file->Get(tname);
168 TBranch *tbranch=tracktree->GetBranch("ITStracks");
169 tbranch->SetAddress(&iotrack);
170
171 gAliceFile->cd();
172 gAlice->GetEvent(currentEvent);
173 gAlice->Particles();
1b446896 174
64e553ff 175 Int_t nentr=(Int_t)tracktree->GetEntries();
b0d56759 176
177 cout<<". Found "<<nentr<<" tracks.";
178 fflush(0);
179
64e553ff 180 for (Int_t i=0; i<nentr; i++)
181 {
1b446896 182
64e553ff 183 tracktree->GetEvent(i);
184 if(!iotrack) continue;
185 Int_t label = iotrack->GetLabel();
186 if (label < 0)
187 {
64e553ff 188 delete iotrack;
6440fe71 189 continue;
64e553ff 190 }
191 TParticle *p = (TParticle*)gAlice->Particle(label);
192 if(!p)
193 {
194 Warning("Read","Can not get particle with label &d",label);
195 continue;
196 }
197 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
198 //if not take next partilce
1b446896 199
64e553ff 200 AliHBTParticle* part = new AliHBTParticle(*p);
201 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
202 //if it does not delete it and take next good track
203
204 Double_t px=iotrack->GetPx();
205 Double_t py=iotrack->GetPy();
206 Double_t pz=iotrack->GetPz();
207 Double_t mass = p->GetMass();
208 Double_t tEtot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);//total energy of the track
209
210 Double_t x= iotrack->GetX();
211 Double_t y= iotrack->GetY();
212 Double_t z= iotrack->GetZ();
213
214 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), px, py , pz, tEtot, x, y, z, 0.);
215 if(Pass(track)) { delete track;continue;}//check if meets all criteria of any of our cuts
216 //if it does not delete it and take next good track
1b446896 217
64e553ff 218 particles->AddParticle(totalNevents,part);//put track and particle on the run
219 tracks->AddParticle(totalNevents,track);
b0d56759 220 naccepted++;
64e553ff 221 }//end loop over tracks in the event
1b446896 222
64e553ff 223 totalNevents++;
b0d56759 224 cout<<" Accepted "<<naccepted<<" tracks"<<endl;
64e553ff 225 }//end of loop over events in current directory
226
227 gAliceFile->Close();
228 delete gAliceFile;
229 gAliceFile = 0;
230
231 file->Close();
232 delete file;
233 file = 0;
429f368c 234 currentdir++;
64e553ff 235 }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
1b446896 236
1b446896 237
64e553ff 238 delete iotrack;
239 fIsRead = kTRUE;
240 return 0;
241
1b446896 242 }
64e553ff 243/********************************************************************/
1b446896 244
64e553ff 245TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
246{
247//opens files to be read for given directoru nomber in fDirs Array
248 const TString& dirname = GetDirName(ndir);
249 if (dirname == "")
250 {
251 Error("OpenGAliceFile","Can not get directory name");
252 return 0x0;
253 }
254 TString filename = dirname + "/" + fITSTracksFileName;
f20bcaaa 255
256 TFile *file = TFile::Open(filename.Data());
64e553ff 257 if (!file)
258 {
259 Error("Read","Can not open file %s",filename.Data());
260 return 0x0;
261 }
262 if (!file->IsOpen())
263 {
264 Error("Read","Can not open file %s",filename.Data());
265 return 0x0;
266 }
267
268 return file;
269}
1b446896 270
271
64e553ff 272/********************************************************************/
273TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
274{
275 const TString& dirname = GetDirName(ndir);
276 if (dirname == "")
277 {
278 Error("OpenGAliceFile","Can not get directory name");
279 return 0x0;
280 }
1b446896 281
64e553ff 282 TString filename = dirname + "/" + fGAliceFileName;
283
284 TFile* gAliceFile = TFile::Open(filename.Data());
285 if ( gAliceFile== 0x0)
1b446896 286 {
64e553ff 287 Error("OpenFiles","Can't open file named %s",filename.Data());
288 return 0x0;
1b446896 289 }
64e553ff 290 if (!gAliceFile->IsOpen())
1b446896 291 {
64e553ff 292 Error("OpenFiles","Can't open file named %s",filename.Data());
293 return 0x0;
1b446896 294 }
1b446896 295
64e553ff 296 if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
297 {
298 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
299 gAliceFile->Close();
300 delete gAliceFile;
301 return 0x0;
302 }
303
304 return gAliceFile;
305}
306
307/********************************************************************/
308/********************************************************************/
309
310