123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282 |
- /// \class MpdTpcKalmanTrack
- ///
- /// Kalman filter track object for the MPD central detector
- /// Track parameters: 0: RPhi - coordinate in R-Phi direction
- /// 1: Z - longitudinal coordinate
- /// 2: Phi - local azimuthal angle
- /// 3: Theta - local dip angle (angle w.r.t. the transverse plane)
- /// 4: q/Pt - signed inverse Pt
- /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
- #include "MpdTpcKalmanTrack.h"
- #include "MpdKalmanFilter.h"
- #include "MpdKalmanGeoScheme.h"
- #include "MpdKalmanHit.h"
- #include "MpdKalmanTrack.h"
- //#include "MpdTpcHit.h" // debug
- //#include "TpcPoint.h" // debug
- #include <TMatrixD.h>
- #include <TMatrixDSym.h>
- #include <TClonesArray.h>
- #include <TVector3.h>
- #include <TMath.h>
- #include <Riostream.h>
- using namespace std;
- //__________________________________________________________________________
- MpdTpcKalmanTrack::MpdTpcKalmanTrack()
- : MpdKalmanTrack(),
- fTrHits(new TClonesArray("MpdKalmanHit"))
- {
- /// Default constructor
- }
- //__________________________________________________________________________
- MpdTpcKalmanTrack::MpdTpcKalmanTrack(MpdKalmanHit *hitOut, MpdKalmanHit *hitIn,
- TVector3 &vertex, Double_t pt, Double_t *params)
- : MpdKalmanTrack(0., vertex),
- fTrHits(new TClonesArray("MpdKalmanHit",70))
- {
- /// Constructor from 2 hits
- Double_t rOut = hitOut->GetPos();
- Double_t phiOut = hitOut->GetMeas(0) / rOut;
- Double_t rIn = hitIn->GetPos();
- Double_t phiIn = hitIn->GetMeas(0) / rIn;
- Double_t parOut[4] = {rOut, phiOut, 0., 0.};
- Double_t parIn[4] = {rIn, phiIn, 0., 0.};
- EvalParams(hitOut, hitIn, parOut, parIn, pt, params);
- }
- //__________________________________________________________________________
- MpdTpcKalmanTrack::MpdTpcKalmanTrack(MpdKalmanHit *hitOut, MpdKalmanHit *hitIn,
- TVector3 &vertex, TVector3 &posOut, TVector3 &posIn,
- Double_t pt, Double_t *params)
- : MpdKalmanTrack(0., vertex),
- fTrHits(new TClonesArray("MpdKalmanHit",70))
- {
- /// Constructor from 2 hits and 2 TVector3's (modular geometry)
- Double_t rOut = posOut.Pt(), phiOut = posOut.Phi();
- Double_t rIn = posIn.Pt(), phiIn = posIn.Phi();
- Double_t parOut[4] = {rOut, phiOut, 0., 0.};
- Double_t parIn[4] = {rIn, phiIn, 0., 0.};
- EvalParams(hitOut, hitIn, parOut, parIn, pt, params);
- }
- //__________________________________________________________________________
- MpdTpcKalmanTrack::MpdTpcKalmanTrack (const MpdTpcKalmanTrack& track)
- : MpdKalmanTrack(track),
- fTrHits(new TClonesArray("MpdKalmanHit",70))
- {
- ///copy constructor
- Int_t nHits = track.fTrHits->GetEntriesFast();
- for (Int_t i = 0; i < nHits; ++i) {
- MpdKalmanHit *hit = (MpdKalmanHit*) (track.fTrHits->UncheckedAt(i));
- new ((*fTrHits)[i]) MpdKalmanHit(*hit);
- }
- if (track.fHits == 0x0) {
- for (Int_t i = 0; i < nHits; ++i) {
- fHits->Add(fTrHits->UncheckedAt(i));
- }
- }
- }
- //__________________________________________________________________________
- MpdTpcKalmanTrack & MpdTpcKalmanTrack::operator=(const MpdTpcKalmanTrack& track)
- {
- /// Asignment operator
- // check assignement to self
- if (this == &track) return *this;
- // base class assignement
- MpdKalmanTrack::operator=(track);
- if (fTrHits) { fTrHits->Delete(); delete fTrHits; }
- fTrHits = new TClonesArray("MpdKalmanHit",70);
- Int_t nHits = track.fTrHits->GetEntriesFast();
- for (Int_t i = 0; i < nHits; ++i) {
- MpdKalmanHit *hit = (MpdKalmanHit*)(track.fTrHits->UncheckedAt(i));
- new ((*fTrHits)[i]) MpdKalmanHit(*hit);
- }
- if (track.fHits == 0x0) {
- for (Int_t i = 0; i < nHits; ++i) {
- fHits->Add(fTrHits->UncheckedAt(i));
- }
- }
- return *this;
- }
- //__________________________________________________________________________
- MpdTpcKalmanTrack::~MpdTpcKalmanTrack()
- {
- /// Destructor
- if (fTrHits) fTrHits->Delete();
- delete fTrHits;
- fTrHits = NULL;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanTrack::Reset()
- {
- /// Reset track
- if (fTrHits) fTrHits->Delete();
- delete fTrHits;
- fTrHits = NULL;
- Clear();
- }
- //__________________________________________________________________________
- void MpdTpcKalmanTrack::EvalParams(MpdKalmanHit *hitOut, const MpdKalmanHit *hitIn,
- Double_t *parOut, Double_t *parIn, Double_t pt, Double_t *params)
- {
- /// Evaluate track params
- Double_t rIn = parIn[0], phiIn = parIn[1];
- Double_t rOut = parOut[0], phiOut = parOut[1];
- Double_t xIn = rIn * TMath::Cos(phiIn);
- Double_t yIn = rIn * TMath::Sin(phiIn);
- Double_t xOut = rOut * TMath::Cos(phiOut);
- Double_t yOut = rOut * TMath::Sin(phiOut);
- parOut[2] = xOut;
- parOut[3] = yOut;
- parIn[2] = xIn;
- parIn[3] = yIn;
- fPos = fPosNew = rOut;
- (*fParam)(0,0) = rOut * phiOut; // R-Phi coordinate
- //(*fParam)(1,0) = posOut.Z(); // Z-coordinate
- (*fParam)(1,0) = hitOut->GetMeas(1); // Z-coordinate
- Double_t phi0 = TMath::ATan2 (yOut-yIn, xOut-xIn); // Track Phi - approximation
- (*fParam)(2,0) = phi0 + TMath::Sign (params[2]/2, pt);
- //(*fParam)(3,0) = TMath::ATan2 (hitOut->GetMeas(1)-hitIn->GetMeas(1), rOut-rIn); // Track Theta
- (*fParam)(3,0) = TMath::ATan2 (hitOut->GetMeas(1)-hitIn->GetMeas(1), TMath::Abs(params[0])); // Track Theta
- (*fParam)(4,0) = 1. / pt; // q/Pt
- // Debugging
- /*
- TClonesArray *hits = (TClonesArray*) FairRootManager::Instance()->GetObject("TpcHit");
- TClonesArray *points = (TClonesArray*) FairRootManager::Instance()->GetObject("TpcPoint");
- MpdTpcHit *h = (MpdTpcHit*) hits->UncheckedAt(hitOut->GetIndex());
- TpcPoint *p = (TpcPoint*) points->UncheckedAt(h->GetRefIndex());
- cout << " Phi angle: " << TMath::ATan2(p->GetPy(),p->GetPx()) << " " << phi0 << " " << (*fParam)(2,0) << endl;
- */
-
- *fParamNew = *fParam;
- fWeight->Zero();
- EvalCovar(hitOut, hitIn, parOut, parIn);
-
- fHits->Add(hitOut);
- fLastLay = hitOut->GetLayer(); // layer number
- }
- //__________________________________________________________________________
- void MpdTpcKalmanTrack::EvalCovar(const MpdKalmanHit *hitOut, const MpdKalmanHit *hitIn,
- Double_t *parOut, Double_t *parIn)
- {
- /// Evaluate covariance matrix for track seed
- (*fWeight)(0,0) = hitOut->GetErr(0) * hitOut->GetErr(0); // <RphiRphi>
- (*fWeight)(0,0) *= 4.; // extra factor of 2
- (*fWeight)(1,1) = hitOut->GetErr(1) * hitOut->GetErr(1); // <zz>
- Double_t phiOut = parOut[1], phiIn = parIn[1];
- Double_t dx = parOut[2] - parIn[2], dy = parOut[3] - parIn[3];
- Double_t dist2 = dx * dx + dy * dy;
- Double_t sinPhi = TMath::Sin ((*fParam)(2,0));
- Double_t cosPhi = TMath::Cos ((*fParam)(2,0));
- Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
- Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
- (*fWeight)(2,2) = (pOut * pOut + pIn * pIn) / dist2 * (*fWeight)(0,0); // <PhiPhi>
- (*fWeight)(2,2) *= 4.; // extra factor of 2
- Double_t tanThe = TMath::Tan((*fParam)(3,0));
- //31.12.2018 Double_t dRad = parOut[0] - parIn[0];
- Double_t dRad = TMath::Sqrt (dist2);
- Double_t denom = dRad * (1.+tanThe*tanThe);
- (*fWeight)(3,3) = (*fWeight)(1,1) * 2. / denom / denom; // <TheThe>
- //(*fWeight)(3,3) = (*fWeight)(1,1) * 20. / denom / denom; // <TheThe>
- //(*fWeight)(4,4) = ((*fParam)(4,0)*0.5) * ((*fParam)(4,0)*0.5); // error 50%
- //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
- (*fWeight)(4,4) = ((*fParam)(4,0)*2.) * ((*fParam)(4,0)*2.); // error 200%
- //for (Int_t jj = 0; jj < 5; ++jj) cout << TMath::Sqrt((*fWeight)(jj,jj)) << " ";
- //cout << endl;
-
- //fWeight->Print();
- //fWeight->Invert(); // weight matrix
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(fWeight->GetMatrixArray(), 5, 5, 5, iok);
- //fWeight->Print();
- }
- //__________________________________________________________________________
- void MpdTpcKalmanTrack::StartBack()
- {
- /// Prepare for back tracing
- //#pragma omp critical
- {
- fHits->Sort(); // in descending order in radius
- }
- fChi2 = 0.;
- //fTrackDir = kInward;
- //if (fLastLay == 50) { cout << " back" << endl; Weight2Cov()->Print(); exit(0); }
- Int_t nHits = GetNofHits();
- if (nHits == 1) return;
- nHits *= nHits;
- //*
- for (Int_t i = 0; i < 5; ++i) {
- for (Int_t j = i; j < 5; ++j) {
- //if (j == i) (*fWeight)(i,j) /= 100.;
- if (j == i) (*fWeight)(i,j) /= nHits;
- else (*fWeight)(i,j) = (*fWeight)(j,i) = 0.;
- }
- }
- //*/
- //fWeight->Zero();
- //EvalCovar((MpdKalmanHitR*)fHits->UncheckedAt(0),(MpdKalmanHitR*)fHits->UncheckedAt(1));
- }
- //__________________________________________________________________________
- Int_t MpdTpcKalmanTrack::Compare(const TObject* track) const
- {
- /// "Compare" function to sort in descending order in pt
- MpdKalmanTrack *trackKF = (MpdKalmanTrack*) track;
- Double_t pt = 1. / TMath::Abs(trackKF->GetParam(4));
- Double_t ptthis = 1. / TMath::Abs((*fParam)(4,0));
- if (ptthis < pt) return 1;
- else if (ptthis > pt) return -1;
- return 0;
- }
- //__________________________________________________________________________
- Bool_t MpdTpcKalmanTrack::GetRecoQuality(Double_t dist, Double_t percentage)
- {
- /// returns kTRUE if number of hits closer to boundaries than dist divided by nHits is larger than percentage
- MpdKalmanHit *hit = NULL;
- Int_t nCloseHits = 0;
- Int_t nTrHits = fTrHits->GetEntries();
- if (nTrHits == 0) {cout << "NO KALMAN HITS ARE FOUND" << endl; return kTRUE;}
- for (Int_t itr = 0; itr < nTrHits; itr++)
- {
- hit = (MpdKalmanHit*) fTrHits->UncheckedAt(itr);
- if (TMath::Abs(hit->GetEdge()) < dist) nCloseHits++;
- }
- if ((Double_t)(nCloseHits/(Double_t)nTrHits) > percentage) return kTRUE;
- else return kFALSE;
- }
- ClassImp(MpdTpcKalmanTrack)
|