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