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