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