1 #include "AliHBTReaderTPC.h"
11 #include <AliTPCtrack.h>
12 #include <AliTPCParam.h>
13 #include <AliTPCtracker.h>
15 #include "AliHBTRun.h"
16 #include "AliHBTEvent.h"
17 #include "AliHBTParticle.h"
18 #include "AliHBTParticleCut.h"
21 ClassImp(AliHBTReaderTPC)
22 //reader for TPC tracking
23 //needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc
25 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
26 //Piotr.Skowronski@cern.ch
29 AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
30 const Char_t* galicefilename):
31 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
32 fGAliceFileName(galicefilename)
36 // trackfilename = "AliTPCtracks.root"
37 // clusterfilename = "AliTPCclusters.root"
38 // galicefilename = "" - this means: Do not open gAlice file -
39 // just leave the global pointer untached
41 fParticles = new AliHBTRun();
42 fTracks = new AliHBTRun();
45 /********************************************************************/
47 AliHBTReaderTPC(TObjArray* dirs,
48 const Char_t* trackfilename, const Char_t* clusterfilename,
49 const Char_t* galicefilename):
50 AliHBTReader(dirs), fTrackFileName(trackfilename),
51 fClusterFileName(clusterfilename),fGAliceFileName(galicefilename)
56 // trackfilename = "AliTPCtracks.root"
57 // clusterfilename = "AliTPCclusters.root"
58 // galicefilename = "" - this means: Do not open gAlice file -
59 // just leave the global pointer untached
61 fParticles = new AliHBTRun();
62 fTracks = new AliHBTRun();
66 /********************************************************************/
68 AliHBTReaderTPC::~AliHBTReaderTPC()
74 /********************************************************************/
76 AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
78 //returns Nth event with simulated particles
80 if(Read(fParticles,fTracks))
82 Error("GetParticleEvent","Error in reading");
85 return fParticles->GetEvent(n);
87 /********************************************************************/
88 AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
90 //returns Nth event with reconstructed tracks
92 if(Read(fParticles,fTracks))
94 Error("GetTrackEvent","Error in reading");
97 return fTracks->GetEvent(n);
99 /********************************************************************/
101 Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
103 //returns number of events of particles
105 if ( Read(fParticles,fTracks))
107 Error("GetNumberOfPartEvents","Error in reading");
110 return fParticles->GetNumberOfEvents();
113 /********************************************************************/
114 Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
116 //returns number of events of tracks
118 if(Read(fParticles,fTracks))
120 Error("GetNumberOfTrackEvents","Error in reading");
123 return fTracks->GetNumberOfEvents();
125 /********************************************************************/
128 Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
130 //reads data and puts put to the particles and tracks objects
131 //reurns 0 if everything is OK
133 cout<<"AliHBTReaderTPC::Read()"<<endl;
134 Int_t i; //iterator and some temprary values
136 Int_t totalNevents = 0;
137 TFile *aTracksFile;//file with tracks
138 TFile *aClustersFile;//file with clusters
139 TFile *aGAliceFile;//!ile name with galice
141 if (!particles) //check if an object is instatiated
143 Error("Read"," particles object must instatiated before passing it to the reader");
145 if (!tracks) //check if an object is instatiated
147 Error("Read"," tracks object must instatiated before passing it to the reader");
149 particles->Reset();//clear runs == delete all old events
152 TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
153 tarray->SetOwner(); //set the ownership of the objects it contains
154 //when array is is deleted or cleared all objects
155 //that it contains are deleted
156 Int_t currentdir = 0;
159 if (fDirs) //if array with directories is supplied by user
161 Ndirs = fDirs->GetEntries(); //get the number if directories
165 Ndirs = 0; //if the array is not supplied read only from current directory
168 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
171 if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
173 Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
178 if (gAlice->TreeE())//check if tree E exists
180 Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
181 cout<<"________________________________________________________\n";
182 cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
183 cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
184 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
187 {//if not return an error
188 Error("Read","Can not find Header tree (TreeE) in gAlice");
192 aClustersFile->cd();//set cluster file active
193 AliTPCParam *TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60");
196 Error("Read","TPC parameters have not been found !\n");
201 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
203 cout<<"Reading Event "<<currentEvent<<endl;
204 /**************************************/
205 /**************************************/
206 /**************************************/
208 aTracksFile->cd();//set track file active
210 Char_t treename[100];
211 sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
215 tracktree=(TTree*)aTracksFile->Get(treename);//get the tree
216 if (!tracktree) //check if we got the tree
217 {//if not return with error
218 Error("Read","Can't get a tree with TPC tracks !\n");
223 TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
224 if (!trackbranch) ////check if we got the branch
225 {//if not return with error
226 Error("Read","Can't get a branch with TPC tracks !\n");
229 Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
230 cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
231 //Copy tracks to array
233 AliTPCtrack *iotrack=0;
235 aClustersFile->cd();//set cluster file active
236 AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
237 if (!tracker) //check if it has created succeffuly
238 {//if not return with error
239 Error("Read","Can't get a tracker !\n");
242 tracker->LoadInnerSectors();
243 tracker->LoadOuterSectors();
245 for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
247 iotrack=new AliTPCtrack; //create new tracks
248 trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
249 tracktree->GetEvent(i); //stream track i to the iotrack
250 tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
251 //which is the label of corresponding simulated particle
252 tarray->AddLast(iotrack); //put the track in the array
255 aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
256 aTracksFile->Delete("tracks");//delete branch from memmory
257 delete tracker; //delete tracker
265 Float_t phi, lam, pt;//angles and transverse momentum
266 Int_t label; //label of the current track
269 gAlice->GetEvent(currentEvent);
273 for (i=0; i<NTPCtracks; i++) //loop over all good tracks
275 iotrack = (AliTPCtrack*)tarray->At(i);
276 label = iotrack->GetLabel();
278 if (label < 0) continue;
280 TParticle *p = (TParticle*)gAlice->Particle(label);
282 if(p == 0x0) continue; //if returned pointer is NULL
283 if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
285 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
286 //if not take next partilce
288 AliHBTParticle* part = new AliHBTParticle(*p);
289 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
290 //if it does not delete it and take next good track
293 iotrack->PropagateToVertex();
295 iotrack->GetExternalParameters(xk,par); //get properties of the track
296 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
297 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
298 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
300 pt=1.0/TMath::Abs(par[4]);
302 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
303 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
304 Double_t tpz = pt * lam; //track z coordinate of momentum
306 Double_t mass = p->GetMass();
307 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
309 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
310 if(Pass(track))//check if meets all criteria of any of our cuts
311 //if it does not delete it and take next good track
317 particles->AddParticle(totalNevents,part);//put track and particle on the run
318 tracks->AddParticle(totalNevents,track);
321 tarray->Clear(); //clear the array
323 /**************************************/
324 /**************************************/
325 /**************************************/
329 //save environment (resouces) --
330 //clean your place after the work
331 CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
333 }while(currentdir < Ndirs);
341 /********************************************************************/
342 Int_t AliHBTReaderTPC::OpenFiles
343 (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
345 //opens all the files
348 const TString& dirname = GetDirName(event);
351 Error("OpenFiles","Can not get directory name");
355 TString filename = dirname +"/"+ fTrackFileName;
356 aTracksFile = TFile::Open(filename.Data());
357 if ( aTracksFile == 0x0 )
359 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
362 if (!aTracksFile->IsOpen())
364 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
368 filename = dirname +"/"+ fClusterFileName;
369 aClustersFile = TFile::Open(filename.Data());
370 if ( aClustersFile == 0x0 )
372 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
375 if (!aClustersFile->IsOpen())
377 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
381 filename = dirname +"/"+ fGAliceFileName;
382 agAliceFile = TFile::Open(filename.Data());
383 if ( agAliceFile== 0x0)
385 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
388 if (!agAliceFile->IsOpen())
390 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
394 if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
396 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
402 /********************************************************************/
404 /********************************************************************/
406 void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
412 clustersFile->Close();
427 /********************************************************************/
429 /********************************************************************/
430 /********************************************************************/
431 /********************************************************************/