Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 : //
16 : //-------------------------------------------------------
17 : // Implementation of the TPC cosmic trackfit
18 : //
19 : // Origin: Xianguo Lu lu@physi.uni-heidelberg.de Xianguo.Lu@cern.ch
20 : //
21 : //=========================================================================================
22 : // Motivation:
23 : //
24 : // In the default reconstruction in the ALICE the cosmic tracks are found as two independent particles.
25 : //
26 : // In general any of subtracks can be used for the physics studies. In order to avoid the double counting,
27 : // the track from the upper hemisphere can be used.
28 : //
29 : // The momentum resolution is determined by the lever arm (1/L^2) and by the number of clusters
30 : // used for the track fitting (1/sqrt(Ncl)).
31 : // Combining/refitting the two segments together significantly better momentum resolution can be obtained.
32 : // sigma_{1/pt} ~ 8x10^-3 - defaul tracking (e.g only upper track)
33 : // sigma_{1/pt} ~ 8x10^-4 - combined tracking
34 : //===========================================================================================
35 : //
36 : // Interface/Implementation:
37 : // This class provide the functionality to combine two ESD tracks and perform the trackfit using all related track points.
38 : // Input: AliESDtrack's (or AliTPCseed's, expert use)
39 : // Output: the trackparams from the combined fit are stored as data member of this class, and they can be accessed via getter or copy.
40 : //
41 : //===========================================================================================
42 : //
43 : // Algorithm:
44 : // 1. Read in two ESDtracks (see function CombineESDtracks())
45 : // TPCseeds are got from ESDfriends
46 : // Cuts are applied on the following quantities: #TPC cluster, relation between global-y coordinate of the two tracks, lever arm.
47 : // 2. AliTPCCosmicUtils::CombinedFit is called and refit is performed using the TPCseeds from both tracks
48 : // 3. The fitting quality is analyzed and stored and can be accessed via function GetStatus(). The detailed judgement of the quality is recorded as "enum CombineStatus"
49 : //
50 : //
51 : //===========================================================================================
52 : // Algorithm numerical debugging:
53 : // The class can be used in the different debug/verbose level (see fDebugLevel)
54 : // Several intermediate variables can be stored in the trees, printout, or draw.
55 : // Given functionality (dumping of variables to the tree) was also used for the tuning of the pair
56 : // selection criterias, and for validation of the fit functionality.
57 : //
58 : // ----- Debug output:
59 : // for (debuglevel & 1)==1 && good fit, the following info saved:
60 : /*
61 : (*fStreamer)<<"TrackProp"<<
62 : "Tup.="<<fTrackparUp<< //AliExternalTrackParam at uppermost cluster obtained by upward propagation
63 : "Tlow.="<<fTrackparLow<< //AliExternalTrackParam at lowermost cluster obtained by downward propagation
64 : "icup.="<<&fInnerClusterUp<< //TVector3 position of the innermost cluster of the upper track
65 : "iclow.="<<&fInnerClusterLow<<
66 : "leverarm="<<fLeverArm<<
67 : "ncl="<<fFitNcls<< //number of clusters used in successful propagation
68 : "nmiss="<<fMissNcls<< //number of clusters failed in propagation, should always be 0 in this case.
69 : "chi2="<<fPreChi2<< //chi2/nfit
70 : "momup="<< momup << //momentum at uppermost cluster with upward propagation
71 : "momlow="<< momlow << //momentum at lowermost cluster with downward propagation
72 : "ptup="<< ptup <<
73 : "ptlow="<< ptlow <<
74 : "\n";
75 : */
76 : // for (debuglevel & 2)==1, debug info in AliTPCCosmicUtils::FitKernel saved
77 : //
78 : //===========================================================================================
79 : // Usage
80 : /*
81 : fCosmicTrackfit = new AliTPCCosmicTrackfit(debuglevel, "anystring");
82 :
83 : //order not important; will be internally ordered (potinters modified due to &) such that track0 is the upper one
84 : //kfit = kTRUE: good fit, kFALSE: bad fit
85 : const Bool_t kfit = fCosmicTrackfit->CombineESDtracks(esdtrack0, esdtrack1);
86 :
87 : //status = 0 for good fit (i.e. kfit=kTRUE), non-0 for bad fit (i.e. kfit=kFALSE), see "enum CombineStatus" definition in header file
88 : const Int_t status = fCosmicTrackfit->GetStatus();
89 : */
90 : //===========================================================================================
91 : // Efficiency
92 : //
93 : // for 2011 Feb. cosmic data nch=2 events, the kfit and status look like:
94 : /*
95 : kfit,status ( 0, 1): 16337 / 443901 = 3.680% //kFailGetTPCseeds
96 : kfit,status ( 0, 2): 3692 / 443901 = 0.832% //not both tracks have ncl > AliTPCCosmicUtils::NclsMin()
97 : kfit,status ( 0, 3): 6527 / 443901 = 1.470% //clusters in two tracks should be clearly separated in y, i.e. lowest cluster of upper track higher than highest cluster of lower track; otherwise fail
98 : kfit,status ( 0, 4): 7033 / 443901 = 1.584% //fLeverArm<CutLeverArm()
99 : kfit,status ( 0, 6): 4398 / 443901 = 0.991% //fail in propagation of at least one cluster
100 : kfit,status ( 0, 7): 508 / 443901 = 0.114% //chi2/nfit > MaxChi2()
101 : kfit,status ( 0, 8): 52 / 443901 = 0.012% //fail in geting impact parameters
102 : kfit,status ( 1, 0): 405354 / 443901 = 91.316% //i.e. 91% of nch=2 events are successfully fitted.
103 : */
104 : //
105 : //===========================================================================================
106 : // Resolution
107 : //
108 : // for muon momentum small than 20 GeV, energy loss in muon filter is visable when compaing fTrackparUp and fTrackparLow; energy loss estimated as 5 MeV/cm.
109 : // particle traversing muon filter can be rejected by requiring "fInnerClusterUp.fZ > -40 && fInnerClusterLow.fZ > -40"
110 : // momentum resolution is estimated by comparing the trackfit result by upward propagation through odd pad rows and that by downward propagation through even pad rows. Number of clusters used in this case is only half of that in normal usage.
111 : // RMS of log10 p = 0.01 at 10 GeV/c, 0.1 at 100 GeV/c, 0.5 at 1 TeV/c.
112 : // muon filter deteriorates momentum resolution by about +0.01 (absolute value).
113 : //
114 :
115 : #include <TAxis.h>
116 : #include <TCanvas.h>
117 : #include <TFile.h>
118 : #include <TGraph.h>
119 : #include <TTreeStream.h>
120 : #include <TVector3.h>
121 :
122 : #include "AliESDtrack.h"
123 : #include "AliESDfriendTrack.h"
124 : #include "AliTPCseed.h"
125 : #include "AliTrackerBase.h"
126 : #include "AliTrackPointArray.h"
127 :
128 : #include "AliTPCCosmicUtils.h"
129 : #include "AliTPCCosmicTrackfit.h"
130 :
131 : AliTPCCosmicTrackfit::AliTPCCosmicTrackfit(const Int_t dlev, const TString tag):
132 16 : fStreamer(0x0), fDebugLevel(dlev)
133 8 : , fSeedUp(0x0), fSeedLow(0x0), fTrackparUp(0x0), fTrackparLow(0x0), fIniTrackparUp(0x0), fIniTrackparLow(0x0)
134 8 : , fStatus(-999)
135 8 : , fKswap(0)
136 8 : , fInnerClusterUp(-999,-999,-999)
137 8 : , fInnerClusterLow(-999,-999,-999)
138 8 : , fLeverArm(-999)
139 8 : , fFitNcls(-999), fMissNcls(-999), fPreChi2(-999), fFitLeverArm(-999), fImpactD(-999), fImpactZ(-999)
140 8 : , fRowStartShift(-999)
141 8 : , fRowStep(-999)
142 8 : , fXMin(-999)
143 8 : , fXMax(-999)
144 16 : {
145 : //
146 : //Constructor
147 : //
148 :
149 8 : if(fDebugLevel>0)
150 0 : fStreamer = new TTreeSRedirector(Form("CosmicTrackfit_%s.root", tag.Data()));
151 :
152 8 : SetRow(0,1);
153 8 : SetX(-1e10, 1e10);
154 16 : }
155 :
156 : AliTPCCosmicTrackfit::~AliTPCCosmicTrackfit()
157 16 : {
158 : //
159 : //Destructor
160 : //
161 :
162 16 : delete fStreamer;
163 :
164 9 : delete fTrackparUp;
165 9 : delete fTrackparLow;
166 9 : delete fIniTrackparUp;
167 9 : delete fIniTrackparLow;
168 16 : }
169 :
170 : Bool_t AliTPCCosmicTrackfit::CombineESDtracks(AliESDtrack * &trk0, AliESDtrack *&trk1)
171 : {
172 : //
173 : //Get TPCseeds from the 2 ESDtracks, swap TPCseeds and ESDTracks (if necessary) according to y (0:upper 1:lower), perform trackfit using TPCseeds
174 : //if fStatus==0, i.e. combine is successful, swap of the ESDtracks is kept since pointer *& is used
175 : //
176 :
177 6 : IniCombineESDtracks();
178 :
179 3 : if(!GetTPCseeds(trk0, trk1)){
180 0 : return kFALSE;
181 : }
182 :
183 3 : CombineTPCseeds();
184 :
185 3 : if(fStatus == 0){
186 0 : if(fKswap){
187 0 : AliESDtrack * tmptrk = trk0;
188 0 : trk0 = trk1;
189 0 : trk1 = tmptrk;
190 0 : }
191 0 : return kTRUE;
192 : }
193 : else
194 3 : return kFALSE;
195 3 : }
196 :
197 : Bool_t AliTPCCosmicTrackfit::CombineTPCseedsFast(AliTPCseed * tpcseeds[], const AliExternalTrackParam * trkpars[])
198 : {
199 : //
200 : //do combined track fit for given TPC seeds and initial trackpar, [0]: upper, [1]: lower
201 : //
202 :
203 0 : IniCombineESDtracks();
204 :
205 0 : fSeedUp = tpcseeds[0];
206 0 : fSeedLow = tpcseeds[1];
207 :
208 0 : fTrackparUp = new AliExternalTrackParam(*(trkpars[0]));
209 0 : fTrackparLow = new AliExternalTrackParam(*(trkpars[1]));
210 :
211 0 : AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
212 0 : const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
213 : TTreeSRedirector * debugstreamer = 0x0;
214 0 : if(fDebugLevel&2){
215 0 : debugstreamer = fStreamer;
216 0 : }
217 :
218 0 : AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
219 :
220 0 : Update();
221 :
222 0 : if(fStatus == 0){
223 0 : return kTRUE;
224 : }
225 : else
226 0 : return kFALSE;
227 0 : }
228 :
229 : Bool_t AliTPCCosmicTrackfit::CombineTPCseeds(AliTPCseed * &seed0, AliTPCseed *&seed1)
230 : {
231 : //
232 : //same as AliTPCCosmicTrackfit::CombineESDtracks, except that the seeds are passed in from outside, which can be still unordered
233 : //if fStatus==0, i.e. combine is successful, swap of the TPCseeds is kept since pointer *& is used
234 : //
235 0 : IniCombineESDtracks();
236 :
237 0 : fSeedUp = seed0;
238 0 : fSeedLow = seed1;
239 :
240 0 : CombineTPCseeds();
241 :
242 0 : if(fStatus==0){
243 0 : if(fKswap){
244 0 : AliTPCseed * tmpseed = seed0;
245 0 : seed0 = seed1;
246 0 : seed1 = tmpseed;
247 0 : }
248 0 : return kTRUE;
249 : }
250 : else
251 0 : return kFALSE;
252 0 : }
253 :
254 : void AliTPCCosmicTrackfit::Print() const
255 : {
256 : //
257 : //print out variable values
258 : //
259 0 : printf("Status %2d NclsU %3d NclsD %3d ZinnerU %7.2f ZinnerD %7.2f LeverArm %7.2f\n", fStatus, fSeedUp->GetNumberOfClusters(), fSeedLow->GetNumberOfClusters(), fInnerClusterUp.Z(), fInnerClusterLow.Z(), fLeverArm);
260 0 : }
261 :
262 : TVector3 AliTPCCosmicTrackfit::ImpactParameter2D() const
263 : {
264 : //
265 : //calculate the 2D-impactparameter from (0,0)
266 : //
267 :
268 0 : const TVector3 p0(0,0,0);
269 0 : const TVector3 l1(fInnerClusterUp.X(), fInnerClusterUp.Y(), 0);
270 0 : const TVector3 l2(fInnerClusterLow.X(), fInnerClusterLow.Y(), 0);
271 :
272 0 : TVector3 vtx;
273 0 : AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
274 : return vtx;
275 0 : }
276 :
277 : TVector3 AliTPCCosmicTrackfit::ImpactParameter3D() const
278 : {
279 : //
280 : //calculate the 3D-impactparameter from (0,0,0)
281 : //
282 :
283 0 : const TVector3 p0(0,0,0);
284 0 : const TVector3 l1(fInnerClusterUp);
285 0 : const TVector3 l2(fInnerClusterLow);
286 :
287 0 : TVector3 vtx;
288 0 : AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
289 :
290 : //==========================
291 : /*
292 : printf("testing... \n");
293 : TVector3 tmpv = ImpactParameter2D();
294 : if(fabs(tmpv.Mag()-vtx.Pt())>1e-6){
295 : printf("strange!!! %e %e %e\n", tmpv.Mag(), vtx.Pt(), fabs(tmpv.Mag()-vtx.Pt()));
296 : vtx.Print();
297 : tmpv.Print();
298 : fInnerClusterUp.Print();
299 : fInnerClusterLow.Print();
300 : }
301 : */
302 : //==========================
303 : return vtx;
304 0 : }
305 :
306 : Double_t AliTPCCosmicTrackfit::MinPhi() const
307 : {
308 : //
309 : //the smaller phi of the two tracks w.r.t. horizon
310 : //
311 0 : Double_t fsp[] = {TMath::Abs(TMath::Sin(fTrackparUp->Phi())), TMath::Abs(TMath::Sin(fTrackparLow->Phi()))};;
312 0 : return TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
313 : }
314 : //===================================================================================================
315 : //===================================================================================================
316 :
317 : void AliTPCCosmicTrackfit::IniCombineESDtracks()
318 : {
319 : //
320 : //initialization, for reuse of the same AliTPCCosmicTrackfit instance
321 : //
322 :
323 6 : fSeedUp = 0x0;
324 3 : fSeedLow = 0x0;
325 4 : delete fTrackparUp;
326 4 : delete fTrackparLow;
327 3 : fTrackparUp = 0x0;
328 3 : fTrackparLow = 0x0;
329 :
330 4 : delete fIniTrackparUp;
331 4 : delete fIniTrackparLow;
332 3 : fIniTrackparUp = 0x0;
333 3 : fIniTrackparLow = 0x0;
334 :
335 3 : fStatus = 0;
336 3 : fKswap = kFALSE;
337 3 : }
338 :
339 : void AliTPCCosmicTrackfit::CombineTPCseeds()
340 : {
341 : //
342 : //do combined trackfit using TPCseeds
343 : //
344 :
345 : if(
346 9 : !CheckNcls()
347 6 : || !AnaSeeds() // fSeedUp/Low swapped here! check by runTest.C
348 6 : || !CheckLeverArm()
349 : )
350 : return;
351 :
352 : //AliExternalTrackParam object created
353 2 : fTrackparUp = AliTPCCosmicUtils::MakeSeed(fSeedUp);
354 2 : fTrackparLow = AliTPCCosmicUtils::MakeSeed(fSeedLow);
355 4 : if(!fTrackparUp || !fTrackparLow){
356 0 : fStatus = kFailMakeSeed;
357 0 : return;
358 : }
359 :
360 4 : fIniTrackparUp = new AliExternalTrackParam(*fTrackparUp);
361 4 : fIniTrackparLow = new AliExternalTrackParam(*fTrackparLow);
362 :
363 2 : AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
364 2 : const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
365 : TTreeSRedirector * debugstreamer = 0x0;
366 2 : if(fDebugLevel&2){
367 0 : debugstreamer = fStreamer;
368 0 : }
369 :
370 2 : AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
371 :
372 2 : Update();
373 :
374 : return;
375 5 : }
376 :
377 : void AliTPCCosmicTrackfit::Update()
378 : {
379 : //
380 : //Update variables depending on the fit result
381 : //
382 :
383 6 : if(fMissNcls || fFitNcls==0){
384 0 : fStatus = kFailPropagation;
385 0 : return;
386 : }
387 :
388 2 : fPreChi2 /= fFitNcls;
389 2 : if(fPreChi2>MaxChi2()){
390 2 : fStatus = kFailChi2;
391 2 : return;
392 : }
393 :
394 0 : if(fImpactD<0){
395 0 : fStatus = kFailImpact;
396 0 : return;
397 : }
398 :
399 0 : if( fStatus == 0 && (fDebugLevel&1) ){
400 0 : Double_t momup = fTrackparUp->P();
401 0 : Double_t momlow = fTrackparLow->P();
402 0 : Double_t ptup = fTrackparUp->Pt();
403 0 : Double_t ptlow = fTrackparLow->Pt();
404 :
405 0 : (*fStreamer)<<"TrackProp"<<
406 0 : "Tup.="<<fTrackparUp<<
407 0 : "Tlow.="<<fTrackparLow<<
408 0 : "icup.="<<&fInnerClusterUp<<
409 0 : "iclow.="<<&fInnerClusterLow<<
410 0 : "leverarm="<<fLeverArm<<
411 0 : "ncl="<<fFitNcls<<
412 0 : "nmiss="<<fMissNcls<<
413 0 : "chi2="<<fPreChi2<<
414 0 : "momup="<< momup <<
415 0 : "momlow="<< momlow <<
416 0 : "ptup="<< ptup <<
417 0 : "ptlow="<< ptlow <<
418 : "\n";
419 0 : }
420 2 : }
421 :
422 : Bool_t AliTPCCosmicTrackfit::CheckLeverArm()
423 : {
424 : //
425 : //if lever arm is too short, no need to use combined track fit.
426 : //On the other hand, short lever arm from two tracks mostly means they are fake pairs.
427 : //lever arm extents over one quadrant, e.g. (0,250)-(250,0): 250*sqrt(2)~350
428 : //
429 6 : if(fLeverArm<CutLeverArm()){
430 1 : fStatus = kFailLeverArm;
431 1 : return kFALSE;
432 : }
433 : else
434 2 : return kTRUE;
435 3 : }
436 :
437 : Bool_t AliTPCCosmicTrackfit::AnaSeeds()
438 : {
439 : //
440 : //swap seeds (if necessary) so that (y of fSeedUp) > (y of fSeedLow)
441 : //
442 :
443 : //---------------------------------- navigate through all clusters ----------------------------------
444 6 : AliTPCseed ** seeds[]={&fSeedUp, &fSeedLow};
445 :
446 : //min, max according to y
447 30 : TVector3 singlemin[2], singlemax[2];
448 18 : for(Int_t ii=0; ii<2; ii++){
449 6 : singlemin[ii].SetXYZ( 1e10, 1e10, 1e10);
450 6 : singlemax[ii].SetXYZ(-1e10, -1e10, -1e10);
451 : }
452 :
453 18 : for(Int_t itrk=0; itrk<2; itrk++){
454 1920 : for(Int_t irow=0; irow<AliTPCCosmicUtils::NRow(); irow++){
455 954 : const AliTPCclusterMI * cls = (*seeds[itrk])->GetClusterPointer(irow);
456 954 : if(!cls)
457 83 : continue;
458 :
459 871 : Float_t xyz[3]={-999,-999,-999};
460 871 : cls->GetGlobalXYZ(xyz);
461 871 : if(xyz[1]<singlemin[itrk].Y()){
462 456 : singlemin[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
463 456 : }
464 871 : if(xyz[1]>singlemax[itrk].Y()){
465 405 : singlemax[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
466 405 : }
467 871 : }
468 : }
469 :
470 : //--------------------------------
471 :
472 : //kpass true if y of the two seeds clearly separate: min of one > max of the other
473 : Bool_t kpass = kFALSE;
474 :
475 3 : fInnerClusterUp.SetXYZ(-999,-999,-999);
476 3 : fInnerClusterLow.SetXYZ(-999,-999,-999);
477 6 : TVector3 combinedmin, combinedmax;
478 3 : if(singlemin[0].Y()> singlemax[1].Y()){
479 0 : fInnerClusterUp = singlemin[0];
480 0 : fInnerClusterLow = singlemax[1];
481 :
482 : //no need to swap
483 0 : fKswap = kFALSE;
484 :
485 : kpass = kTRUE;
486 :
487 0 : combinedmax = singlemax[0];
488 0 : combinedmin = singlemin[1];
489 0 : }
490 3 : else if(singlemin[1].Y()> singlemax[0].Y()){
491 3 : fInnerClusterUp = singlemin[1];
492 3 : fInnerClusterLow = singlemax[0];
493 :
494 : //have to be swapped
495 3 : fKswap = kTRUE;
496 3 : AliTPCseed *tmp=*(seeds[0]);
497 3 : *(seeds[0])=*(seeds[1]);
498 3 : *(seeds[1])=tmp;
499 :
500 : kpass = kTRUE;
501 :
502 3 : combinedmax = singlemax[1];
503 3 : combinedmin = singlemin[0];
504 3 : }
505 : else
506 : kpass = kFALSE;
507 :
508 3 : const TVector3 comdelta = combinedmax-combinedmin;
509 6 : fLeverArm = comdelta.Pt();
510 :
511 3 : if(!kpass){
512 0 : fStatus = kFailSwapSeeds;
513 0 : return kFALSE;
514 : }
515 : else
516 3 : return kTRUE;
517 21 : }
518 :
519 : Bool_t AliTPCCosmicTrackfit::CheckNcls()
520 : {
521 : //
522 : //check number of clusters in TPCseed, for too small number MakeSeed will fail
523 : //
524 9 : if( fSeedUp->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() || fSeedLow->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() ){
525 0 : fStatus = kFailNclsMin;
526 0 : return kFALSE;
527 : }
528 : else
529 3 : return kTRUE;
530 3 : }
531 :
532 : Bool_t AliTPCCosmicTrackfit::GetTPCseeds(const AliESDtrack *trk0, const AliESDtrack *trk1)
533 : {
534 : //
535 : //Get TPC seeds from ESDfriendTrack
536 : //
537 6 : fSeedUp = AliTPCCosmicUtils::GetTPCseed(trk0);
538 3 : fSeedLow = AliTPCCosmicUtils::GetTPCseed(trk1);
539 :
540 6 : if(!fSeedUp || !fSeedLow){
541 0 : fStatus = kFailGetTPCseeds;
542 0 : return kFALSE;
543 : }
544 :
545 3 : return kTRUE;
546 3 : }
547 :
|
