Adaption to new fluka common blocks (E. Futo)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
CommitLineData
1b446896 1#include "AliHBTReaderTPC.h"
8fe7288a 2
1b446896 3#include <TTree.h>
4#include <TFile.h>
16701d1b 5#include <TParticle.h>
1b446896 6
16701d1b 7#include <AliRun.h>
8#include <AliMagF.h>
1b446896 9#include <AliTPCtrack.h>
10#include <AliTPCParam.h>
1b446896 11
1b446896 12#include "AliHBTRun.h"
13#include "AliHBTEvent.h"
14#include "AliHBTParticle.h"
15#include "AliHBTParticleCut.h"
16
17
18ClassImp(AliHBTReaderTPC)
8fe7288a 19//______________________________________________
20//
21// class AliHBTReaderTPC
22//
1b446896 23//reader for TPC tracking
98295f4b 24//needs galice.root, AliTPCtracks.root
1b446896 25//
26//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
27//Piotr.Skowronski@cern.ch
8fe7288a 28//
1b446896 29
98295f4b 30
8fe7288a 31AliHBTReaderTPC:: AliHBTReaderTPC(const Char_t* trackfilename,
32 const Char_t* clusterfilename,
33 const Char_t* galicefilename):
98295f4b 34 fParticles(new AliHBTRun()),
35 fTracks(new AliHBTRun()),
8fe7288a 36 fTrackFileName(trackfilename),
37 fClusterFileName(clusterfilename),
98295f4b 38 fGAliceFileName(galicefilename),
39 fIsRead(kFALSE),
40 fMagneticField(0.4),
41 fUseMagFFromRun(kTRUE)
1b446896 42{
16701d1b 43 //constructor,
1b446896 44 //Defaults:
45 // trackfilename = "AliTPCtracks.root"
46 // clusterfilename = "AliTPCclusters.root"
1b446896 47 // galicefilename = "" - this means: Do not open gAlice file -
48 // just leave the global pointer untached
98295f4b 49 //this class is not supposed to be written, that is why I let myself and user
50 //for comfort of having default constructor which allocates dynamic memmory
51 //in case it is going to be written, it should be changed
1b446896 52
16701d1b 53}
54/********************************************************************/
8fe7288a 55AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs,
df9a1768 56 const Char_t* trackfilename, const Char_t* clusterfilename,
57 const Char_t* galicefilename):
8fe7288a 58 AliHBTReader(dirs),
98295f4b 59 fParticles(new AliHBTRun()),
60 fTracks(new AliHBTRun()),
8fe7288a 61 fTrackFileName(trackfilename),
62 fClusterFileName(clusterfilename),
98295f4b 63 fGAliceFileName(galicefilename),
64 fIsRead(kFALSE),
65 fMagneticField(0.4),
66 fUseMagFFromRun(kTRUE)
16701d1b 67{
68 //constructor,
69 //Defaults:
70 // trackfilename = "AliTPCtracks.root"
71 // clusterfilename = "AliTPCclusters.root"
72 // galicefilename = "" - this means: Do not open gAlice file -
73 // just leave the global pointer untached
1b446896 74}
75/********************************************************************/
76
77AliHBTReaderTPC::~AliHBTReaderTPC()
78 {
79 //desctructor
80 delete fParticles;
81 delete fTracks;
82 }
83/********************************************************************/
84
85AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
86 {
87 //returns Nth event with simulated particles
16701d1b 88 if (!fIsRead)
89 if(Read(fParticles,fTracks))
90 {
91 Error("GetParticleEvent","Error in reading");
92 return 0x0;
93 }
1b446896 94 return fParticles->GetEvent(n);
95 }
96/********************************************************************/
8fe7288a 97
1b446896 98AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
99 {
100 //returns Nth event with reconstructed tracks
16701d1b 101 if (!fIsRead)
102 if(Read(fParticles,fTracks))
103 {
104 Error("GetTrackEvent","Error in reading");
105 return 0x0;
106 }
1b446896 107 return fTracks->GetEvent(n);
108 }
109/********************************************************************/
110
111Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
112 {
113 //returns number of events of particles
16701d1b 114 if (!fIsRead)
115 if ( Read(fParticles,fTracks))
116 {
117 Error("GetNumberOfPartEvents","Error in reading");
118 return 0;
119 }
1b446896 120 return fParticles->GetNumberOfEvents();
121 }
122
123/********************************************************************/
124Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
125 {
126 //returns number of events of tracks
16701d1b 127 if (!fIsRead)
128 if(Read(fParticles,fTracks))
129 {
130 Error("GetNumberOfTrackEvents","Error in reading");
131 return 0;
132 }
1b446896 133 return fTracks->GetNumberOfEvents();
134 }
135/********************************************************************/
136
137
138Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
139 {
140 //reads data and puts put to the particles and tracks objects
141 //reurns 0 if everything is OK
142 //
8fe7288a 143 Info("Read","");
1b446896 144 Int_t i; //iterator and some temprary values
16701d1b 145 Int_t Nevents = 0;
146 Int_t totalNevents = 0;
147 TFile *aTracksFile;//file with tracks
148 TFile *aClustersFile;//file with clusters
149 TFile *aGAliceFile;//!ile name with galice
150
1b446896 151 if (!particles) //check if an object is instatiated
152 {
7da9bdd2 153 Error("Read"," particles object must instatiated before passing it to the reader");
1b446896 154 }
155 if (!tracks) //check if an object is instatiated
156 {
7da9bdd2 157 Error("Read"," tracks object must instatiated before passing it to the reader");
1b446896 158 }
159 particles->Reset();//clear runs == delete all old events
160 tracks->Reset();
16701d1b 161
162 TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
163 tarray->SetOwner(); //set the ownership of the objects it contains
164 //when array is is deleted or cleared all objects
165 //that it contains are deleted
166 Int_t currentdir = 0;
b71a5edb 167
168 Int_t Ndirs;
169 if (fDirs) //if array with directories is supplied by user
170 {
171 Ndirs = fDirs->GetEntries(); //get the number if directories
172 }
173 else
174 {
175 Ndirs = 0; //if the array is not supplied read only from current directory
176 }
177
16701d1b 178 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
1b446896 179 {
16701d1b 180
181 if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
182 {
7da9bdd2 183 Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
5ee8b891 184 currentdir++;
185 continue;
16701d1b 186 }
1b446896 187
188
16701d1b 189 if (gAlice->TreeE())//check if tree E exists
190 {
191 Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
8fe7288a 192 Info("Read","________________________________________________________");
193 Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
98295f4b 194 Float_t mf;
195 if (fUseMagFFromRun)
196 {
197 mf = gAlice->Field()->SolenoidField();
198 Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
199 }
200 else
201 {
202 Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
203 mf = fMagneticField*10.;
204 }
205 AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
16701d1b 206 }
207 else
208 {//if not return an error
7da9bdd2 209 Error("Read","Can not find Header tree (TreeE) in gAlice");
5ee8b891 210 currentdir++;
211 continue;
16701d1b 212 }
1b446896 213
98295f4b 214 aGAliceFile->cd();//set cluster file active
215 AliTPCParam *TPCParam= (AliTPCParam*)aGAliceFile->Get("75x40_100x60_150x60");
16701d1b 216 if (!TPCParam)
217 {
98295f4b 218 TPCParam= (AliTPCParam*)aGAliceFile->Get("75x40_100x60");
2f1287d2 219 if (!TPCParam)
220 {
221 Error("Read","TPC parameters have not been found !\n");
222 currentdir++;
223 continue;
224 }
16701d1b 225 }
1b446896 226
1b446896 227
16701d1b 228 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
229 {
8fe7288a 230 Info("Read","Reading Event %d",currentEvent);
16701d1b 231 /**************************************/
232 /**************************************/
233 /**************************************/
1b446896 234
16701d1b 235 aTracksFile->cd();//set track file active
236
1b446896 237 Char_t treename[100];
238 sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
239
240 TTree *tracktree=0;
241
16701d1b 242 tracktree=(TTree*)aTracksFile->Get(treename);//get the tree
1b446896 243 if (!tracktree) //check if we got the tree
244 {//if not return with error
7da9bdd2 245 Error("Read","Can't get a tree with TPC tracks !\n");
5ee8b891 246 continue;
1b446896 247 }
248
249 TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
250 if (!trackbranch) ////check if we got the branch
251 {//if not return with error
7da9bdd2 252 Error("Read","Can't get a branch with TPC tracks !\n");
5ee8b891 253 continue;
1b446896 254 }
255 Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
8fe7288a 256 Info("Read","Found %d TPC tracks.",NTPCtracks);
1b446896 257 //Copy tracks to array
258
259 AliTPCtrack *iotrack=0;
260
98295f4b 261// aClustersFile->cd();//set cluster file active
1b446896 262
263 for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
264 {
265 iotrack=new AliTPCtrack; //create new tracks
266 trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
267 tracktree->GetEvent(i); //stream track i to the iotrack
1b446896 268 tarray->AddLast(iotrack); //put the track in the array
269 }
270
16701d1b 271 aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
272 aTracksFile->Delete("tracks");//delete branch from memmory
1b446896 273
1b446896 274 trackbranch = 0x0;
275 tracktree = 0x0;
1b446896 276
277 Double_t xk;
278 Double_t par[5];
279 Float_t phi, lam, pt;//angles and transverse momentum
280 Int_t label; //label of the current track
16701d1b 281
282 aGAliceFile->cd();
283 gAlice->GetEvent(currentEvent);
284
285 gAlice->Particles();
286
287 for (i=0; i<NTPCtracks; i++) //loop over all good tracks
1b446896 288 {
16701d1b 289 iotrack = (AliTPCtrack*)tarray->At(i);
290 label = iotrack->GetLabel();
8fe7288a 291
16701d1b 292 if (label < 0) continue;
1b446896 293
8fe7288a 294
16701d1b 295 TParticle *p = (TParticle*)gAlice->Particle(label);
1b446896 296
5ee8b891 297 if(p == 0x0) continue; //if returned pointer is NULL
298 if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
f1e2da22 299
16701d1b 300 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
301 //if not take next partilce
8089a083 302 AliHBTParticle* part = new AliHBTParticle(*p,i);
1b446896 303 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
304 //if it does not delete it and take next good track
305
16701d1b 306 iotrack->PropagateToVertex();
7da9bdd2 307
1b446896 308 iotrack->GetExternalParameters(xk,par); //get properties of the track
309 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
310 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
311 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
312 lam=par[3];
313 pt=1.0/TMath::Abs(par[4]);
314
315 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
316 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
317 Double_t tpz = pt * lam; //track z coordinate of momentum
318
16701d1b 319 Double_t mass = p->GetMass();
1b446896 320 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
321
8089a083 322 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
5ee8b891 323 if(Pass(track))//check if meets all criteria of any of our cuts
7da9bdd2 324 //if it does not delete it and take next good track
325 {
326 delete track;
327 delete part;
328 continue;
329 }
16701d1b 330 particles->AddParticle(totalNevents,part);//put track and particle on the run
331 tracks->AddParticle(totalNevents,track);
1b446896 332 }
333 tarray->Clear(); //clear the array
16701d1b 334
335 /**************************************/
1b446896 336 /**************************************/
16701d1b 337 /**************************************/
338 totalNevents++;
339 }
16701d1b 340 //save environment (resouces) --
341 //clean your place after the work
342 CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
343 currentdir++;
b71a5edb 344 }while(currentdir < Ndirs);
16701d1b 345
346
1b446896 347 delete tarray;
1b446896 348 fIsRead = kTRUE;
349 return 0;
350 }
351
352/********************************************************************/
16701d1b 353Int_t AliHBTReaderTPC::OpenFiles
354(TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
1b446896 355{
356 //opens all the files
16701d1b 357
358
359 const TString& dirname = GetDirName(event);
360 if (dirname == "")
361 {
362 Error("OpenFiles","Can not get directory name");
363 return 4;
364 }
365
366 TString filename = dirname +"/"+ fTrackFileName;
367 aTracksFile = TFile::Open(filename.Data());
368 if ( aTracksFile == 0x0 )
1b446896 369 {
16701d1b 370 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
1b446896 371 return 1;
372 }
16701d1b 373 if (!aTracksFile->IsOpen())
374 {
375 Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
376 return 1;
377 }
98295f4b 378
379/*
16701d1b 380 filename = dirname +"/"+ fClusterFileName;
381 aClustersFile = TFile::Open(filename.Data());
382 if ( aClustersFile == 0x0 )
1b446896 383 {
16701d1b 384 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
1b446896 385 return 2;
386 }
16701d1b 387 if (!aClustersFile->IsOpen())
388 {
389 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
390 return 2;
391 }
98295f4b 392*/
16701d1b 393 filename = dirname +"/"+ fGAliceFileName;
394 agAliceFile = TFile::Open(filename.Data());
395 if ( agAliceFile== 0x0)
396 {
397 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
398 return 3;
399 }
400 if (!agAliceFile->IsOpen())
401 {
402 Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
403 return 3;
404 }
405
406 if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
407 {
408 Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
409 return 5;
410 }
1b446896 411
16701d1b 412 return 0;
1b446896 413}
16701d1b 414/********************************************************************/
1b446896 415
416/********************************************************************/
417
16701d1b 418void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
1b446896 419{
420 //closes the files
16701d1b 421 tracksFile->Close();
7da9bdd2 422 delete tracksFile;
16701d1b 423 tracksFile = 0x0;
98295f4b 424
425// clustersFile->Close();
426// delete clustersFile;
427// clustersFile = 0x0;
7da9bdd2 428
429 delete gAlice;
430 gAlice = 0;
431
432 if (gAliceFile)
433 {
434 gAliceFile->Close();
435 delete gAliceFile;
436 gAliceFile = 0x0;
437 }
1b446896 438}
439
440/********************************************************************/
441
442/********************************************************************/
443/********************************************************************/
444/********************************************************************/
445