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