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