d4c8e6bd2f7882ef2f0e7863732072d0c86e2328
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
1 #include "AliHBTReaderTPC.h"
2 //______________________________________________
3 //
4 // class AliHBTReaderTPC
5 //
6 // reader for TPC tracks
7 // needs galice.root
8 // just to shut up coding conventions checker
9 // 
10 // more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
11 // Piotr.Skowronski@cern.ch
12 //
13 ///////////////////////////////////////////////////////////////////////////
14 #include <TTree.h>
15 #include <TParticle.h>
16 #include <TH1.h>
17
18
19 #include <AliRun.h>
20 #include <AliLoader.h>
21 #include <AliStack.h>
22 #include <AliMagF.h>
23 #include <AliTPCtrack.h>
24 #include <AliTPCLoader.h>
25
26 #include "AliHBTEvent.h"
27 #include "AliHBTParticle.h"
28 #include "AliHBTTrackPoints.h"
29 #include "AliHBTClusterMap.h"
30
31
32 ClassImp(AliHBTReaderTPC)
33
34 AliHBTReaderTPC::AliHBTReaderTPC():
35  fFileName("galice.root"),
36  fRunLoader(0x0),
37  fTPCLoader(0x0),
38  fMagneticField(0.0),
39  fUseMagFFromRun(kTRUE),
40  fNTrackPoints(0),
41  fdR(0.0),
42  fClusterMap(kFALSE),
43  fNClustMin(0),
44  fNClustMax(150),
45  fNChi2PerClustMin(0.0),
46  fNChi2PerClustMax(10e5),
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),
55  fC44Min(0.0),
56  fC44Max(10e5)
57 {
58   //constructor
59   
60 }
61 /********************************************************************/
62
63 AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
64  fFileName(galicefilename),
65  fRunLoader(0x0),
66  fTPCLoader(0x0),
67  fMagneticField(0.0),
68  fUseMagFFromRun(kTRUE),
69  fNTrackPoints(0),
70  fdR(0.0),
71  fNClustMin(0),
72  fNClustMax(150),
73  fNChi2PerClustMin(0.0),
74  fNChi2PerClustMax(10e5),
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),
83  fC44Min(0.0),
84  fC44Max(10e5)
85 {
86   //constructor, 
87   //Defaults:
88 }
89 /********************************************************************/
90
91 AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
92  AliHBTReader(dirs), 
93  fFileName(galicefilename),
94  fRunLoader(0x0),
95  fTPCLoader(0x0),
96  fMagneticField(0.0),
97  fUseMagFFromRun(kTRUE),
98  fNTrackPoints(0),
99  fdR(0.0),
100  fClusterMap(kFALSE),
101  fNClustMin(0),
102  fNClustMax(150),
103  fNChi2PerClustMin(0.0),
104  fNChi2PerClustMax(10e5),
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),
113  fC44Min(0.0),
114  fC44Max(10e5)
115 {
116   //constructor, 
117   //Defaults:
118   //  galicefilename = ""  - this means: Do not open gAlice file - 
119   //                         just leave the global pointer untached
120   
121 }
122 /********************************************************************/
123 AliHBTReaderTPC::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 /********************************************************************/
150
151 AliHBTReaderTPC::~AliHBTReaderTPC()
152 {
153  //desctructor
154    if (AliHBTParticle::GetDebug()) 
155     {
156       Info("~AliHBTReaderTPC","deleting Run Loader");
157       AliLoader::SetDebug(AliHBTParticle::GetDebug());
158     }
159    
160    delete fRunLoader;
161    
162    if (AliHBTParticle::GetDebug()) 
163     {
164       Info("~AliHBTReaderTPC","deleting Run Loader Done");
165     }
166 }
167 /********************************************************************/
168
169 AliHBTReaderTPC& 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;
196  return *this; 
197
198 /********************************************************************/
199
200 void AliHBTReaderTPC::Rewind()
201 {
202 //Rewind reading to the beginning
203   delete fRunLoader;
204   fRunLoader = 0x0;
205   fCurrentDir = 0;
206   fNEventsRead= 0;
207   if (fTrackCounter) fTrackCounter->Reset();
208 }
209 /********************************************************************/
210
211 Int_t AliHBTReaderTPC::ReadNext()
212  {
213  //reads data and puts put to the particles and tracks objects
214  //reurns 0 if everything is OK
215  //
216   Info("Read","");
217   
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
222
223   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
224   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
225
226   fParticlesEvent->Reset();
227   fTracksEvent->Reset();
228
229   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
230    {
231       
232     if (fRunLoader == 0x0) 
233       if (OpenNextSession()) continue;//directory counter is increased inside in case of error
234
235     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
236      {
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
242      }
243      
244     Info("ReadNext","Reading Event %d",fCurrentEvent);
245     
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"); 
251          fCurrentEvent++;//go to next dir
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"); 
259         fCurrentEvent++;//go to next dir
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
269      {
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
274      }
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)
283      {
284        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
285        fCurrentEvent++;
286        continue;
287      }
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;
296        
297        if (CheckTrack(iotrack)) continue;//Checks the cuts on track parameters cov. mtx etc
298        
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
304        if(Rejected(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
305                                    //if not take next partilce
306
307        AliHBTParticle* part = new AliHBTParticle(*p,i);
308        if(Rejected(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
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
316        
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.);
331        if(Rejected(track))//check if meets all criteria of any of our cuts
332                     //if it does not delete it and take next good track
333         { 
334           delete track;
335           delete part;
336           continue;
337         }
338         
339        if (fNTrackPoints > 0) 
340         {
341           AliHBTTrackPoints* tpts = new AliHBTTrackPoints(fNTrackPoints,iotrack,fdR);
342           track->SetTrackPoints(tpts);
343         }
344        if (  fClusterMap ) 
345         {
346           AliHBTClusterMap* cmap = new AliHBTClusterMap(iotrack);
347           track->SetClusterMap(cmap);
348         }
349     
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);
357     fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
358     fCurrentEvent++;
359     fNEventsRead++;
360     delete tarray;
361     return 0;
362    }while(fCurrentDir < GetNumberOfDirs());
363
364   delete tarray;
365   return 1;
366  }
367 /********************************************************************/
368
369 Int_t AliHBTReaderTPC::OpenNextSession()
370 {
371 //Opens session thats from fCurrentDir 
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;
379   fRunLoader = AliRunLoader::Open(filename,AliConfig::GetDefaultEventFolderName());
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
422       }
423      mf = fMagneticField*10.;
424    }
425   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
426
427
428   if (fTPCLoader->LoadTracks())
429    {
430      DoOpenError("Error occured while loading TPC tracks.");
431      return 1;
432    }
433   
434   fCurrentEvent = 0;
435   return 0;
436 }
437 /********************************************************************/
438
439 void 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 }
453
454 /********************************************************************/
455 Bool_t AliHBTReaderTPC::CheckTrack(AliTPCtrack* t) const
456 {
457   //Performs check of the track
458   if ( (t->GetNumberOfClusters() > fNClustMax) || (t->GetNumberOfClusters() < fNClustMin) ) return kTRUE;
459
460   Float_t chisqpercl = t->GetChi2()/((Double_t)t->GetNumberOfClusters());
461   if ( (chisqpercl < fNChi2PerClustMin) || (chisqpercl > fNChi2PerClustMax) ) return kTRUE;
462   
463   Double_t cc[15];
464   t->GetExternalCovariance(cc);
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;
471   
472   
473   return kFALSE;
474   
475 }
476 /********************************************************************/
477
478 void 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
486 void 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
494 void 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
503 void 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
512 void 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
521 void 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
530 void 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
538 /********************************************************************/
539 /********************************************************************/
540 /********************************************************************/
541
542 /*
543 void (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