diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 468e2f35e1b8cb744e4fc8c9261b6c2771fe572f..c707cea2e325b5c5013754013d629695f2bd5b50 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -124,7 +124,7 @@ L1Algo/L1TrackParFit.cxx
L1Algo/L1Event.cxx
L1Algo/L1EventMatch.cxx
L1Algo/L1MCEvent.cxx
-
+L1Algo/L1Fit.cxx
CbmL1MCTrack.cxx
ParticleFinder/CbmL1PFFitter.cxx
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index bce8567cdc70c55de9848d87317cce47c0db26cc..e579a0976dfa4c65cd8a29ac375d28e853507c7f 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1502,7 +1502,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
StsHitsLocal.push_back(algo->vRecoHits[start_hit++]);
}
- t.mass = 0.1395679; // pion mass
+ t.mass = algo->fDefaultMass; // pion mass
t.is_electron = 0;
t.SetId(vRTracksCur.size());
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index e83a706d826cec3808a833196447de6dfab56fcf..2aba2178443dcd689c8533d3f010c6d6c5118549 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -13,33 +13,23 @@
*
*====================================================================
*/
-#include "CbmL1.h"
-
+#include "CbmKF.h"
+#include "CbmKFMath.h"
#include "CbmKFTrack.h" // for vertex pulls
+#include "CbmL1.h"
#include "CbmL1Constants.h"
-#include "FairTrackParam.h" // for vertex pulls
-#include "L1Algo/L1AddMaterial.h" // for vertex pulls
-#include "L1Algo/L1Algo.h"
-#include "L1Algo/L1Extrapolation.h" // for vertex pulls
-
+#include "CbmL1Counters.h"
+#include "CbmMatch.h"
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
-
-#include "CbmTrdHit.h"
-#include "CbmTrdPoint.h"
-
-#include "CbmTofHit.h"
-#include "CbmTofPoint.h"
-
-#include "CbmKF.h"
-#include "CbmKFMath.h"
-
#include "CbmStsSetup.h"
#include "CbmStsStation.h"
+#include "CbmTofHit.h"
+#include "CbmTofPoint.h"
+#include "CbmTrdHit.h"
+#include "CbmTrdPoint.h"
-#include "CbmMatch.h"
-
-#include "CbmL1Counters.h"
+#include "FairTrackParam.h" // for vertex pulls
#include <TFile.h>
@@ -48,6 +38,10 @@
#include <map>
#include <vector>
+#include "L1Algo/L1Algo.h"
+#include "L1Algo/L1Extrapolation.h" // for vertex pulls
+#include "L1Algo/L1Fit.h" // for vertex pulls
+
using std::cout;
using std::endl;
using std::ios;
@@ -1547,6 +1541,10 @@ void CbmL1::TrackFitPerformance() {
static bool first_call = 1;
+ L1Fit fit;
+ fit.SetParticleMass(algo->GetDefaultParticleMass());
+
+
if (first_call) {
first_call = 0;
@@ -1861,9 +1859,8 @@ void CbmL1::TrackFitPerformance() {
&& (dir * (mc.z - algo->vStations[iSta].z[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
- L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp);
- if (iSta + dir == NMvdStations - 1)
- L1AddPipeMaterial(trPar, trPar.qp);
+ fit.L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp, 1);
+ if (iSta + dir == NMvdStations - 1) fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
}
}
if (mc.z != trPar.z[0]) continue;
@@ -1941,21 +1938,13 @@ void CbmL1::TrackFitPerformance() {
trPar.x - trPar.tx * dz, trPar.y - trPar.ty * dz, B[0]);
fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
- fvec mass2 = 0.1395679f * 0.1395679f;
L1Extrapolate(trPar, algo->vStations[iSta].z[0], trPar.qp, fld);
- L1AddMaterial(trPar,
- algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y),
- trPar.qp);
- EnergyLossCorrection(
- trPar,
- mass2,
- algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y),
- trPar.qp,
- fvec(1.f));
+ fit.L1AddMaterial(trPar, algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y), trPar.qp, 1);
+ fit.EnergyLossCorrection(trPar, algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y), trPar.qp, fvec(1.f),
+ fvec(1.f));
if (iSta + dir == NMvdStations - 1) {
- L1AddPipeMaterial(trPar, trPar.qp);
- EnergyLossCorrection(
- trPar, mass2, PipeRadThick, trPar.qp, fvec(1.f));
+ fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
+ fit.EnergyLossCorrection(trPar, fit.PipeRadThick, trPar.qp, fvec(1.f), fvec(1.f));
}
B[2] = B[1];
z[2] = z[1];
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 89c594f930bf89313ce5ea5212064fe26b6a3156..215fa11a5e49ee69f52872e5661e4ec67a153269 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -244,9 +244,22 @@ public:
L1Algo(const L1Algo&) = delete;
L1Algo operator=(const L1Algo&) = delete;
+ /// set a default particle mass for the track fit
+ /// it is used during reconstruction
+ /// for the multiple scattering and energy loss estimation
+ void SetDefaultParticleMass(float mass) { fDefaultMass = mass; }
+
+ /// get default particle mass
+ float GetDefaultParticleMass() const { return fDefaultMass; }
+
+ /// get default particle mass squared
+ float GetDefaultParticleMass2() const { return fDefaultMass * fDefaultMass; }
+
static const int nTh = 1;
static const int nSta = 25;
+ float fDefaultMass = 0.10565800; // muon mass
+
L1Vector<L1Triplet> TripletsLocal1[nSta][nTh];
L1Vector<L1Branch> CandidatesTrack[nTh];
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 978b851f231ccdaab4c54ff4fc0e134e7a5f0c5b..64466fa09f2111c2fe04a7b1f6e52ae780dbc1f5 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -14,11 +14,11 @@
*
*/
-#include "L1AddMaterial.h"
#include "L1Algo.h"
#include "L1Branch.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
+#include "L1Fit.h"
#include "L1HitPoint.h"
#include "L1Track.h"
#include "L1TrackPar.h"
@@ -142,6 +142,15 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
L1FieldValue m_B,
m_B_global _fvecalignment; // field for the next extrapolations
+ L1Fit fit;
+
+ if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
+ fit.SetParticleMass(fDefaultMass); // muon
+ }
+ else {
+ fit.SetParticleMass(0.000511f); // electron
+ }
+
for (int i1_V = 0; i1_V < n1_V; i1_V++) {
L1TrackPar& T = T_1[i1_V];
L1FieldRegion& fld1 =
@@ -308,23 +317,12 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
for (int ista = 0; ista <= istal - 1; ista++) {
- if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom);
+ fit.L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom, 1);
#else
- L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom);
+ fit.L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom, 1);
#endif
- if (ista == NMvdStations - 1) L1AddPipeMaterial(T, MaxInvMom);
- }
- else {
-#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f);
-
-#else
- L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f);
-#endif
- if (ista == NMvdStations - 1) L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f);
- }
+ if (ista == NMvdStations - 1) fit.L1AddPipeMaterial(T, MaxInvMom, 1);
}
// add left hit
@@ -356,23 +354,13 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
FilterTime(T, time, timeEr);
- if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom);
+ fit.L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom, 1);
#else
- L1AddMaterial(T, stal.materialInfo, MaxInvMom);
+ fit.L1AddMaterial(T, stal.materialInfo, MaxInvMom, 1);
#endif
- if ((istam >= NMvdStations) && (istal <= NMvdStations - 1)) L1AddPipeMaterial(T, MaxInvMom);
- }
- else {
-#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f);
-#else
- L1AddMaterial(T, stal.materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f);
-#endif
- if ((istam >= NMvdStations) && (istal <= NMvdStations - 1))
- L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f);
- }
+ if ((istam >= NMvdStations) && (istal <= NMvdStations - 1)) { fit.L1AddPipeMaterial(T, MaxInvMom, 1); }
+
fvec dz = zstam - T.z;
L1ExtrapolateTime(T, dz);
@@ -572,6 +560,9 @@ inline void L1Algo::f30( // input
fvec fvec_0;
L1TrackPar L1TrackPar_0;
+ L1Fit fit;
+ fit.SetParticleMass(fDefaultMass);
+
Tindex n3_V = 0, n3_4 = 0;
T_3.push_back(L1TrackPar_0);
@@ -669,23 +660,12 @@ inline void L1Algo::f30( // input
FilterTime(T2, timeM, timeMEr);
- if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
-#ifdef USE_RL_TABLE
- L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp);
-#else
- L1AddMaterial(T2, stam.materialInfo, T2.qp);
-#endif
- if ((istar >= NMvdStations) && (istam <= NMvdStations - 1)) L1AddPipeMaterial(T2, T2.qp);
- }
- else {
#ifdef USE_RL_TABLE
- L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1, 0.000511f * 0.000511f);
+ fit.L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1);
#else
- L1AddMaterial(T2, stam.materialInfo, T2.qp, 1, 0.000511f * 0.000511f);
+ fit.L1AddMaterial(T2, stam.materialInfo, T2.qp, 1);
#endif
- if ((istar >= NMvdStations) && (istam <= NMvdStations - 1))
- L1AddPipeMaterial(T2, T2.qp, 1, 0.000511f * 0.000511f);
- }
+ if ((istar >= NMvdStations) && (istam <= NMvdStations - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); }
fvec dz2 = star.z - T2.z;
L1ExtrapolateTime(T2, dz2);
@@ -925,6 +905,9 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction
Tindex n3, int istal, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, vector<THitI>& hitsm_3,
vector<THitI>& hitsr_3, int nIterations)
{
+ L1Fit fit;
+ fit.SetParticleMass(fDefaultMass);
+
const int NHits = 3; // triplets
// prepare data
@@ -997,11 +980,11 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction
for (int ih = 1; ih < NHits; ++ih) {
L1Extrapolate(T, z[ih], T.qp, fld);
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
+ fit.L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1);
#else
- L1AddMaterial(T, sta[ih].materialInfo, T.qp);
+ fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1);
#endif
- if (ista[ih] == NMvdStations - 1) L1AddPipeMaterial(T, T.qp);
+ if (ista[ih] == NMvdStations - 1) fit.L1AddPipeMaterial(T, T.qp, 1);
L1Filter(T, sta[ih].frontInfo, u[ih]);
L1Filter(T, sta[ih].backInfo, v[ih]);
@@ -1032,11 +1015,11 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction
for (ih = NHits - 2; ih >= 0; ih--) {
L1Extrapolate(T, z[ih], T.qp, fld);
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
+ fit.L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1);
#else
- L1AddMaterial(T, sta[ih].materialInfo, T.qp);
+ fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1);
#endif
- if (ista[ih] == NMvdStations) L1AddPipeMaterial(T, T.qp);
+ if (ista[ih] == NMvdStations) fit.L1AddPipeMaterial(T, T.qp, 1);
L1Filter(T, sta[ih].frontInfo, u[ih]);
L1Filter(T, sta[ih].backInfo, v[ih]);
@@ -1063,11 +1046,11 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction
for (ih = 1; ih < NHits; ++ih) {
L1Extrapolate(T, z[ih], T.qp, fld);
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
+ fit.L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1);
#else
- L1AddMaterial(T, sta[ih].materialInfo, T.qp);
+ fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1);
#endif
- if (ista[ih] == NMvdStations + 1) L1AddPipeMaterial(T, T.qp);
+ if (ista[ih] == NMvdStations + 1) fit.L1AddPipeMaterial(T, T.qp, 1);
L1Filter(T, sta[ih].frontInfo, u[ih]);
L1Filter(T, sta[ih].backInfo, v[ih]);
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f7fcb723df6794f70d565aa3678f2effa97b26b3
--- /dev/null
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -0,0 +1,8 @@
+/// \file L1Fit.cxx
+/// \brief Implementation of the L1Fit class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \date 15.03.2021
+
+#include "L1Fit.h"
+
+//ClassImp(L1Fit);
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
new file mode 100644
index 0000000000000000000000000000000000000000..cd10e2e12ad05802c0ab111aa05b2f41c58d14e1
--- /dev/null
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -0,0 +1,82 @@
+/// \file L1Fit.h
+/// \brief Definition of the L1Fit class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \date 15.03.2021
+
+#ifndef L1Fit_H
+#define L1Fit_H
+
+#include "CbmL1Def.h"
+
+#include "L1MaterialInfo.h"
+#include "L1TrackPar.h"
+
+///
+/// A collection of fit routines used by the CA tracker
+///
+class L1Fit {
+
+public:
+ /// constructor
+ L1Fit() {};
+
+ /// destructor
+ ~L1Fit() {};
+
+ /// set particle mass for the fit
+ void SetParticleMass(fvec mass) { fMass2 = mass * mass; }
+
+ /// get the particle mass squared
+ fvec GetParticleMass2() const { return fMass2; }
+
+ /// === routines for material
+
+ fvec BetheBlochIron(const float qp);
+
+ fvec BetheBlochCarbon(const float qp);
+
+ fvec BetheBlochAl(const float qp);
+
+ static fvec ApproximateBetheBloch(const fvec& bg2);
+
+ static fvec ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
+ const fvec& kp4);
+
+
+ float CalcQpAfterEloss(float qp, float eloss, float mass2);
+
+ void EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
+
+ void EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
+
+
+ void EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
+
+
+ void EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w);
+
+ void L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w);
+
+ void L1AddMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0, fvec w);
+
+
+ void L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
+
+ void L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0);
+
+ void L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w);
+
+ void L1AddTargetMaterial(L1TrackPar& T, fvec qp0, fvec w);
+
+ const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium
+ const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold
+
+private:
+ fvec fMass2 = 0.10565800 * 0.10565800; // muon mass by default (pion 0.13957000f )
+
+ //ClassDefNV(L1Fit, 0)
+};
+
+#include "L1FitMaterial.h"
+
+#endif
diff --git a/reco/L1/L1Algo/L1AddMaterial.h b/reco/L1/L1Algo/L1FitMaterial.h
similarity index 63%
rename from reco/L1/L1Algo/L1AddMaterial.h
rename to reco/L1/L1Algo/L1FitMaterial.h
index 87aa9075f843e1b950cb7bb55322908b6e0adbad..2483c44802b741ea619214a6b31bd735e9863c69 100644
--- a/reco/L1/L1Algo/L1AddMaterial.h
+++ b/reco/L1/L1Algo/L1FitMaterial.h
@@ -1,18 +1,20 @@
-#ifndef L1AddMaterial_h
-#define L1AddMaterial_h
+#ifndef L1FitMaterial_h
+#define L1FitMaterial_h
#include "CbmL1Def.h"
+
#include "L1MaterialInfo.h"
#include "L1TrackPar.h"
//#define cnst static const fvec
#define cnst const fvec
-const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium
-const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold
+//const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium
+//const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold
-inline fvec BetheBlochIron(const float qp) {
+inline fvec L1Fit::BetheBlochIron(const float qp)
+{
float K = 0.000307075; // [GeV*cm^2/g]
float z = (qp > 0.) ? 1 : -1.;
@@ -31,8 +33,7 @@ inline fvec BetheBlochIron(const float qp) {
float me = 0.000511; // GeV
float ratio = me / M;
- float Tmax =
- (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
+ float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
// density correction
float dc = 0.;
@@ -42,12 +43,11 @@ inline fvec BetheBlochIron(const float qp) {
dc = log(hwp / I) + log(beta * gamma) - 0.5;
}
- return K * z * z * (Z / A) * (1. / betaSq)
- * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq
- - dc);
+ return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
}
-inline fvec BetheBlochCarbon(const float qp) {
+inline fvec L1Fit::BetheBlochCarbon(const float qp)
+{
float K = 0.000307075; // GeV * g^-1 * cm^2
float z = (qp > 0.) ? 1 : -1.;
@@ -66,8 +66,7 @@ inline fvec BetheBlochCarbon(const float qp) {
float me = 0.000511; // GeV
float ratio = me / M;
- float Tmax =
- (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
+ float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
// density correction
float dc = 0.;
@@ -77,12 +76,11 @@ inline fvec BetheBlochCarbon(const float qp) {
dc = log(hwp / I) + log(beta * gamma) - 0.5;
}
- return K * z * z * (Z / A) * (1. / betaSq)
- * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq
- - dc);
+ return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
}
-inline fvec BetheBlochAl(const float qp) {
+inline fvec L1Fit::BetheBlochAl(const float qp)
+{
float K = 0.000307075; // GeV * g^-1 * cm^2
float z = (qp > 0.) ? 1 : -1.;
@@ -101,8 +99,7 @@ inline fvec BetheBlochAl(const float qp) {
float me = 0.000511; // GeV
float ratio = me / M;
- float Tmax =
- (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
+ float Tmax = (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
// density correction
float dc = 0.;
@@ -112,12 +109,11 @@ inline fvec BetheBlochAl(const float qp) {
dc = log(hwp / I) + log(beta * gamma) - 0.5;
}
- return K * z * z * (Z / A) * (1. / betaSq)
- * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq
- - dc);
+ return K * z * z * (Z / A) * (1. / betaSq) * (0.5 * log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq - dc);
}
-inline fvec ApproximateBetheBloch(const fvec& bg2) {
+inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
+{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
//
@@ -157,20 +153,15 @@ inline fvec ApproximateBetheBloch(const fvec& bg2) {
d2 = fvec(init & (lhwI + x - 0.5f));
const fvec& r = (x1 - x) / (x1 - x0);
init = (x > x0) & (x1 > x);
- d2 = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r))
- + fvec((!init) & d2);
+ d2 = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r)) + fvec((!init) & d2);
return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2)
- - d2);
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
}
-inline fvec ApproximateBetheBloch(const fvec& bg2,
- const fvec& kp0,
- const fvec& kp1,
- const fvec& kp2,
- const fvec& kp3,
- const fvec& kp4) {
+inline fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2,
+ const fvec& kp3, const fvec& kp4)
+{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
//
@@ -210,41 +201,35 @@ inline fvec ApproximateBetheBloch(const fvec& bg2,
d2 = fvec(init & (lhwI + x - 0.5f));
const fvec& r = (x1 - x) / (x1 - x0);
init = (x > x0) & (x1 > x);
- d2 = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r))
- + fvec((!init) & d2);
+ d2 = fvec(init & (lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r)) + fvec((!init) & d2);
return mK * mZA * (fvec(1.f) + bg2) / bg2
- * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2)
- - d2);
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
}
-inline float CalcQpAfterEloss(float qp, float eloss, float mass2) {
+inline float L1Fit::CalcQpAfterEloss(float qp, float eloss, float mass2)
+{
float p = fabs(1. / qp);
float E = sqrt(p * p + mass2);
float q = (qp > 0) ? 1. : -1.;
float Enew = E + eloss;
- // float pnew = (Enew > sqrt(mass2)) ? std::sqrt(Enew * Enew - mass2) : 0.;
+ // float pnew = (Enew > sqrt(fMass2)) ? std::sqrt(Enew * Enew - fMass2) : 0.;
float pnew = std::sqrt(Enew * Enew - mass2);
// if (pnew!=0)
return q / pnew;
// else { return 1e5; }
}
-inline void EnergyLossCorrection(L1TrackPar& T,
- const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1,
- fvec /*dE_*/ = 0) {
+inline void L1Fit::EnergyLossCorrection(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (T.qp * T.qp);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
// fvec qp =T.qp;
- const fvec& bethe = ApproximateBetheBloch(p2 / mass2);
+ const fvec& bethe = ApproximateBetheBloch(p2 / fMass2);
// fvec bethe2 = BetheBlochAl(qp[0]);
@@ -257,15 +242,14 @@ inline void EnergyLossCorrection(L1TrackPar& T,
// dE = fabs(dE_);//dE2;//bethe * radThick*tr * 9.34961f*2.33f ;
- const fvec& E2Corrected =
- (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
- fvec init = fvec(!(corr == corr)) | fvec(w < 1);
- corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
+ const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
+ fvec init = fvec(!(corr == corr)) | fvec(w < 1);
+ corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
qp0 *= corr;
- // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], mass2[0]);
+ // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], fMass2[0]);
// qp0 = dqp;
// T.qp = dqp;
T.qp *= corr;
@@ -276,15 +260,10 @@ inline void EnergyLossCorrection(L1TrackPar& T,
T.C44 *= corr * corr;
}
-inline void EnergyLossCorrectionAl(L1TrackPar& T,
- const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1,
- fvec /*dE_*/ = 0) {
+inline void L1Fit::EnergyLossCorrectionAl(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (T.qp * T.qp);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
fvec qp = T.qp;
@@ -300,15 +279,13 @@ inline void EnergyLossCorrectionAl(L1TrackPar& T,
fvec i;
- if (atomicZ < 13)
- i = (12. * atomicZ + 7.) * 1.e-9;
+ if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
- // const fvec& bethe = ApproximateBetheBloch( p2/mass2 );
+ // const fvec& bethe = ApproximateBetheBloch( p2/fMass2 );
// fvec bethe2 = BetheBlochAl(qp[0]);
@@ -319,14 +296,13 @@ inline void EnergyLossCorrectionAl(L1TrackPar& T,
fvec dE = bethe * radThick * tr * radLen * rho;
- const fvec& E2Corrected =
- (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
- fvec init = fvec(!(corr == corr)) | fvec(w < 1);
- corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
+ const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
+ fvec init = fvec(!(corr == corr)) | fvec(w < 1);
+ corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
qp0 *= corr;
- // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], mass2[0]);
+ // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], fMass2[0]);
// qp0 = dqp;
// T.qp = dqp;
T.qp *= corr;
@@ -341,7 +317,7 @@ inline void EnergyLossCorrectionAl(L1TrackPar& T,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2);
- fvec GAMMA = sqrt(E2) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -349,7 +325,7 @@ inline void EnergyLossCorrectionAl(L1TrackPar& T,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
@@ -374,15 +350,10 @@ inline void EnergyLossCorrectionAl(L1TrackPar& T,
}
-inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
- const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1,
- fvec /*dE_*/ = 0) {
+inline void L1Fit::EnergyLossCorrectionCarbon(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (T.qp * T.qp);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
fvec qp = T.qp;
@@ -393,16 +364,14 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
float radLen = 18.76f;
fvec i;
- if (atomicZ < 13)
- i = (12. * atomicZ + 7.) * 1.e-9;
+ if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
- // const fvec& bethe = ApproximateBetheBloch( p2/mass2 );
+ // const fvec& bethe = ApproximateBetheBloch( p2/fMass2 );
// fvec bethe2 = BetheBlochCarbon(qp[0]);
@@ -413,11 +382,10 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
fvec dE = bethe * radThick * tr * rho * radLen;
// fvec dE2 = bethe2 * radThick*tr*2.70f* 18.76f;
- const fvec& E2Corrected =
- (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
- fvec init = fvec(!(corr == corr)) | fvec(w < 1);
- corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
+ const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
+ fvec init = fvec(!(corr == corr)) | fvec(w < 1);
+ corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
qp0 *= corr;
T.qp *= corr;
@@ -433,7 +401,7 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -441,7 +409,7 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
@@ -450,7 +418,7 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
fvec SDEDX = ((E2) *DEDX2) / std::pow(P, 6);
- // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], mass2[0]);
+ // fvec dqp = CalcQpAfterEloss(qp[0], (direction*dE)[0], fMass2[0]);
// qp0 = dqp;
// T.qp = dqp;
@@ -470,15 +438,10 @@ inline void EnergyLossCorrectionCarbon(L1TrackPar& T,
}
-inline void EnergyLossCorrectionIron(L1TrackPar& T,
- const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1,
- fvec /*dE_*/ = 0) {
+inline void L1Fit::EnergyLossCorrectionIron(L1TrackPar& T, const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (T.qp * T.qp);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
fvec qp = T.qp;
@@ -488,13 +451,11 @@ inline void EnergyLossCorrectionIron(L1TrackPar& T,
float radLen = 1.758f;
fvec i;
- if (atomicZ < 13)
- i = (12. * atomicZ + 7.) * 1.e-9;
+ if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
// fvec bethe2 = BetheBlochIron(T.qp[0]);
@@ -505,11 +466,10 @@ inline void EnergyLossCorrectionIron(L1TrackPar& T,
fvec dE = bethe * radThick * tr * radLen * rho;
//fvec dE2 = bethe2 * radThick*tr * 1.758f*2.33;
- const fvec& E2Corrected =
- (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
- fvec init = fvec(!(corr == corr)) | fvec(w < 1);
- corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
+ const fvec& E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
+ fvec init = fvec(!(corr == corr)) | fvec(w < 1);
+ corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
qp0 *= corr;
T.qp *= corr;
@@ -524,7 +484,7 @@ inline void EnergyLossCorrectionIron(L1TrackPar& T,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -532,7 +492,7 @@ inline void EnergyLossCorrectionIron(L1TrackPar& T,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
@@ -556,11 +516,8 @@ inline void EnergyLossCorrectionIron(L1TrackPar& T,
}
-inline void L1AddMaterial(L1TrackPar& T,
- fvec radThick,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.10565f * 0.10565f) {
+inline void L1Fit::L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w)
+{
cnst ONE = 1.;
fvec tx = T.tx;
@@ -573,24 +530,19 @@ inline void L1AddMaterial(L1TrackPar& T,
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 =
- (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
+ fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*radThick*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
T.C22 += w * txtx1 * a;
T.C32 += w * tx * ty * a;
T.C33 += w * (ONE + tyty) * a;
}
-inline void L1AddMaterial(L1TrackPar& T,
- L1MaterialInfo& info,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.10565f * 0.10565f) {
+inline void L1Fit::L1AddMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0, fvec w)
+{
cnst ONE = 1.f;
fvec tx = T.tx;
@@ -603,20 +555,18 @@ inline void L1AddMaterial(L1TrackPar& T,
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 =
- (c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
+ fvec s0 = (c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*info.RadThick*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
T.C22 += w * txtx1 * a;
T.C32 += w * tx * ty * a;
T.C33 += w * (ONE + tyty) * a;
}
-// inline void L1AddThickMaterial( L1TrackPar &T, fvec radThick, fvec qp0, fvec thickness=0, fvec w = 1, fvec mass2 = 0.10565f*0.10565f, bool fDownstream = 1 )
+// inline void L1Fit::L1AddThickMaterial( L1TrackPar &T, fvec radThick, fvec qp0, fvec thickness=0, fvec w = 1, fvec mass2 = 0.10565f*0.10565f, bool fDownstream = 1 )
// {
// cnst ONE = 1.;
//
@@ -658,14 +608,8 @@ inline void L1AddMaterial(L1TrackPar& T,
//
// }
-inline void L1AddThickMaterial(L1TrackPar& T,
- fvec radThick,
- fvec qp0,
- fvec w,
- fvec thickness = 0,
- bool fDownstream = 1) {
- fvec mass2 = 0.10565f * 0.10565f;
-
+inline void L1Fit::L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
+{
cnst ONE = 1.;
fvec tx = T.tx;
@@ -678,13 +622,11 @@ inline void L1AddThickMaterial(L1TrackPar& T,
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 =
- (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
+ fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*radThick*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
fvec D = (fDownstream) ? 1. : -1.;
fvec T23 = (thickness * thickness) / 3.0;
@@ -706,10 +648,9 @@ inline void L1AddThickMaterial(L1TrackPar& T,
}
-inline void L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0) {
- cnst ONE = 1.;
- cnst mass2 = 0.10565f * 0.10565f;
-
+inline void L1Fit::L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0)
+{
+ cnst ONE = 1.f;
fvec tx = T.tx;
fvec ty = T.ty;
fvec txtx = tx * tx;
@@ -720,24 +661,19 @@ inline void L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0) {
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 = (c1 + c2 * (info.logRadThick + log(0.5)) + c3 * h
- + h2 * (c4 + c5 * h + c6 * h2))
- * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*0.5*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * info.RadThick * 0.5 * s0 * s0);
+ fvec s0 = (c1 + c2 * (info.logRadThick + log(0.5)) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*info.RadThick*0.5*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * 0.5 * s0 * s0);
T.C22 += txtx1 * a;
T.C32 += tx * ty * a;
T.C33 += (ONE + tyty) * a;
}
-inline void L1AddPipeMaterial(L1TrackPar& T,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.10565f * 0.10565f) {
+inline void L1Fit::L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w)
+{
cnst ONE = 1.f;
// static const fscal RadThick=0.0009f;//0.5/18.76;
@@ -755,23 +691,18 @@ inline void L1AddPipeMaterial(L1TrackPar& T,
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 =
- (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2))
- * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*RadThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ fvec s0 = (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*RadThick*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
T.C22 += w * txtx1 * a;
T.C32 += w * tx * ty * a;
T.C33 += w * (ONE + tyty) * a;
}
-inline void L1AddTargetMaterial(L1TrackPar& T,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.10565f * 0.10565f) {
+inline void L1Fit::L1AddTargetMaterial(L1TrackPar& T, fvec qp0, fvec w)
+{
cnst ONE = 1.f;
const fscal logRadThick = log(TargetRadThick[0]);
@@ -785,13 +716,10 @@ inline void L1AddTargetMaterial(L1TrackPar& T,
fvec h2 = h * h;
fvec qp0t = qp0 * t;
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
- c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
- fvec s0 =
- (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2))
- * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*RadThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * TargetRadThick * s0 * s0);
+ cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
+ fvec s0 = (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
+ //fvec a = ( (ONE+fMass2*qp0*qp0t)*RadThick*s0*s0 );
+ fvec a = ((t + fMass2 * qp0 * qp0t) * TargetRadThick * s0 * s0);
T.C22 += w * txtx1 * a;
T.C32 += w * tx * ty * a;
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index 8a31f5310d41570e68607b917e29a43b8cc83e74..6a7986ac370f908d40bc645464d9e1b5ec165317 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -1,15 +1,15 @@
-#include "L1AddMaterial.h"
+#include <iostream>
+
#include "L1Algo.h"
#include "L1Branch.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
+#include "L1Fit.h"
#include "L1HitArea.h"
#include "L1HitPoint.h"
#include "L1Track.h"
#include "L1TrackPar.h"
-#include <iostream>
-
// using namespace std;
using std::cout;
using std::endl;
@@ -29,6 +29,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t,
const bool initParams) {
L1_assert(t.NHits >= 3);
+ L1Fit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
+
// get hits of current track
const std::vector<THitI>& hits =
t.StsHits; // array of indeses of hits of current track
@@ -138,11 +141,11 @@ void L1Algo::BranchFitterFast(const L1Branch& t,
#endif
L1ExtrapolateTime(T, dz);
- L1AddMaterial(T, sta.materialInfo, qp0);
+ fit.L1AddMaterial(T, sta.materialInfo, qp0, 1);
if ((step * ista <= step * (NMvdStations + (step + 1) / 2 - 1))
&& (step * ista_prev
>= step * (NMvdStations + (step + 1) / 2 - 1 - step)))
- L1AddPipeMaterial(T, qp0);
+ fit.L1AddPipeMaterial(T, qp0, 1);
fvec u = hit.u;
fvec v = hit.v;
@@ -204,6 +207,8 @@ void L1Algo::FindMoreHits(L1Branch& t,
std::vector<THitI> newHits;
newHits.clear();
newHits.reserve(5);
+ L1Fit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
const signed short int step =
-2 * static_cast<int>(dir) + 1; // increment for station index
@@ -338,7 +343,7 @@ void L1Algo::FindMoreHits(L1Branch& t,
L1ExtrapolateTime(T, dz1);
L1ExtrapolateLine(T, z);
- L1AddMaterial(T, sta.materialInfo, qp0);
+ fit.L1AddMaterial(T, sta.materialInfo, qp0, 1);
L1UMeasurementInfo info = sta.frontInfo;
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 4f002816b8a0dc1fc3c7e7543965a1dab932d3d4..8298f692fcd0fd7d78b453735ca479d0c4ff36f1 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -12,16 +12,16 @@
*/
-#include "L1AddMaterial.h"
+#include <iostream>
+#include <vector>
+
#include "L1Algo.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h" // for KFTrackFitter_simple
+#include "L1Fit.h"
#include "L1TrackPar.h"
#include "L1TrackParFit.h"
-#include <iostream>
-#include <vector>
-
using std::cout;
using std::endl;
using std::vector;
@@ -55,6 +55,9 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
L1TrackPar T; // fitting parametr coresponding to current track
+ L1Fit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
+
fvec qp0 = 0.25;
//fvec qp0 = 2./t.Momentum;
for (int iter = 0; iter < 3; iter++) {
@@ -144,12 +147,12 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
// T.L1Extrapolate( sta.z, qp0, fld );
// L1Extrapolate( T, hit.z, qp0, fld );
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, ONE);
+ fit.L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, ONE);
#else
- L1AddMaterial(T, sta.materialInfo, qp0, ONE);
+ fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE);
#endif
- // if (ista==NMvdStations-1) L1AddPipeMaterial( T, qp0);
+ // if (ista==NMvdStations-1) fit.L1AddPipeMaterial( T, qp0, 1);
fvec u = hit.u;
fvec v = hit.v;
@@ -283,11 +286,11 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
// T.L1Extrapolate( sta.z, qp0, fld );
// L1Extrapolate( T, hit.z, qp0, fld );
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, ONE);
+ fit.L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, ONE);
#else
- L1AddMaterial(T, sta.materialInfo, qp0, ONE);
+ fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE);
#endif
- // if (ista==NMvdStations) L1AddPipeMaterial( T, qp0);
+ // if (ista==NMvdStations) fit.L1AddPipeMaterial( T, qp0, 1);
L1Filter(T, sta.frontInfo, u);
L1Filter(T, sta.backInfo, v);
@@ -351,6 +354,10 @@ void L1Algo::L1KFTrackFitter() {
L1TrackPar T; // fitting parametr coresponding to current track
L1TrackParFit T1; // fitting parametr coresponding to current track
+ T1.SetParticleMass(GetDefaultParticleMass());
+
+ L1Fit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
L1Track* t[fvecLen];
@@ -374,7 +381,6 @@ void L1Algo::L1KFTrackFitter() {
}
unsigned short N_vTracks = NTracksIsecAll;
- const fvec mass2 = 0.1395679f * 0.1395679f;
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
@@ -488,7 +494,7 @@ void L1Algo::L1KFTrackFitter() {
T1.Filter(time[i], timeEr[i], w_time[i]);
- // L1AddMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddMaterial( T, sta[i].materialInfo, qp0, 1 );
fz1 = z[i];
@@ -525,23 +531,21 @@ void L1Algo::L1KFTrackFitter() {
if (i == NMvdStations - 1) {
- L1AddPipeMaterial(T, qp0, wIn);
- EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(1.f), wIn);
+ fit.L1AddPipeMaterial(T, qp0, wIn);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
T1.L1AddPipeMaterial(qp01, wIn);
- T1.EnergyLossCorrection(mass2, PipeRadThick, qp01, fvec(1.f), wIn);
+ T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(1.f), wIn);
}
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
- EnergyLossCorrection(
- T, mass2, fRadThick[i].GetRadThick(T.x, T.y), qp0, fvec(1.f), wIn);
+ fit.L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, fvec(1.f), wIn);
T1.L1AddMaterial(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, wIn);
- T1.EnergyLossCorrection(
- mass2, fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(1.f), wIn);
+ T1.EnergyLossCorrection(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(1.f), wIn);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
T1.L1AddMaterial(sta[i].materialInfo, qp01, wIn);
#endif
L1UMeasurementInfo info = sta[i].frontInfo;
@@ -564,7 +568,7 @@ void L1Algo::L1KFTrackFitter() {
fB1 = fB0;
fz1 = fz0;
}
- // L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->TFirst[0] = T1.fx[iVec];
@@ -660,7 +664,7 @@ void L1Algo::L1KFTrackFitter() {
T1.Filter(time[i], timeEr[i], w_time[i]);
- // L1AddMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddMaterial( T, sta[i].materialInfo, qp0, 1 );
qp0 = T.qp;
qp01 = T1.fqp;
@@ -695,22 +699,20 @@ void L1Algo::L1KFTrackFitter() {
// T1.ExtrapolateLine( z[i]);
if (i == NMvdStations) {
- L1AddPipeMaterial(T, qp0, wIn);
- EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(-1.f), wIn);
+ fit.L1AddPipeMaterial(T, qp0, wIn);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
T1.L1AddPipeMaterial(qp01, wIn);
- T1.EnergyLossCorrection(mass2, PipeRadThick, qp01, fvec(-1.f), wIn);
+ T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(-1.f), wIn);
}
#ifdef USE_RL_TABLE
- L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
- EnergyLossCorrection(
- T, mass2, fRadThick[i].GetRadThick(T.x, T.y), qp0, fvec(-1.f), wIn);
+ fit.L1AddMaterial(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, fRadThick[i].GetRadThick(T.x, T.y), qp0, fvec(-1.f), wIn);
T1.L1AddMaterial(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, wIn);
- T1.EnergyLossCorrection(
- mass2, fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(-1.f), wIn);
+ T1.EnergyLossCorrection(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(-1.f), wIn);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
#endif
L1Filter(T, sta[i].frontInfo, u[i], w1);
L1Filter(T, sta[i].backInfo, v[i], w1);
@@ -732,7 +734,7 @@ void L1Algo::L1KFTrackFitter() {
fB1 = fB0;
fz1 = fz0;
}
- // L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->TLast[0] = T1.fx[iVec];
@@ -791,6 +793,10 @@ void L1Algo::L1KFTrackFitterMuch() {
L1TrackPar T; // fitting parametr coresponding to current track
L1TrackParFit T1; // fitting parametr coresponding to current track
+ T1.SetParticleMass(GetDefaultParticleMass());
+
+ L1Fit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
L1Track* t[fvecLen];
@@ -814,7 +820,6 @@ void L1Algo::L1KFTrackFitterMuch() {
}
unsigned short N_vTracks = NTracksIsecAll;
- const fvec mass2 = 0.10565f * 0.10565f;
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
@@ -979,11 +984,11 @@ void L1Algo::L1KFTrackFitterMuch() {
T1.Extrapolate(z[i], qp01, fld, &w1);
if (i == NMvdStations) {
- L1AddPipeMaterial(T, qp0, wIn);
- EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(-1.f), wIn);
+ fit.L1AddPipeMaterial(T, qp0, wIn);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
T1.L1AddPipeMaterial(qp01, wIn);
- T1.EnergyLossCorrection(mass2, PipeRadThick, qp01, fvec(-1.f), wIn);
+ T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(-1.f), wIn);
}
fB2 = fB1;
@@ -991,11 +996,7 @@ void L1Algo::L1KFTrackFitterMuch() {
fB1 = fB0;
fz1 = fz0;
#ifdef USE_RL_TABLE
- T1.EnergyLossCorrection(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy),
- qp01,
- fvec(-1.f),
- wIn);
+ T1.EnergyLossCorrection(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(-1.f), wIn);
T1.L1AddThickMaterial(fRadThick[i].GetRadThick(T1.fx, T1.fy),
qp01,
@@ -1073,27 +1074,17 @@ void L1Algo::L1KFTrackFitterMuch() {
// L1ExtrapolateLine( T, z_last);
#ifdef USE_RL_TABLE
if (i == 11 || i == 14 || i == 17)
- T1.EnergyLossCorrectionIron(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy)
- / (nofSteps + 1),
- qp01,
- fvec(-1.f),
+ T1.EnergyLossCorrectionIron(fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1), qp01, fvec(-1.f),
wIn);
if (i == 8)
T1.EnergyLossCorrectionCarbon(
- mass2,
fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1),
qp01,
fvec(-1.f),
wIn);
if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16
|| i == 18 || i == 19)
- T1.EnergyLossCorrectionAl(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy)
- / (nofSteps + 1),
- qp01,
- fvec(-1.f),
- wIn);
+ T1.EnergyLossCorrectionAl(fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1), qp01, fvec(-1.f), wIn);
T1.L1AddThickMaterial(
fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + fvec(1)),
@@ -1105,7 +1096,7 @@ void L1Algo::L1KFTrackFitterMuch() {
wIn = wIn & (one - mask1);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
#endif
}
@@ -1124,7 +1115,7 @@ void L1Algo::L1KFTrackFitterMuch() {
T1.Filter(time[i], timeEr[i], w1);
}
}
- // L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
@@ -1187,7 +1178,7 @@ void L1Algo::L1KFTrackFitterMuch() {
// fvec initialised = fvec(z[i] < z_end) & fvec(z_start <= z[i]);
// fvec wIn = (ONE & (initialised));
- // T1.EnergyLossCorrectionAl( mass2, fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(1.f), wIn);
+ // T1.EnergyLossCorrectionAl( fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(1.f), wIn);
//
// T1.L1AddThickMaterial(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, ONE, sta[i].materialInfo.thick, 0);
@@ -1240,27 +1231,16 @@ void L1Algo::L1KFTrackFitterMuch() {
#ifdef USE_RL_TABLE
if (i == 11 || i == 14 || i == 17)
- T1.EnergyLossCorrectionIron(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy)
- / (nofSteps + 1),
- qp01,
- fvec(1.f),
- w2);
+ T1.EnergyLossCorrectionIron(fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1), qp01, fvec(1.f), w2);
if (i == 8)
T1.EnergyLossCorrectionCarbon(
- mass2,
fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1),
qp01,
fvec(1.f),
w2);
if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16
|| i == 18)
- T1.EnergyLossCorrectionAl(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy)
- / (nofSteps + 1),
- qp01,
- fvec(1.f),
- w2);
+ T1.EnergyLossCorrectionAl(fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + 1), qp01, fvec(1.f), w2);
T1.L1AddThickMaterial(
fRadThick[i].GetRadThick(T1.fx, T1.fy) / (nofSteps + fvec(1)),
@@ -1272,7 +1252,7 @@ void L1Algo::L1KFTrackFitterMuch() {
w2 = w2 & (one - mask1);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, w2);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, w2);
#endif
}
@@ -1333,11 +1313,7 @@ void L1Algo::L1KFTrackFitterMuch() {
#ifdef USE_RL_TABLE
- T1.EnergyLossCorrection(mass2,
- fRadThick[i].GetRadThick(T1.fx, T1.fy),
- qp01,
- fvec(1.f),
- wIn);
+ T1.EnergyLossCorrection(fRadThick[i].GetRadThick(T1.fx, T1.fy), qp01, fvec(1.f), wIn);
T1.L1AddThickMaterial(fRadThick[i].GetRadThick(T1.fx, T1.fy),
qp01,
wIn,
@@ -1345,7 +1321,7 @@ void L1Algo::L1KFTrackFitterMuch() {
0);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
#endif
L1UMeasurementInfo info = sta[i].frontInfo;
@@ -1365,7 +1341,7 @@ void L1Algo::L1KFTrackFitterMuch() {
fz1 = fz0;
}
}
- // L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 6d8ff30623684966201fc4e12daff3456a3ec96b..eb00bb70cc0dc829b9c0b16b22a446e53ef72f5e 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -1,5 +1,10 @@
#include "L1TrackParFit.h"
+#include "L1Fit.h"
+
+const fvec PipeRadThick = 7.87e-3f; // 0.7 mm Aluminium
+const fvec TargetRadThick = 3.73e-2f * 2; // 250 mum Gold
+
#define cnst const fvec
@@ -430,7 +435,7 @@ void L1TrackParFit::
+ (-tx * B[step][1] + 2.f * ty * B[step][0]) * tx2ty2qp;
fvec p2 = 1.f / (qp0 * qp0);
- fvec m2 = 0.1395679f * 0.1395679f;
+ fvec m2 = fMass2;
fvec v = 29.9792458f * sqrt(p2 / (m2 + p2));
At[step] = sqrt(1.f + tx * tx + ty * ty) / v;
At_tx[step] = tx / sqrt(1.f + tx * tx + ty * ty) / v;
@@ -761,7 +766,8 @@ void L1TrackParFit::
+ ((!initialised) & C55);
}
-void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w, fvec mass2) {
+void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w)
+{
cnst ONE = 1.f;
// static const fscal RadThick=0.0009f;//0.5/18.76;
@@ -785,7 +791,7 @@ void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w, fvec mass2) {
(c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2))
* qp0t;
//fvec a = ( (ONE+mass2*qp0*qp0t)*RadThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
+ fvec a = ((t + fMass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
C22 += w * txtx1 * a;
C32 += w * tx * ty * a;
@@ -793,7 +799,8 @@ void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w, fvec mass2) {
}
-void L1TrackParFit::L1AddMaterial(fvec radThick, fvec qp0, fvec w, fvec mass2) {
+void L1TrackParFit::L1AddMaterial(fvec radThick, fvec qp0, fvec w)
+{
cnst ONE = 1.;
fvec tx = ftx;
@@ -812,7 +819,7 @@ void L1TrackParFit::L1AddMaterial(fvec radThick, fvec qp0, fvec w, fvec mass2) {
fvec s0 =
(c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
//fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
+ fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
C22 += w * txtx1 * a;
C32 += w * tx * ty * a;
@@ -822,7 +829,6 @@ void L1TrackParFit::L1AddMaterial(fvec radThick, fvec qp0, fvec w, fvec mass2) {
void L1TrackParFit::L1AddThickMaterial(fvec radThick,
fvec qp0,
fvec w,
- fvec mass2,
fvec thickness,
bool fDownstream) {
cnst ONE = 1.;
@@ -843,7 +849,7 @@ void L1TrackParFit::L1AddThickMaterial(fvec radThick,
fvec s0 =
(c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
//fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
+ fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
fvec D = (fDownstream) ? 1. : -1.;
fvec T23 = (thickness * thickness) / 3.0;
@@ -864,10 +870,8 @@ void L1TrackParFit::L1AddThickMaterial(fvec radThick,
}
-void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info,
- fvec qp0,
- fvec w,
- fvec mass2) {
+void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info, fvec qp0, fvec w)
+{
cnst ONE = 1.f;
fvec tx = ftx;
@@ -887,7 +891,7 @@ void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info,
fvec s0 =
(c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
//fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*s0*s0 );
- fvec a = ((t + mass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
+ fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
C22 += w * txtx1 * a;
C32 += w * tx * ty * a;
@@ -895,15 +899,12 @@ void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info,
}
-void L1TrackParFit::EnergyLossCorrection(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w) {
+void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (fqp * fqp);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
- const fvec& bethe = ApproximateBetheBloch(p2 / mass2);
+ const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2);
fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
@@ -911,7 +912,7 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& mass2,
const fvec& E2Corrected =
(sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
fvec init = fvec(!(corr == corr)) | fvec(w < 1);
corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -925,13 +926,10 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& mass2,
C44 *= corr * corr;
}
-void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w) {
+void L1TrackParFit::EnergyLossCorrectionIron(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (qp0 * qp0);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
int atomicZ = 26;
float atomicA = 55.845f;
@@ -944,8 +942,7 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
@@ -954,7 +951,7 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
const fvec& E2Corrected =
(sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
fvec init = fvec(!(corr == corr)) | fvec(w < 1);
corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -971,7 +968,7 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2Corrected) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -979,7 +976,7 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
@@ -996,13 +993,10 @@ void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
C44 += fabs(SDEDX);
}
-void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w) {
+void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (qp0 * qp0);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
int atomicZ = 6;
float atomicA = 12.011f;
@@ -1015,8 +1009,7 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
@@ -1025,7 +1018,7 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
const fvec& E2Corrected =
(sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
fvec init = fvec(!(corr == corr)) | fvec(w < 1);
corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -1042,7 +1035,7 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2Corrected) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -1050,7 +1043,7 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
@@ -1067,13 +1060,10 @@ void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
C44 += fabs(SDEDX);
}
-void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w) {
+void L1TrackParFit::EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
+{
const fvec& p2 = 1.f / (qp0 * qp0);
- const fvec& E2 = mass2 + p2;
+ const fvec& E2 = fMass2 + p2;
int atomicZ = 13;
float atomicA = 26.981f;
@@ -1086,8 +1076,7 @@ void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec& bethe =
- ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec& bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
@@ -1096,7 +1085,7 @@ void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
const fvec& E2Corrected =
(sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - mass2));
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
fvec init = fvec(!(corr == corr)) | fvec(w < 1);
corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
@@ -1113,7 +1102,7 @@ void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
fvec EMASS = 0.511 * 1e-3; // GeV
fvec BETA = P / sqrt(E2Corrected);
- fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
+ fvec GAMMA = sqrt(E2Corrected) / sqrt(fMass2);
// Calculate xi factor (KeV).
fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
@@ -1121,7 +1110,7 @@ void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
// Maximum energy transfer to atomic electron (KeV).
fvec ETA = BETA * GAMMA;
fvec ETASQ = ETA * ETA;
- fvec RATIO = EMASS / sqrt(mass2);
+ fvec RATIO = EMASS / sqrt(fMass2);
fvec F1 = 2. * EMASS * ETASQ;
fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
fvec EMAX = 1e6 * F1 / F2;
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 9118ac6d2bda5479f73facd62af72cc83b6c9841..ee50616d00a94e2edefd082ca94715833149f036 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -1,13 +1,11 @@
#ifndef L1TrackParFit_h
#define L1TrackParFit_h
-#include "L1AddMaterial.h"
+#include "../CbmL1Def.h"
#include "L1Field.h"
-#include "L1UMeasurementInfo.h"
-
+#include "L1MaterialInfo.h"
#include "L1TrackPar.h"
-
-#include "../CbmL1Def.h"
+#include "L1UMeasurementInfo.h"
class L1TrackParFit {
@@ -16,6 +14,7 @@ public:
C32, C33, C40, C41, C42, C43, C44, C50, C51, C52, C53, C54, C55, chi2, NDF;
// fvec n;
+ fvec fMass2 = 0.000511f * 0.000511f; // muon mass
L1TrackParFit()
: fx(0)
@@ -87,6 +86,13 @@ public:
void Print(int i = -1);
//Fit functionality
+
+ /// set particle mass for the fit
+ void SetParticleMass(float mass) { fMass2 = mass * mass; }
+
+ /// get the particle mass squared
+ fvec GetParticleMass2() const { return fMass2; }
+
void Filter(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
void Filter(fvec t0, fvec dt0, fvec w = 1.);
void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
@@ -94,42 +100,16 @@ public:
void ExtrapolateLine(fvec z_out, fvec* w = 0);
void ExtrapolateLine1(fvec z_out, fvec* w = 0, fvec v = 0);
void Compare(L1TrackPar& T);
- void EnergyLossCorrection(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1);
- void L1AddMaterial(L1MaterialInfo& info,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.1395679f * 0.1395679f);
- void L1AddMaterial(fvec radThick,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.1395679f * 0.1395679f);
- void L1AddThickMaterial(fvec radThick,
- fvec qp0,
- fvec w = 1,
- fvec mass2 = 0.1395679f * 0.1395679f,
- fvec thickness = 0,
- bool fDownstream = 1);
- void
- L1AddPipeMaterial(fvec qp0, fvec w = 1, fvec mass2 = 0.1395679f * 0.1395679f);
- void EnergyLossCorrectionIron(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1);
- void EnergyLossCorrectionCarbon(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1);
- void EnergyLossCorrectionAl(const fvec& mass2,
- const fvec& radThick,
- fvec& qp0,
- fvec direction,
- fvec w = 1);
+ void EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
+ void L1AddMaterial(L1MaterialInfo& info, fvec qp0, fvec w = 1);
+
+ void L1AddMaterial(fvec radThick, fvec qp0, fvec w = 1);
+
+ void L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
+ void L1AddPipeMaterial(fvec qp0, fvec w = 1);
+ void EnergyLossCorrectionIron(const fvec& radThick, fvec& qp0, fvec direction, fvec w = 1);
+ void EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w = 1);
+ void EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w = 1);
// void L1Extrapolate
// (
// // L1TrackParFit &T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
diff --git a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
index e193a079ada319e9ff8153402229e7457a434ed6..579a94ae7435fd847d59ca2a16fcd16932e12cbd 100644
--- a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
+++ b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
@@ -15,21 +15,23 @@
*/
#include "CbmL1StsTrackFinder.h"
-#include "L1Algo/L1Algo.h"
-
#include "CbmEvent.h"
#include "CbmKFMath.h"
#include "CbmStsHit.h"
#include "CbmStsTrack.h"
+
#include "FairHit.h"
#include "FairMCPoint.h"
#include "FairRootManager.h"
#include "TClonesArray.h"
+#include "TDatabasePDG.h"
#include <iostream>
#include <vector>
+#include "L1Algo/L1Algo.h"
+
using std::cout;
using std::endl;
using std::vector;
@@ -136,3 +138,23 @@ Int_t CbmL1StsTrackFinder::FindTracks(CbmEvent* event) {
return nTracks;
}
// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void SetDefaultParticlePDG(int pdg = 211)
+{
+ /// set a default particle mass for the track fit
+ /// it is used during reconstruction for the multiple scattering estimation
+ CbmL1* l1 = CbmL1::Instance();
+ if (!l1 || !l1->algo) {
+ LOG(fatal) << "L1 instance doesn't exist or is not initialised";
+ return;
+ }
+ auto* p = TDatabasePDG::Instance()->GetParticle(pdg);
+ if (!p) {
+ LOG(fatal) << "Particle with pdg " << pdg << " doesn't exist";
+ return;
+ }
+ l1->algo->SetDefaultParticleMass(p->Mass());
+}
+// -------------------------------------------------------------------------
diff --git a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.h b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.h
index 221097e65a2058c97014a3d13a3f1c358af6ba0b..38e2f8d4c7202715cc2c45786a371e600ecbb918 100644
--- a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.h
+++ b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.h
@@ -51,6 +51,10 @@ public:
**/
virtual Int_t FindTracks(CbmEvent* event);
+ /// set a default particle mass for the track fit
+ /// it is used during reconstruction
+ /// for the multiple scattering and energy loss estimation
+ void SetDefaultParticlePDG(int pdg = 13); // muon
private:
/** Copy the tracks from the L1-internal format and array
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 118780a367dad36b3ea900106845bebb4baf17b7..90b73a83d9656b2dff5d6b1767252c08a085bbfa 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -22,21 +22,22 @@
#include "TClonesArray.h"
//L1Algo tools
+#include "CbmKFVertex.h"
#include "CbmL1Track.h"
-#include "L1AddMaterial.h"
+
+#include "FairRootManager.h"
+
+#include "TDatabasePDG.h"
+
+#include "KFParticleDatabase.h"
#include "L1Algo.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
+#include "L1Fit.h"
#include "L1MaterialInfo.h"
#include "L1Station.h"
#include "L1TrackPar.h"
-#include "FairRootManager.h"
-#include "TDatabasePDG.h"
-
-#include "CbmKFVertex.h"
-#include "KFParticleDatabase.h"
-
using std::vector;
@@ -96,6 +97,9 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
int nTracks_SIMD = fvecLen;
L1TrackPar T; // fitting parametr coresponding to current track
+ L1Fit fit;
+ fit.SetParticleMass(0.000511f); // muon
+
CbmStsTrack* t[fvecLen];
int ista;
@@ -120,7 +124,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
Tracks[itrack].SetPidHypo(pidHypo[itrack]);
}
- fvec mass, mass2;
+ fvec mass = 0.000511f;
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
@@ -155,7 +159,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
// mass[i] = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
mass[i] = KFParticleDatabase::Instance()->GetMass(pid);
}
- mass2 = mass * mass;
+ fit.SetParticleMass(mass);
+
// get hits of current track
for (i = 0; i < nHits; i++) {
w[i] = ZERO;
@@ -247,24 +252,14 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
L1Extrapolate(T, z[i], qp0, fld, &w1);
if (i == NMvdStations) {
- L1AddPipeMaterial(T, qp0, wIn);
- EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(-1.f), wIn);
+ fit.L1AddPipeMaterial(T, qp0, wIn);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
}
#ifdef USE_RL_TABLE
- L1AddMaterial(T,
- CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
- qp0,
- wIn,
- mass2);
- EnergyLossCorrection(
- T,
- mass2,
- CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
- qp0,
- -1,
- wIn);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, -1, wIn);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn, mass2);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
#endif
L1Filter(T, sta[i].frontInfo, u[i], w1);
L1Filter(T, sta[i].backInfo, v[i], w1);
@@ -332,24 +327,14 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
L1Extrapolate(T, z[i], qp0, fld, &w1);
if (i == NMvdStations - 1) {
- L1AddPipeMaterial(T, qp0, wIn);
- EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(1.f), wIn);
+ fit.L1AddPipeMaterial(T, qp0, wIn);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
}
#ifdef USE_RL_TABLE
- L1AddMaterial(T,
- CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
- qp0,
- wIn,
- mass2);
- EnergyLossCorrection(
- T,
- mass2,
- CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
- qp0,
- 1,
- wIn);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, 1, wIn);
#else
- L1AddMaterial(T, sta[i].materialInfo, qp0, wIn, mass2);
+ fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
#endif
L1Filter(T, sta[i].frontInfo, u[i], w1);
L1Filter(T, sta[i].backInfo, v[i], w1);
@@ -409,6 +394,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks,
int nTracks_SIMD = fvecLen;
L1TrackPar T; // fitting parametr coresponding to current track
+ L1Fit fit;
CbmStsTrack* t[fvecLen];
@@ -513,29 +499,18 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks,
L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
if (iSt == NMvdStations - 1) {
- L1AddPipeMaterial(T, T.qp, w, mass2);
- EnergyLossCorrection(T, mass2, PipeRadThick, T.qp, fvec(1.f), w);
+ fit.L1AddPipeMaterial(T, T.qp, w);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, T.qp, fvec(1.f), w);
}
- L1AddMaterial(
- T,
- CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y),
- T.qp,
- w,
- mass2);
- EnergyLossCorrection(
- T,
- mass2,
- CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y),
- T.qp,
- fvec(1.f),
- w);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, w);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, fvec(1.f), w);
}
if (NMvdStations <= 0) {
- L1AddPipeMaterial(T, T.qp, ONE, mass2);
- EnergyLossCorrection(T, mass2, PipeRadThick, T.qp, fvec(1.f), ONE);
+ fit.L1AddPipeMaterial(T, T.qp, ONE);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, T.qp, fvec(1.f), ONE);
}
L1Extrapolate(T, primVtx.GetRefZ(), T.qp, fld);
- L1AddTargetMaterial(T, T.qp);
+ fit.L1AddTargetMaterial(T, T.qp, 1);
Double_t Cv[3] = {primVtx.GetCovMatrix()[0],
primVtx.GetCovMatrix()[1],