2 #include "AliHBTReaderITSv2.h"
7 #include <TObjString.h>
10 #include <TParticle.h>
14 #include <AliITStrackV2.h>
15 //#include <AliITSParam.h>
16 #include <AliITStrackerV2.h>
17 #include <AliITSgeom.h>
19 #include "AliHBTRun.h"
20 #include "AliHBTEvent.h"
21 #include "AliHBTParticle.h"
22 #include "AliHBTParticleCut.h"
25 ClassImp(AliHBTReaderITSv2)
28 AliHBTReaderITSv2(const Char_t* trackfilename, const Char_t* clusterfilename,
29 const Char_t* galicefilename)
30 :fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
31 fGAliceFileName(galicefilename)
35 // trackfilename = "AliITStracksV2.root"
36 // clusterfilename = "AliITSclustersV2.root"
37 // galicefilename = "galice.root"
38 fParticles = new AliHBTRun();
39 fTracks = new AliHBTRun();
42 /********************************************************************/
45 AliHBTReaderITSv2(TObjArray* dirs, const Char_t* trackfilename,
46 const Char_t* clusterfilename, const Char_t* galicefilename)
48 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
49 fGAliceFileName(galicefilename)
53 // trackfilename = "AliITStracksV2.root"
54 // clusterfilename = "AliITSclustersV2.root"
55 // galicefilename = "galice.root"
57 fParticles = new AliHBTRun();
58 fTracks = new AliHBTRun();
61 /********************************************************************/
63 AliHBTReaderITSv2::~AliHBTReaderITSv2()
65 if (fParticles) delete fParticles;
66 if (fTracks) delete fTracks;
68 /********************************************************************/
69 /********************************************************************/
71 AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
73 //returns Nth event with simulated particles
75 if(Read(fParticles,fTracks))
77 Error("GetParticleEvent","Error in reading");
81 return fParticles->GetEvent(n);
83 /********************************************************************/
85 AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
87 //returns Nth event with reconstructed tracks
89 if(Read(fParticles,fTracks))
91 Error("GetTrackEvent","Error in reading");
94 return fTracks->GetEvent(n);
96 /********************************************************************/
98 Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
100 //returns number of events of particles
102 if(Read(fParticles,fTracks))
104 Error("GetNumberOfPartEvents","Error in reading");
107 return fParticles->GetNumberOfEvents();
110 /********************************************************************/
111 Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
113 //returns number of events of tracks
115 if(Read(fParticles,fTracks))
117 Error("GetNumberOfTrackEvents","Error in reading");
120 return fTracks->GetNumberOfEvents();
124 /********************************************************************/
125 /********************************************************************/
126 Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
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
135 TFile *aTracksFile;//file with tracks
136 TFile *aClustersFile;//file with clusters
137 TFile *aGAliceFile;//file name with galice
139 AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
140 TTree *tracktree; // tree for tracks
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
147 char tname[100]; //buffer for tree name
148 AliITStrackV2 *iotrack= new AliITStrackV2(); //buffer track for reading data from tree
150 if (!particles) //check if an object is instatiated
152 Error("Read"," particles object must instatiated before passing it to the reader");
154 if (!tracks) //check if an object is instatiated
156 Error("Read"," tracks object must instatiated before passing it to the reader");
158 particles->Reset();//clear runs == delete all old events
161 if (fDirs) //if array with directories is supplied by user
163 Ndirs = fDirs->GetEntries(); //get the number if directories
167 Ndirs = 0; //if the array is not supplied read only from current directory
170 // cout<<"Found "<<Ndirs<<" directory entries"<<endl;
172 do //do while is good even if Ndirs==0 (than read from current directory)
174 if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
176 Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
181 if (gAlice->TreeE())//check if tree E exists
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());
190 {//if not return an error
191 Error("Read","Can not find Header tree (TreeE) in gAlice");
196 AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom");
199 Error("Read","Can't get the ITS geometry!");
204 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
206 cout<<"Reading Event "<<currentEvent<<endl;
209 gAlice->GetEvent(currentEvent);
210 TParticle * part = gAlice->Particle(0);
211 Double_t orz=part->Vz();
214 tracker = new AliITStrackerV2(geom,currentEvent,orz);
216 sprintf(tname,"TreeT_ITS_%d",currentEvent);
218 tracktree=(TTree*)aTracksFile->Get(tname);
221 Error("Read","Can't get a tree with ITS tracks");
226 TBranch *tbranch=tracktree->GetBranch("tracks");
228 Ntracks=(Int_t)tracktree->GetEntries();
233 for (i=0; i<Ntracks; i++) //loop over all tpc tracks
235 if(i%100 == 0)cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<"\r";
237 tbranch->SetAddress(&iotrack);
238 tracktree->GetEvent(i);
240 label=iotrack->GetLabel();
246 tracker->CookLabel(iotrack,0.);
247 Int_t itsLabel=iotrack->GetLabel();
248 if (itsLabel != label)
254 TParticle *p = (TParticle*)gAlice->Particle(label);
255 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
256 //if not take next partilce
258 AliHBTParticle* part = new AliHBTParticle(*p);
259 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
260 //if it does not delete it and take next good track
263 iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
264 iotrack->PropagateToVertex();
266 iotrack->GetExternalParameters(xk,par); //get properties of the track
267 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
268 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
269 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
271 pt=1.0/TMath::Abs(par[4]);
273 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
274 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
275 Double_t tpz = pt * lam; //track z coordinate of momentum
277 Double_t mass = p->GetMass();
278 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
280 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
281 if(Pass(track))//check if meets all criteria of any of our cuts
282 //if it does not delete it and take next good track
288 particles->AddParticle(totalNevents,part);//put track and particle on the run
289 tracks->AddParticle(totalNevents,track);
291 }//end of loop over tracks in the event
293 aTracksFile->Delete(tname);
294 aTracksFile->Delete("tracks");
298 CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
299 cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<endl;
301 }//end of loop over events in current directory
303 }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
310 /********************************************************************/
311 Int_t AliHBTReaderITSv2::OpenFiles
312 (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
314 //opens all the files
317 const TString& dirname = GetDirName(event);
320 Error("OpenFiles","Can not get directory name");
324 TString filename = dirname +"/"+ fTrackFileName;
325 aTracksFile = TFile::Open(filename.Data());
326 if ( aTracksFile == 0x0 )
328 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
331 if (!aTracksFile->IsOpen())
333 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
337 filename = dirname +"/"+ fClusterFileName;
338 aClustersFile = TFile::Open(filename.Data());
339 if ( aClustersFile == 0x0 )
341 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
344 if (!aClustersFile->IsOpen())
346 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
350 filename = dirname +"/"+ fGAliceFileName;
351 agAliceFile = TFile::Open(filename.Data());
352 if ( agAliceFile== 0x0)
354 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
357 if (!agAliceFile->IsOpen())
359 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
363 if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
365 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
371 /********************************************************************/
373 /********************************************************************/
375 void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
382 clustersFile->Close();
397 /********************************************************************/