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