Coding Violations Corrected.
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
CommitLineData
1b446896 1#include "AliHBTReaderTPC.h"
3f745d47 2//______________________________________________
3//
4// class AliHBTReaderTPC
5//
6// reader for TPC tracking
7// needs galice.root, AliTPCtracks.root, AliTPCclusters.root
8//
9// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
10// Piotr.Skowronski@cern.ch
11//
12///////////////////////////////////////////////////////////////////////////
1b446896 13#include <TTree.h>
14#include <TFile.h>
16701d1b 15#include <TParticle.h>
bfb09ece 16#include <TH1.h>
1b446896 17
bed069a4 18#include <Varargs.h>
19
16701d1b 20#include <AliRun.h>
88cb7938 21#include <AliLoader.h>
22#include <AliStack.h>
16701d1b 23#include <AliMagF.h>
1b446896 24#include <AliTPCtrack.h>
25#include <AliTPCParam.h>
88cb7938 26#include <AliTPCtracker.h>
bed069a4 27#include <AliTPCLoader.h>
1b446896 28
1b446896 29#include "AliHBTRun.h"
30#include "AliHBTEvent.h"
31#include "AliHBTParticle.h"
32#include "AliHBTParticleCut.h"
f5de2f09 33#include "AliHBTTrackPoints.h"
66d1d1a4 34#include "AliHBTClusterMap.h"
bfb09ece 35
36
1b446896 37ClassImp(AliHBTReaderTPC)
3f745d47 38
bed069a4 39AliHBTReaderTPC::AliHBTReaderTPC():
40 fFileName("galice.root"),
41 fRunLoader(0x0),
42 fTPCLoader(0x0),
43 fMagneticField(0.0),
3f745d47 44 fUseMagFFromRun(kTRUE),
f5de2f09 45 fNTrackPoints(0),
46 fdR(0.0),
66d1d1a4 47 fClusterMap(kFALSE),
3f745d47 48 fNClustMin(0),
f5de2f09 49 fNClustMax(150),
3f745d47 50 fNChi2PerClustMin(0.0),
51 fNChi2PerClustMax(10e5),
f5de2f09 52 fC00Min(0.0),
53 fC00Max(10e5),
54 fC11Min(0.0),
55 fC11Max(10e5),
56 fC22Min(0.0),
57 fC22Max(10e5),
58 fC33Min(0.0),
59 fC33Max(10e5),
3f745d47 60 fC44Min(0.0),
61 fC44Max(10e5)
88cb7938 62{
f5de2f09 63 //constructor
88cb7938 64
88cb7938 65}
f5de2f09 66/********************************************************************/
98295f4b 67
bed069a4 68AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
69 fFileName(galicefilename),
70 fRunLoader(0x0),
71 fTPCLoader(0x0),
72 fMagneticField(0.0),
3f745d47 73 fUseMagFFromRun(kTRUE),
f5de2f09 74 fNTrackPoints(0),
75 fdR(0.0),
3f745d47 76 fNClustMin(0),
f5de2f09 77 fNClustMax(150),
3f745d47 78 fNChi2PerClustMin(0.0),
79 fNChi2PerClustMax(10e5),
f5de2f09 80 fC00Min(0.0),
81 fC00Max(10e5),
82 fC11Min(0.0),
83 fC11Max(10e5),
84 fC22Min(0.0),
85 fC22Max(10e5),
86 fC33Min(0.0),
87 fC33Max(10e5),
3f745d47 88 fC44Min(0.0),
89 fC44Max(10e5)
1b446896 90{
16701d1b 91 //constructor,
1b446896 92 //Defaults:
16701d1b 93}
94/********************************************************************/
f5de2f09 95
88cb7938 96AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
bed069a4 97 AliHBTReader(dirs),
98 fFileName(galicefilename),
99 fRunLoader(0x0),
100 fTPCLoader(0x0),
101 fMagneticField(0.0),
3f745d47 102 fUseMagFFromRun(kTRUE),
f5de2f09 103 fNTrackPoints(0),
104 fdR(0.0),
66d1d1a4 105 fClusterMap(kFALSE),
3f745d47 106 fNClustMin(0),
f5de2f09 107 fNClustMax(150),
3f745d47 108 fNChi2PerClustMin(0.0),
109 fNChi2PerClustMax(10e5),
f5de2f09 110 fC00Min(0.0),
111 fC00Max(10e5),
112 fC11Min(0.0),
113 fC11Max(10e5),
114 fC22Min(0.0),
115 fC22Max(10e5),
116 fC33Min(0.0),
117 fC33Max(10e5),
3f745d47 118 fC44Min(0.0),
119 fC44Max(10e5)
16701d1b 120{
121 //constructor,
122 //Defaults:
16701d1b 123 // galicefilename = "" - this means: Do not open gAlice file -
124 // just leave the global pointer untached
88cb7938 125
1b446896 126}
127/********************************************************************/
128
129AliHBTReaderTPC::~AliHBTReaderTPC()
bed069a4 130{
1b446896 131 //desctructor
3f745d47 132 if (AliHBTParticle::GetDebug())
133 {
134 Info("~AliHBTReaderTPC","deleting Run Loader");
135 AliLoader::SetDebug(AliHBTParticle::GetDebug());
136 }
137
bed069a4 138 delete fRunLoader;
3f745d47 139
140 if (AliHBTParticle::GetDebug())
141 {
142 Info("~AliHBTReaderTPC","deleting Run Loader Done");
143 }
bed069a4 144}
1b446896 145/********************************************************************/
bed069a4 146
147void AliHBTReaderTPC::Rewind()
148{
149 delete fRunLoader;
150 fRunLoader = 0x0;
151 fCurrentDir = 0;
152 fNEventsRead= 0;
bfb09ece 153 if (fTrackCounter) fTrackCounter->Reset();
bed069a4 154}
1b446896 155/********************************************************************/
156
bed069a4 157Int_t AliHBTReaderTPC::ReadNext()
1b446896 158 {
159 //reads data and puts put to the particles and tracks objects
160 //reurns 0 if everything is OK
161 //
8fe7288a 162 Info("Read","");
bfb09ece 163
16701d1b 164 TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
165 tarray->SetOwner(); //set the ownership of the objects it contains
166 //when array is is deleted or cleared all objects
167 //that it contains are deleted
b71a5edb 168
bed069a4 169 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
170 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
171
172 fParticlesEvent->Reset();
173 fTracksEvent->Reset();
b71a5edb 174
16701d1b 175 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
1b446896 176 {
bed069a4 177
178 if (fRunLoader == 0x0)
179 if (OpenNextSession()) continue;//directory counter is increased inside in case of error
180
181 if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
88cb7938 182 {
bed069a4 183 //read next directory
184 delete fRunLoader;//close current session
185 fRunLoader = 0x0;//assure pointer is null
186 fCurrentDir++;//go to next dir
187 continue;//directory counter is increased inside in case of error
16701d1b 188 }
bed069a4 189
190 Info("ReadNext","Reading Event %d",fCurrentEvent);
1b446896 191
bed069a4 192 fRunLoader->GetEvent(fCurrentEvent);
193 TTree *tracktree = fTPCLoader->TreeT();//get the tree
194 if (!tracktree) //check if we got the tree
195 {//if not return with error
196 Error("ReadNext","Can't get a tree with TPC tracks !\n");
8041e7cf 197 fCurrentEvent++;//go to next dir
bed069a4 198 continue;
199 }
200
201 TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
202 if (!trackbranch) ////check if we got the branch
203 {//if not return with error
204 Error("ReadNext","Can't get a branch with TPC tracks !\n");
8041e7cf 205 fCurrentEvent++;//go to next dir
bed069a4 206 continue;
207 }
208 Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
209 Info("ReadNext","Found %d TPC tracks.",ntpctracks);
210 //Copy tracks to array
211
212 AliTPCtrack *iotrack=0;
213
214 for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
88cb7938 215 {
bed069a4 216 iotrack=new AliTPCtrack; //create new tracks
217 trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
218 tracktree->GetEvent(i); //stream track i to the iotrack
219 tarray->AddLast(iotrack); //put the track in the array
88cb7938 220 }
bed069a4 221
222 Double_t xk;
223 Double_t par[5];
224 Float_t phi, lam, pt;//angles and transverse momentum
225 Int_t label; //label of the current track
226
227 AliStack* stack = fRunLoader->Stack();
228 if (stack == 0x0)
16701d1b 229 {
bed069a4 230 Error("ReadNext","Can not get stack for current event",fCurrentEvent);
231 fCurrentEvent++;
5ee8b891 232 continue;
16701d1b 233 }
bed069a4 234 stack->Particles();
235
236 for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
237 {
238 iotrack = (AliTPCtrack*)tarray->At(i);
239 label = iotrack->GetLabel();
240
241 if (label < 0) continue;
3f745d47 242
243 if (CheckTrack(iotrack)) continue;
244
bed069a4 245 TParticle *p = (TParticle*)stack->Particle(label);
246
247 if(p == 0x0) continue; //if returned pointer is NULL
248 if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
249
250 if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
251 //if not take next partilce
252
253 AliHBTParticle* part = new AliHBTParticle(*p,i);
254 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
255 //if it does not delete it and take next good track
256
257// iotrack->PropagateTo(3.,0.0028,65.19);
258// iotrack->PropagateToVertex(36.66,1.2e-3);//orig
259 iotrack->PropagateToVertex(50.,0.0353);
260
261 iotrack->GetExternalParameters(xk,par); //get properties of the track
3f745d47 262
bed069a4 263 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
264 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
265 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
266 lam=par[3];
267 pt=1.0/TMath::Abs(par[4]);
268
269 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
270 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
271 Double_t tpz = pt * lam; //track z coordinate of momentum
272
273 Double_t mass = p->GetMass();
274 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
275
276 AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
277 if(Pass(track))//check if meets all criteria of any of our cuts
278 //if it does not delete it and take next good track
279 {
280 delete track;
281 delete part;
282 continue;
283 }
f5de2f09 284
285 if (fNTrackPoints > 0)
286 {
287 AliHBTTrackPoints* tpts = new AliHBTTrackPoints(fNTrackPoints,iotrack,fdR);
288 track->SetTrackPoints(tpts);
289 }
66d1d1a4 290 if ( fClusterMap )
291 {
292 AliHBTClusterMap* cmap = new AliHBTClusterMap(iotrack);
293 track->SetClusterMap(cmap);
294 }
295
bed069a4 296 fParticlesEvent->AddParticle(part);
297 fTracksEvent->AddParticle(track);
298 }
299
300 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
301 fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
302 fNEventsRead,fCurrentEvent,fCurrentDir);
bfb09ece 303 fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
bed069a4 304 fCurrentEvent++;
305 fNEventsRead++;
306 delete tarray;
307 return 0;
308 }while(fCurrentDir < GetNumberOfDirs());
309
310 delete tarray;
311 return 1;
312 }
313/********************************************************************/
314
315Int_t AliHBTReaderTPC::OpenNextSession()
316{
317 TString filename = GetDirName(fCurrentDir);
318 if (filename.IsNull())
319 {
320 DoOpenError("Can not get directory name");
321 return 1;
322 }
323 filename = filename +"/"+ fFileName;
324 fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
325 if( fRunLoader == 0x0)
326 {
327 DoOpenError("Can not open session.");
328 return 1;
329 }
330
331 fRunLoader->LoadHeader();
332 if (fRunLoader->GetNumberOfEvents() <= 0)
333 {
334 DoOpenError("There is no events in this directory.");
335 return 1;
336 }
337
338 if (fRunLoader->LoadKinematics())
339 {
340 DoOpenError("Error occured while loading kinematics.");
341 return 1;
342 }
343 fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
344 if ( fTPCLoader == 0x0)
345 {
346 DoOpenError("Exiting due to problems with opening files.");
347 return 1;
348 }
349 Info("OpenNextSession","________________________________________________________");
350 Info("OpenNextSession","Found %d event(s) in directory %s",
351 fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
352 Float_t mf;
353 if (fUseMagFFromRun)
354 {
355 fRunLoader->LoadgAlice();
356 mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
357 Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
358 fRunLoader->UnloadgAlice();
359 }
360 else
361 {
362 Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
363 if (fMagneticField == 0x0)
364 {
365 Fatal("OpenNextSession","Magnetic field can not be 0.");
366 return 1;//pro forma
16701d1b 367 }
bed069a4 368 mf = fMagneticField*10.;
369 }
370 AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
1b446896 371
16701d1b 372
bed069a4 373 fRunLoader->CdGAFile();
374 AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
375 if (!TPCParam)
376 {
377 TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
378 if (!TPCParam)
379 {
380 DoOpenError("TPC parameters have not been found !\n");
381 return 1;
382 }
383 }
384
385 if (fTPCLoader->LoadTracks())
386 {
387 DoOpenError("Error occured while loading TPC tracks.");
388 return 1;
389 }
88cb7938 390
bed069a4 391 fCurrentEvent = 0;
1b446896 392 return 0;
bed069a4 393}
394/********************************************************************/
1b446896 395
bed069a4 396void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
397{
398 // Does error display and clean-up in case error caught on Open Next Session
399
400 va_list ap;
401 va_start(ap,va_(fmt));
402 Error("OpenNextSession", va_(fmt), ap);
403 va_end(ap);
404
405 delete fRunLoader;
406 fRunLoader = 0x0;
407 fTPCLoader = 0x0;
408 fCurrentDir++;
409}
1b446896 410
3f745d47 411/********************************************************************/
412Bool_t AliHBTReaderTPC::CheckTrack(AliTPCtrack* t)
413{
414 //Performs check of the track
415 if ( (t->GetNumberOfClusters() > fNClustMax) || (t->GetNumberOfClusters() < fNClustMin) ) return kTRUE;
416
f5de2f09 417 Float_t chisqpercl = t->GetChi2()/((Double_t)t->GetNumberOfClusters());
418 if ( (chisqpercl < fNChi2PerClustMin) || (chisqpercl > fNChi2PerClustMax) ) return kTRUE;
419
3f745d47 420 Double_t cc[15];
421 t->GetCovariance(cc);
f5de2f09 422
423 if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE;
424 if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE;
425 if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE;
426 if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE;
427 if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
3f745d47 428
3f745d47 429
430 return kFALSE;
431
432}
433/********************************************************************/
434
435void AliHBTReaderTPC::SetNClustersRange(Int_t min,Int_t max)
436{
437 //sets range of Number Of Clusters that tracks have to have
438 fNClustMin = min;
439 fNClustMax = max;
440}
441/********************************************************************/
442
443void AliHBTReaderTPC::SetChi2PerCluserRange(Float_t min, Float_t max)
444{
445 //sets range of Chi2 per Cluster
446 fNChi2PerClustMin = min;
447 fNChi2PerClustMax = max;
448}
449/********************************************************************/
450
f5de2f09 451void AliHBTReaderTPC::SetC00Range(Float_t min, Float_t max)
452{
453 //Sets range of C00 parameter of covariance matrix of the track
454 //it defines uncertainty of the momentum
455 fC00Min = min;
456 fC00Max = max;
457}
458/********************************************************************/
459
460void AliHBTReaderTPC::SetC11Range(Float_t min, Float_t max)
461{
462 //Sets range of C11 parameter of covariance matrix of the track
463 //it defines uncertainty of the momentum
464 fC11Min = min;
465 fC11Max = max;
466}
467/********************************************************************/
468
469void AliHBTReaderTPC::SetC22Range(Float_t min, Float_t max)
470{
471 //Sets range of C22 parameter of covariance matrix of the track
472 //it defines uncertainty of the momentum
473 fC22Min = min;
474 fC22Max = max;
475}
476/********************************************************************/
477
478void AliHBTReaderTPC::SetC33Range(Float_t min, Float_t max)
479{
480 //Sets range of C33 parameter of covariance matrix of the track
481 //it defines uncertainty of the momentum
482 fC33Min = min;
483 fC33Max = max;
484}
485/********************************************************************/
486
3f745d47 487void AliHBTReaderTPC::SetC44Range(Float_t min, Float_t max)
488{
489 //Sets range of C44 parameter of covariance matrix of the track
490 //it defines uncertainty of the momentum
491 fC44Min = min;
492 fC44Max = max;
493}
494
1b446896 495/********************************************************************/
496/********************************************************************/
497/********************************************************************/
498
3f745d47 499/*
500void (AliTPCtrack* iotrack, Double_t curv)
501{
502 Double_t x[5];
503 iotrack->GetExternalPara
504 //xk is a
505 Double_t fP4=iotrack->GetC();
506 Double_t fP2=iotrack->GetEta();
507
508 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fP0, z1=fP1;
509 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
510 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
511
512 fP0 += dx*(c1+c2)/(r1+r2);
513 fP1 += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
514
515}
516*/
517