2 #include "AliHBTReaderITSv1.h"
3 #include "AliHBTEvent.h"
5 #include "AliHBTParticle.h"
6 #include "AliHBTParticleCut.h"
14 #include <TObjArray.h>
15 #include <TParticle.h>
17 #include <TObjString.h>
21 #include <AliKalmanTrack.h>
22 #include <AliITSIOTrack.h>
24 ClassImp(AliHBTReaderITSv1)
25 /********************************************************************/
28 AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
29 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
31 fParticles = new AliHBTRun();
32 fTracks = new AliHBTRun();
35 /********************************************************************/
38 AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
40 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
42 fParticles = new AliHBTRun();
43 fTracks = new AliHBTRun();
46 /********************************************************************/
48 AliHBTReaderITSv1::~AliHBTReaderITSv1()
53 /********************************************************************/
55 AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
57 //returns Nth event with simulated particles
59 if(Read(fParticles,fTracks))
61 Error("GetParticleEvent","Error in reading");
65 return fParticles->GetEvent(n);
67 /********************************************************************/
69 AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
71 //returns Nth event with reconstructed tracks
73 if(Read(fParticles,fTracks))
75 Error("GetTrackEvent","Error in reading");
78 return fTracks->GetEvent(n);
80 /********************************************************************/
82 Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
84 //returns number of events of particles
86 if(Read(fParticles,fTracks))
88 Error("GetNumberOfPartEvents","Error in reading");
91 return fParticles->GetNumberOfEvents();
94 /********************************************************************/
95 Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
97 //returns number of events of tracks
99 if(Read(fParticles,fTracks))
101 Error("GetNumberOfTrackEvents","Error in reading");
104 return fTracks->GetNumberOfEvents();
106 /********************************************************************/
108 Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
110 cout<<"AliHBTReaderITSv1::Read()"<<endl;
112 AliITSIOTrack *iotrack=new AliITSIOTrack;
113 Int_t currentdir = 0;
115 Int_t totalNevents = 0;
119 Ndirs = fDirs->GetEntries();
126 do //do while is good even if
128 TFile* gAliceFile = OpenGAliceFile(currentdir);
129 if(gAliceFile == 0x0)
131 Error("Read","Can not open the file with gAlice");
135 if (gAlice->TreeE())//check if tree E exists
137 Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
138 cout<<"________________________________________________________\n";
139 cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
140 cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
141 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
144 {//if not return an error
145 Error("Read","Can not find Header tree (TreeE) in gAlice");
150 TFile *file = OpenTrackFile(currentdir);
153 Error("Read","Can not open the file with ITS tracks V1");
161 for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
163 cout<<"Reading Event "<<currentEvent;
165 sprintf(tname,"TreeT%d",currentEvent);
167 TTree *tracktree=(TTree*)file->Get(tname);
168 TBranch *tbranch=tracktree->GetBranch("ITStracks");
169 tbranch->SetAddress(&iotrack);
172 gAlice->GetEvent(currentEvent);
175 Int_t nentr=(Int_t)tracktree->GetEntries();
177 cout<<". Found "<<nentr<<" tracks.";
180 for (Int_t i=0; i<nentr; i++)
183 tracktree->GetEvent(i);
184 if(!iotrack) continue;
185 Int_t label = iotrack->GetLabel();
191 TParticle *p = (TParticle*)gAlice->Particle(label);
194 Warning("Read","Can not get particle with label &d",label);
197 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
198 //if not take next partilce
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
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
210 Double_t x= iotrack->GetX();
211 Double_t y= iotrack->GetY();
212 Double_t z= iotrack->GetZ();
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
218 particles->AddParticle(totalNevents,part);//put track and particle on the run
219 tracks->AddParticle(totalNevents,track);
221 }//end loop over tracks in the event
224 cout<<" Accepted "<<naccepted<<" tracks"<<endl;
225 }//end of loop over events in current directory
235 }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
243 /********************************************************************/
245 TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
247 //opens files to be read for given directoru nomber in fDirs Array
248 const TString& dirname = GetDirName(ndir);
251 Error("OpenGAliceFile","Can not get directory name");
254 TString filename = dirname + "/" + fITSTracksFileName;
256 TFile *file = TFile::Open(filename.Data());
259 Error("Read","Can not open file %s",filename.Data());
264 Error("Read","Can not open file %s",filename.Data());
272 /********************************************************************/
273 TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
275 const TString& dirname = GetDirName(ndir);
278 Error("OpenGAliceFile","Can not get directory name");
282 TString filename = dirname + "/" + fGAliceFileName;
284 TFile* gAliceFile = TFile::Open(filename.Data());
285 if ( gAliceFile== 0x0)
287 Error("OpenFiles","Can't open file named %s",filename.Data());
290 if (!gAliceFile->IsOpen())
292 Error("OpenFiles","Can't open file named %s",filename.Data());
296 if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
298 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
307 /********************************************************************/
308 /********************************************************************/