1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------
17 // Implementation of the TPC cosmic trackfit
19 // Origin: Xianguo Lu lu@physi.uni-heidelberg.de Xianguo.Lu@cern.ch
21 //=========================================================================================
24 // In the default reconstruction in the ALICE the cosmic tracks are found as two independent particles.
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.
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 //===========================================================================================
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.
41 //===========================================================================================
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"
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.
58 // ----- Debug output:
59 // for (debuglevel & 1)==1 && good fit, the following info saved:
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
76 // for (debuglevel & 2)==1, debug info in AliTPCCosmicUtils::FitKernel saved
78 //===========================================================================================
81 fCosmicTrackfit = new AliTPCCosmicTrackfit(debuglevel, "anystring");
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);
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();
90 //===========================================================================================
93 // for 2011 Feb. cosmic data nch=2 events, the kfit and status look like:
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.
105 //===========================================================================================
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).
119 #include <TTreeStream.h>
120 #include <TVector3.h>
122 #include "AliESDtrack.h"
123 #include "AliESDfriendTrack.h"
124 #include "AliTPCseed.h"
125 #include "AliTrackerBase.h"
126 #include "AliTrackPointArray.h"
128 #include "AliTPCCosmicUtils.h"
129 #include "AliTPCCosmicTrackfit.h"
131 AliTPCCosmicTrackfit::AliTPCCosmicTrackfit(const Int_t dlev, const TString tag):
132 fStreamer(0x0), fDebugLevel(dlev)
133 , fSeedUp(0x0), fSeedLow(0x0), fTrackparUp(0x0), fTrackparLow(0x0), fIniTrackparUp(0x0), fIniTrackparLow(0x0)
136 , fInnerClusterUp(-999,-999,-999)
137 , fInnerClusterLow(-999,-999,-999)
139 , fFitNcls(-999), fMissNcls(-999), fPreChi2(-999), fFitLeverArm(-999), fImpactD(-999), fImpactZ(-999)
140 , fRowStartShift(-999)
150 fStreamer = new TTreeSRedirector(Form("CosmicTrackfit_%s.root", tag.Data()));
156 AliTPCCosmicTrackfit::~AliTPCCosmicTrackfit()
166 delete fIniTrackparUp;
167 delete fIniTrackparLow;
170 Bool_t AliTPCCosmicTrackfit::CombineESDtracks(AliESDtrack * &trk0, AliESDtrack *&trk1)
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
177 IniCombineESDtracks();
179 if(!GetTPCseeds(trk0, trk1)){
187 AliESDtrack * tmptrk = trk0;
197 Bool_t AliTPCCosmicTrackfit::CombineTPCseedsFast(AliTPCseed * tpcseeds[], const AliExternalTrackParam * trkpars[])
200 //do combined track fit for given TPC seeds and initial trackpar, [0]: upper, [1]: lower
203 IniCombineESDtracks();
205 fSeedUp = tpcseeds[0];
206 fSeedLow = tpcseeds[1];
208 fTrackparUp = new AliExternalTrackParam(*(trkpars[0]));
209 fTrackparLow = new AliExternalTrackParam(*(trkpars[1]));
211 AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
212 const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
213 TTreeSRedirector * debugstreamer = 0x0;
215 debugstreamer = fStreamer;
218 AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
229 Bool_t AliTPCCosmicTrackfit::CombineTPCseeds(AliTPCseed * &seed0, AliTPCseed *&seed1)
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
235 IniCombineESDtracks();
244 AliTPCseed * tmpseed = seed0;
254 void AliTPCCosmicTrackfit::Print() const
257 //print out variable values
259 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);
262 TVector3 AliTPCCosmicTrackfit::ImpactParameter2D() const
265 //calculate the 2D-impactparameter from (0,0)
268 const TVector3 p0(0,0,0);
269 const TVector3 l1(fInnerClusterUp.X(), fInnerClusterUp.Y(), 0);
270 const TVector3 l2(fInnerClusterLow.X(), fInnerClusterLow.Y(), 0);
273 AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
277 TVector3 AliTPCCosmicTrackfit::ImpactParameter3D() const
280 //calculate the 3D-impactparameter from (0,0,0)
283 const TVector3 p0(0,0,0);
284 const TVector3 l1(fInnerClusterUp);
285 const TVector3 l2(fInnerClusterLow);
288 AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
290 //==========================
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()));
298 fInnerClusterUp.Print();
299 fInnerClusterLow.Print();
302 //==========================
306 Double_t AliTPCCosmicTrackfit::MinPhi() const
309 //the smaller phi of the two tracks w.r.t. horizon
311 Double_t fsp[] = {TMath::Abs(TMath::Sin(fTrackparUp->Phi())), TMath::Abs(TMath::Sin(fTrackparLow->Phi()))};;
312 return TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
314 //===================================================================================================
315 //===================================================================================================
317 void AliTPCCosmicTrackfit::IniCombineESDtracks()
320 //initialization, for reuse of the same AliTPCCosmicTrackfit instance
330 delete fIniTrackparUp;
331 delete fIniTrackparLow;
332 fIniTrackparUp = 0x0;
333 fIniTrackparLow = 0x0;
339 void AliTPCCosmicTrackfit::CombineTPCseeds()
342 //do combined trackfit using TPCseeds
347 || !AnaSeeds() // fSeedUp/Low swapped here! check by runTest.C
352 //AliExternalTrackParam object created
353 fTrackparUp = AliTPCCosmicUtils::MakeSeed(fSeedUp);
354 fTrackparLow = AliTPCCosmicUtils::MakeSeed(fSeedLow);
355 if(!fTrackparUp || !fTrackparLow){
356 fStatus = kFailMakeSeed;
360 fIniTrackparUp = new AliExternalTrackParam(*fTrackparUp);
361 fIniTrackparLow = new AliExternalTrackParam(*fTrackparLow);
363 AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
364 const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
365 TTreeSRedirector * debugstreamer = 0x0;
367 debugstreamer = fStreamer;
370 AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
377 void AliTPCCosmicTrackfit::Update()
380 //Update variables depending on the fit result
383 if(fMissNcls || fFitNcls==0){
384 fStatus = kFailPropagation;
388 fPreChi2 /= fFitNcls;
389 if(fPreChi2>MaxChi2()){
395 fStatus = kFailImpact;
399 if( fStatus == 0 && (fDebugLevel&1) ){
400 Double_t momup = fTrackparUp->P();
401 Double_t momlow = fTrackparLow->P();
402 Double_t ptup = fTrackparUp->Pt();
403 Double_t ptlow = fTrackparLow->Pt();
405 (*fStreamer)<<"TrackProp"<<
406 "Tup.="<<fTrackparUp<<
407 "Tlow.="<<fTrackparLow<<
408 "icup.="<<&fInnerClusterUp<<
409 "iclow.="<<&fInnerClusterLow<<
410 "leverarm="<<fLeverArm<<
412 "nmiss="<<fMissNcls<<
415 "momlow="<< momlow <<
422 Bool_t AliTPCCosmicTrackfit::CheckLeverArm()
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
429 if(fLeverArm<CutLeverArm()){
430 fStatus = kFailLeverArm;
437 Bool_t AliTPCCosmicTrackfit::AnaSeeds()
440 //swap seeds (if necessary) so that (y of fSeedUp) > (y of fSeedLow)
443 //---------------------------------- navigate through all clusters ----------------------------------
444 AliTPCseed ** seeds[]={&fSeedUp, &fSeedLow};
446 //min, max according to y
447 TVector3 singlemin[2], singlemax[2];
448 for(Int_t ii=0; ii<2; ii++){
449 singlemin[ii].SetXYZ( 1e10, 1e10, 1e10);
450 singlemax[ii].SetXYZ(-1e10, -1e10, -1e10);
453 for(Int_t itrk=0; itrk<2; itrk++){
454 for(Int_t irow=0; irow<AliTPCCosmicUtils::NRow(); irow++){
455 const AliTPCclusterMI * cls = (*seeds[itrk])->GetClusterPointer(irow);
459 Float_t xyz[3]={-999,-999,-999};
460 cls->GetGlobalXYZ(xyz);
461 if(xyz[1]<singlemin[itrk].Y()){
462 singlemin[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
464 if(xyz[1]>singlemax[itrk].Y()){
465 singlemax[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
470 //--------------------------------
472 //kpass true if y of the two seeds clearly separate: min of one > max of the other
473 Bool_t kpass = kFALSE;
475 fInnerClusterUp.SetXYZ(-999,-999,-999);
476 fInnerClusterLow.SetXYZ(-999,-999,-999);
477 TVector3 combinedmin, combinedmax;
478 if(singlemin[0].Y()> singlemax[1].Y()){
479 fInnerClusterUp = singlemin[0];
480 fInnerClusterLow = singlemax[1];
487 combinedmax = singlemax[0];
488 combinedmin = singlemin[1];
490 else if(singlemin[1].Y()> singlemax[0].Y()){
491 fInnerClusterUp = singlemin[1];
492 fInnerClusterLow = singlemax[0];
496 AliTPCseed *tmp=*(seeds[0]);
497 *(seeds[0])=*(seeds[1]);
502 combinedmax = singlemax[1];
503 combinedmin = singlemin[0];
508 const TVector3 comdelta = combinedmax-combinedmin;
509 fLeverArm = comdelta.Pt();
512 fStatus = kFailSwapSeeds;
519 Bool_t AliTPCCosmicTrackfit::CheckNcls()
522 //check number of clusters in TPCseed, for too small number MakeSeed will fail
524 if( fSeedUp->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() || fSeedLow->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() ){
525 fStatus = kFailNclsMin;
532 Bool_t AliTPCCosmicTrackfit::GetTPCseeds(const AliESDtrack *trk0, const AliESDtrack *trk1)
535 //Get TPC seeds from ESDfriendTrack
537 fSeedUp = AliTPCCosmicUtils::GetTPCseed(trk0);
538 fSeedLow = AliTPCCosmicUtils::GetTPCseed(trk1);
540 if(!fSeedUp || !fSeedLow){
541 fStatus = kFailGetTPCseeds;