diff --git a/reco/KF/CbmKF.cxx b/reco/KF/CbmKF.cxx index 101bb865a9e0e42c8922ee58af8510e4954fd0c9..5f95e6e42e408ca153f4cd5939c4c32c66af9783 100644 --- a/reco/KF/CbmKF.cxx +++ b/reco/KF/CbmKF.cxx @@ -76,7 +76,9 @@ CbmKF::CbmKF(const char* name, Int_t iVerbose) , fMethod(1) , fMaterialID2IndexMap() { - if (!fInstance) fInstance = this; + if (!fInstance) { + fInstance = this; + } } CbmKF::~CbmKF() { fInstance = nullptr; } @@ -126,12 +128,19 @@ InitStatus CbmKF::Init() FairRunAna* Run = FairRunAna::Instance(); //FairRuntimeDb *RunDB = Run->GetRuntimeDb(); (VF) not used - if (fVerbose) cout << "KALMAN FILTER : === INIT MAGNETIC FIELD ===" << endl; + if (fVerbose) { + cout << "KALMAN FILTER : === INIT MAGNETIC FIELD ===" << endl; + } fMagneticField = reinterpret_cast<FairField*>(Run->GetField()); - if (fVerbose && fMagneticField) cout << "Magnetic field done" << endl; - if (!fMagneticField) cout << "No Magnetic Field Found" << endl; + if (fVerbose && fMagneticField) { + cout << "Magnetic field done" << endl; + } + + if (!fMagneticField) { + cout << "No Magnetic Field Found" << endl; + } /** */ @@ -145,7 +154,9 @@ InitStatus CbmKF::Init() if (useMVD) { auto mvdInterface = CbmMvdTrackingInterface::Instance(); - if (fVerbose) cout << "KALMAN FILTER : === READ MVD MATERIAL ===" << endl; + if (fVerbose) { + cout << "KALMAN FILTER : === READ MVD MATERIAL ===" << endl; + } int NStations = mvdInterface->GetNtrackingStations(); @@ -170,17 +181,20 @@ InitStatus CbmKF::Init() vMvdMaterial.push_back(tube); MvdStationIDMap.insert(pair<Int_t, Int_t>(tube.ID, ist)); - if (fVerbose) + if (fVerbose) { cout << " Mvd material ( id, z, dz, r, R, RadL, dz/RadL )= ( " << tube.ID << ", " << tube.z << ", " << tube.dz << ", " << tube.r << ", " << tube.R << ", " << tube.RadLength << ", " << tube.dz / tube.RadLength << " )" << endl; + } } } //=== Sts stations === - if (fVerbose) cout << "KALMAN FILTER : === READ STS MATERIAL ===" << endl; + if (fVerbose) { + cout << "KALMAN FILTER : === READ STS MATERIAL ===" << endl; + } auto stsInterface = CbmStsTrackingInterface::Instance(); int NStations = stsInterface->GetNtrackingStations(); @@ -205,10 +219,11 @@ InitStatus CbmKF::Init() vStsMaterial.push_back(tube); StsStationIDMap.insert(pair<Int_t, Int_t>(tube.ID, ist)); - if (fVerbose) + if (fVerbose) { cout << " Sts material ( id, z, dz, r, R, RadL, dz/RadL )= ( " << tube.ID << ", " << tube.z << ", " << tube.dz << ", " << tube.r << ", " << tube.R << ", " << tube.RadLength << ", " << tube.dz / tube.RadLength << " )" << endl; + } } // FU 05.03.2020 @@ -357,7 +372,9 @@ Int_t CbmKF::GetMaterialIndex(Int_t uid) Int_t CbmKF::ReadTube(CbmKFTube& tube, FairGeoNode* node) { - if (!node) return 1; + if (!node) { + return 1; + } TString name = node->getName(); TString Sname = node->getShapePointer()->GetName(); @@ -418,10 +435,18 @@ Int_t CbmKF::ReadTube(CbmKFTube& tube, FairGeoNode* node) double z1 = P->At(3 + i * 3 + 0); double r1 = P->At(3 + i * 3 + 1); double R1 = P->At(3 + i * 3 + 2); - if (r1 < r) r = r1; - if (R1 > R) R = R1; - if (z1 < z) z = z1; - if (z1 > Z) Z = z1; + if (r1 < r) { + r = r1; + } + if (R1 > R) { + R = R1; + } + if (z1 < z) { + z = z1; + } + if (z1 > Z) { + Z = z1; + } } tube.r = r; @@ -436,10 +461,18 @@ Int_t CbmKF::ReadTube(CbmKFTube& tube, FairGeoNode* node) double z1 = P->At(4 + i * 3 + 0); double r1 = P->At(4 + i * 3 + 1); double R1 = P->At(4 + i * 3 + 2); - if (r1 < r) r = r1; - if (R1 > R) R = R1; - if (z1 < z) z = z1; - if (z1 > Z) Z = z1; + if (r1 < r) { + r = r1; + } + if (R1 > R) { + R = R1; + } + if (z1 < z) { + z = z1; + } + if (z1 > Z) { + Z = z1; + } } tube.r = r; tube.R = R; @@ -471,7 +504,9 @@ Int_t CbmKF::ReadTube(CbmKFTube& tube, FairGeoNode* node) CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) { - if (!node) return nullptr; + if (!node) { + return nullptr; + } TString name = node->getName(); TString Sname = node->getShapePointer()->GetName(); @@ -480,7 +515,9 @@ CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) FairGeoNode* nxt = node; while ((nxt = nxt->getMotherNode())) { FairGeoTransform* tm = nxt->getLabTransform(); - if (!tm) break; + if (!tm) { + break; + } trans.transFrom(*tm); } @@ -521,10 +558,18 @@ CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) double z1 = P->At(3 + i * 3 + 0); double r1 = P->At(3 + i * 3 + 1); double R1 = P->At(3 + i * 3 + 2); - if (r1 < r) r = r1; - if (R1 > R) R = R1; - if (z1 < z) z = z1; - if (z1 > Z) Z = z1; + if (r1 < r) { + r = r1; + } + if (R1 > R) { + R = R1; + } + if (z1 < z) { + z = z1; + } + if (z1 > Z) { + Z = z1; + } } CbmKFTube tube(ID, x0, y0, z0 + 0.5 * (z + Z), (Z - z), r, R, RadLength); vPassiveTube.push_back(tube); @@ -537,10 +582,18 @@ CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) double z1 = P->At(4 + i * 3 + 0); double r1 = P->At(4 + i * 3 + 1); double R1 = P->At(4 + i * 3 + 2); - if (r1 < r) r = r1; - if (R1 > R) R = R1; - if (z1 < z) z = z1; - if (z1 > Z) Z = z1; + if (r1 < r) { + r = r1; + } + if (R1 > R) { + R = R1; + } + if (z1 < z) { + z = z1; + } + if (z1 > Z) { + Z = z1; + } } CbmKFTube tube(ID, x0, y0, z0 + 0.5 * (z + Z), (Z - z), r, R, RadLength); vPassiveTube.push_back(tube); @@ -557,8 +610,12 @@ CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) for (int i = 0; i < np; i++) { FairGeoVector* v = node->getPoint(i); for (int j = 0; j < 3; j++) { - if (vMin(j) > (*v)(j)) (&vMin.X())[j] = (*v)(j); - if (vMax(j) < (*v)(j)) (&vMax.X())[j] = (*v)(j); + if (vMin(j) > (*v)(j)) { + (&vMin.X())[j] = (*v)(j); + } + if (vMax(j) < (*v)(j)) { + (&vMax.X())[j] = (*v)(j); + } } } FairGeoVector v0 = (vMin + vMax); @@ -582,14 +639,18 @@ CbmKFMaterial* CbmKF::ReadPassive(FairGeoNode* node) Int_t CbmKF::Propagate(Double_t* T, Double_t* C, Double_t z_out, Double_t QP0) { Bool_t err = 0; - if (fabs(T[5] - z_out) < 1.e-5) return 0; + if (fabs(T[5] - z_out) < 1.e-5) { + return 0; + } if (!fMagneticField || (300 <= z_out && 300 <= T[5])) { CbmKFFieldMath::ExtrapolateLine(T, C, T, C, z_out); return 0; } Double_t zz = z_out; - if (z_out < 300. && 300 <= T[5]) CbmKFFieldMath::ExtrapolateLine(T, C, T, C, 300.); + if (z_out < 300. && 300 <= T[5]) { + CbmKFFieldMath::ExtrapolateLine(T, C, T, C, 300.); + } if (T[5] < 300. && 300. < z_out) { zz = 300.; @@ -598,8 +659,9 @@ Int_t CbmKF::Propagate(Double_t* T, Double_t* C, Double_t z_out, Double_t QP0) while (!err && repeat) { const Double_t max_step = 5.; Double_t zzz; - if (fabs(T[5] - zz) > max_step) + if (fabs(T[5] - zz) > max_step) { zzz = T[5] + ((zz > T[5]) ? max_step : -max_step); + } else { zzz = zz; repeat = 0; @@ -641,7 +703,9 @@ Int_t CbmKF::Propagate(Double_t* T, Double_t* C, Double_t z_out, Double_t QP0) */ } } - if (T[5] != z_out) CbmKFFieldMath::ExtrapolateLine(T, C, T, C, z_out); + if (T[5] != z_out) { + CbmKFFieldMath::ExtrapolateLine(T, C, T, C, z_out); + } return err; } diff --git a/reco/KF/CbmKFFieldMath.cxx b/reco/KF/CbmKFFieldMath.cxx index 81cf94c736025e0db3ade50b01447b20296b3212..64a240e98eb0b503a7db7fac6cb8b04dc998ff08 100644 --- a/reco/KF/CbmKFFieldMath.cxx +++ b/reco/KF/CbmKFFieldMath.cxx @@ -17,7 +17,9 @@ ClassImp(CbmKFFieldMath) void CbmKFFieldMath::ExtrapolateLine(const Double_t T_in[], const Double_t C_in[], Double_t T_out[], Double_t C_out[], Double_t z_out) { - if (!T_in) return; + if (!T_in) { + return; + } Double_t dz = z_out - T_in[5]; @@ -131,8 +133,9 @@ Int_t CbmKFFieldMath::ExtrapolateRK4(const Double_t T_in[], // input track para Double_t point[3] = {x[0], x[1], z_in + a[step] * h}; Double_t B[3]; - if (MF) + if (MF) { MF->GetFieldValue(point, B); + } else { B[0] = B[1] = B[2] = 0.; } @@ -143,7 +146,9 @@ Int_t CbmKFFieldMath::ExtrapolateRK4(const Double_t T_in[], // input track para Double_t ty2 = ty * ty; Double_t txty = tx * ty; Double_t tx2ty21 = 1.0 + tx2 + ty2; - if (tx2ty21 > 1.e4) return 1; + if (tx2ty21 > 1.e4) { + return 1; + } Double_t I_tx2ty21 = 1.0 / tx2ty21 * qp0; Double_t tx2ty2 = sqrt(tx2ty21); // Double_t I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ??? @@ -298,15 +303,16 @@ Int_t CbmKFFieldMath::ExtrapolateRK4(const Double_t T_in[], // input track para Double_t dqp = qp_in - qp0; - { - for (Int_t ip = 0; ip < 4; ip++) { - T_out[ip] += J[5 * 4 + ip] * dqp; - } + for (Int_t ip = 0; ip < 4; ip++) { + T_out[ip] += J[5 * 4 + ip] * dqp; } // covariance matrix transport - if (C_in && C_out) CbmKFMath::multQtSQ(5, J, C_in, C_out); + if (C_in && C_out) { + CbmKFMath::multQtSQ(5, J, C_in, C_out); + } + return 0; } // end of RK4 @@ -356,51 +362,63 @@ void CbmKFFieldMath::IntegrateField(CbmMagField* MF, Double_t r0[], Double_t r1[ for (Int_t x = 0; x < 3; x++) { if (si) { si[x] = 0; - for (Int_t n = 0; n < 3; n++) + for (Int_t n = 0; n < 3; n++) { si[x] += c1[n] * B[n][x]; + } si[x] *= dz / 6.; } if (Si) { Si[x] = 0; - for (Int_t n = 0; n < 3; n++) + for (Int_t n = 0; n < 3; n++) { Si[x] += C1[n] * B[n][x]; + } Si[x] *= dz * dz / 6.; } for (Int_t y = 0; y < 3; y++) { if (sii) { sii[x][y] = 0; - for (Int_t n = 0; n < 3; n++) - for (Int_t m = 0; m < 3; m++) + for (Int_t n = 0; n < 3; n++) { + for (Int_t m = 0; m < 3; m++) { sii[x][y] += c2[n][m] * B[n][x] * B[m][y]; + } + } sii[x][y] *= dz * dz / 360.; } if (Sii) { Sii[x][y] = 0; - for (Int_t n = 0; n < 3; n++) - for (Int_t m = 0; m < 3; m++) + for (Int_t n = 0; n < 3; n++) { + for (Int_t m = 0; m < 3; m++) { Sii[x][y] += C2[n][m] * B[n][x] * B[m][y]; + } + } Sii[x][y] *= dz * dz * dz / 2520.; } for (Int_t z = 0; z < 3; z++) { if (siii) { siii[x][y][z] = 0; - for (Int_t n = 0; n < 3; n++) - for (Int_t m = 0; m < 3; m++) - for (Int_t k = 0; k < 3; k++) + for (Int_t n = 0; n < 3; n++) { + for (Int_t m = 0; m < 3; m++) { + for (Int_t k = 0; k < 3; k++) { siii[x][y][z] += c3[n][m][k] * B[n][x] * B[m][y] * B[k][z]; + } + } + } siii[x][y][z] *= dz * dz * dz / 45360.; } if (Siii) { Siii[x][y][z] = 0; - for (Int_t n = 0; n < 3; n++) - for (Int_t m = 0; m < 3; m++) - for (Int_t k = 0; k < 3; k++) + for (Int_t n = 0; n < 3; n++) { + for (Int_t m = 0; m < 3; m++) { + for (Int_t k = 0; k < 3; k++) { Siii[x][y][z] += C3[n][m][k] * B[n][x] * B[m][y] * B[k][z]; + } + } + } Siii[x][y][z] *= dz * dz * dz * dz / 90720.; } } @@ -525,38 +543,44 @@ void CbmKFFieldMath::GetCoefficients(const Double_t x, const Double_t y, Double_ } } if (Xii) { - for (Int_t i = 0; i < 3; i++) + for (Int_t i = 0; i < 3; i++) { for (Int_t j = 0; j < 3; j++) { Xii[0][i][j] = X2[i][j]; Xii[1][i][j] = X2x[i][j]; Xii[2][i][j] = X2y[i][j]; } + } } if (Yii) { - for (Int_t i = 0; i < 3; i++) + for (Int_t i = 0; i < 3; i++) { for (Int_t j = 0; j < 3; j++) { Yii[0][i][j] = Y2[i][j]; Yii[1][i][j] = Y2x[i][j]; Yii[2][i][j] = Y2y[i][j]; } + } } if (Xiii) { - for (Int_t i = 0; i < 3; i++) - for (Int_t j = 0; j < 3; j++) + for (Int_t i = 0; i < 3; i++) { + for (Int_t j = 0; j < 3; j++) { for (Int_t k = 0; k < 3; k++) { Xiii[0][i][j][k] = X3[i][j][k]; Xiii[1][i][j][k] = X3x[i][j][k]; Xiii[2][i][j][k] = X3y[i][j][k]; } + } + } } if (Yiii) { - for (Int_t i = 0; i < 3; i++) - for (Int_t j = 0; j < 3; j++) + for (Int_t i = 0; i < 3; i++) { + for (Int_t j = 0; j < 3; j++) { for (Int_t k = 0; k < 3; k++) { Yiii[0][i][j][k] = Y3[i][j][k]; Yiii[1][i][j][k] = Y3x[i][j][k]; Yiii[2][i][j][k] = Y3y[i][j][k]; } + } + } } } @@ -661,7 +685,7 @@ void CbmKFFieldMath::ExtrapolateAnalytic(const Double_t T_in[], // input track SB1x += S1[n] * B1[1][n]; SB1y += S1[n] * B1[2][n]; - if (order > 1) + if (order > 1) { for (Int_t m = 0; m < 3; m++) { sA2 += s2[n][m] * A2[0][n][m]; sA2x += s2[n][m] * A2[1][n][m]; @@ -677,7 +701,7 @@ void CbmKFFieldMath::ExtrapolateAnalytic(const Double_t T_in[], // input track SB2x += S2[n][m] * B2[1][n][m]; SB2y += S2[n][m] * B2[2][n][m]; - if (order > 2) + if (order > 2) { for (Int_t k = 0; k < 3; k++) { sA3 += s3[n][m][k] * A3[0][n][m][k]; sA3x += s3[n][m][k] * A3[1][n][m][k]; @@ -693,7 +717,9 @@ void CbmKFFieldMath::ExtrapolateAnalytic(const Double_t T_in[], // input track SB3x += S3[n][m][k] * B3[1][n][m][k]; SB3y += S3[n][m][k] * B3[2][n][m][k]; } + } } + } } } @@ -756,15 +782,15 @@ void CbmKFFieldMath::ExtrapolateAnalytic(const Double_t T_in[], // input track Double_t dqp = qp_in - qp0; - { - for (Int_t i = 0; i < 4; i++) { - T_out[i] += J[5 * 4 + i] * dqp; - } + for (Int_t i = 0; i < 4; i++) { + T_out[i] += J[5 * 4 + i] * dqp; } // covariance matrix transport - if (C_in && C_out) CbmKFMath::multQtSQ(5, J, C_in, C_out); + if (C_in && C_out) { + CbmKFMath::multQtSQ(5, J, C_in, C_out); + } } @@ -846,23 +872,23 @@ void CbmKFFieldMath::ExtrapolateACentral(const Double_t T_in[], // input track Double_t IX1x = 0, IY1x = 0; Double_t IX1y = 0, IY1y = 0; - { - for (Int_t n = 0; n < 3; n++) { - if (n != 1) continue; - iX1 += i1[n] * X1[n]; - iX1x += i1[n] * X1x[n]; - iX1y += i1[n] * X1y[n]; - iY1 += i1[n] * Y1[n]; - iY1x += i1[n] * Y1x[n]; - iY1y += i1[n] * Y1y[n]; - - IX1 += I1[n] * X1[n]; - IX1x += I1[n] * X1x[n]; - IX1y += I1[n] * X1y[n]; - IY1 += I1[n] * Y1[n]; - IY1x += I1[n] * Y1x[n]; - IY1y += I1[n] * Y1y[n]; + for (Int_t n = 0; n < 3; n++) { + if (n != 1) { + continue; } + iX1 += i1[n] * X1[n]; + iX1x += i1[n] * X1x[n]; + iX1y += i1[n] * X1y[n]; + iY1 += i1[n] * Y1[n]; + iY1x += i1[n] * Y1x[n]; + iY1y += i1[n] * Y1y[n]; + + IX1 += I1[n] * X1[n]; + IX1x += I1[n] * X1x[n]; + IX1y += I1[n] * X1y[n]; + IY1 += I1[n] * Y1[n]; + IY1x += I1[n] * Y1x[n]; + IY1y += I1[n] * Y1y[n]; } Double_t t2 = 1. + xx + yy, t = sqrt(t2), h = qp0 * c_light, ht = h * t; @@ -925,15 +951,15 @@ void CbmKFFieldMath::ExtrapolateACentral(const Double_t T_in[], // input track Double_t dqp = qp_in - qp0; - { - for (Int_t i = 0; i < 4; i++) { - T_out[i] += J[5 * 4 + i] * dqp; - } + for (Int_t i = 0; i < 4; i++) { + T_out[i] += J[5 * 4 + i] * dqp; } // covariance matrix transport - if (C_in && C_out) CbmKFMath::multQtSQ(5, J, C_in, C_out); + if (C_in && C_out) { + CbmKFMath::multQtSQ(5, J, C_in, C_out); + } } #endif @@ -958,17 +984,22 @@ Int_t CbmKFFieldMath::ExtrapolateALight(const Double_t T_in[], // input track p // { bool ok = 1; - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { ok = ok && std::isfinite(T_in[i]) && (T_in[i] < 1.e5); - if (C_in) - for (int i = 0; i < 15; i++) + } + if (C_in) { + for (int i = 0; i < 15; i++) { ok = ok && std::isfinite(C_in[i]); + } + } if (!ok) { - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { T_out[i] = 0; + } if (C_out) { - for (int i = 0; i < 15; i++) + for (int i = 0; i < 15; i++) { C_out[i] = 0; + } C_out[0] = C_out[2] = C_out[5] = C_out[9] = C_out[14] = 100.; } return 1; @@ -1004,7 +1035,9 @@ Int_t CbmKFFieldMath::ExtrapolateALight(const Double_t T_in[], // input track p Double_t t2 = 1. + xx + yy; - if (t2 > 1.e4) return 1; + if (t2 > 1.e4) { + return 1; + } Double_t t = sqrt(t2), h = qp0 * c_light, ht = h * t; Double_t sx = 0, sy = 0, sz = 0, syy = 0, syz = 0, syyy = 0, Sx = 0, Sy = 0, Sz = 0, Syy = 0, Syz = 0, Syyy = 0; @@ -1058,11 +1091,12 @@ Int_t CbmKFFieldMath::ExtrapolateALight(const Double_t T_in[], // input track p Double_t c2[3][3] = {{5, -4, -1}, {44, 80, -4}, {11, 44, 5}}; // /=360. Double_t C2[3][3] = {{38, 8, -4}, {148, 208, -20}, {3, 36, 3}}; // /=2520. - for (Int_t n = 0; n < 3; n++) + for (Int_t n = 0; n < 3; n++) { for (Int_t m = 0; m < 3; m++) { syz += c2[n][m] * B[n][1] * B[m][2]; Syz += C2[n][m] * B[n][1] * B[m][2]; } + } syz *= dz * dz / 360.; Syz *= dz * dz * dz / 2520.; @@ -1170,14 +1204,15 @@ Int_t CbmKFFieldMath::ExtrapolateALight(const Double_t T_in[], // input track p Double_t dqp = qp_in - qp0; - { - for (Int_t i = 0; i < 4; i++) { - T_out[i] += J[5 * 4 + i] * dqp; - } + for (Int_t i = 0; i < 4; i++) { + T_out[i] += J[5 * 4 + i] * dqp; } // covariance matrix transport - if (C_in && C_out) CbmKFMath::multQtSQ(5, J, C_in, C_out); + if (C_in && C_out) { + CbmKFMath::multQtSQ(5, J, C_in, C_out); + } + return 0; } diff --git a/reco/KF/CbmKFMaterial.cxx b/reco/KF/CbmKFMaterial.cxx index 2ac3f91968ed400ffac90229e6d59e70631b183e..560d3b49f5e8064fd99f8d2afa87499597db28bd 100644 --- a/reco/KF/CbmKFMaterial.cxx +++ b/reco/KF/CbmKFMaterial.cxx @@ -29,13 +29,19 @@ Int_t CbmKFMaterial::Pass(Double_t ZCross, Double_t ZThick, CbmKFTrackInterface& Double_t* T = track.GetTrack(); Double_t* C = track.GetCovMatrix(); err = err || CbmKF::Instance()->Propagate(T, C, ZCross, QP0); - if (err) return err; - if (IsOutside(T[0], T[1])) return 0; + if (err) { + return err; + } + if (IsOutside(T[0], T[1])) { + return 0; + } Double_t Q5, Q8, Q9, Ecor; err = err || CbmKFMath::GetNoise(ZThick / RadLength, F, Fe, T[2], T[3], QP0, track.GetMass(), track.IsElectron(), downstream, &Q5, &Q8, &Q9, &Ecor); - if (err) return err; + if (err) { + return err; + } C[5] += Q5; C[8] += Q8; C[9] += Q9; @@ -121,19 +127,30 @@ Int_t CbmKFCone::Pass(Double_t ZCross, Double_t /*ZThick*/, CbmKFTrackInterface& Double_t zthick = ZThickness, zcross = ZReference; Double_t T_tmp[6]; - for (Int_t j = 0; j < 6; j++) + for (Int_t j = 0; j < 6; j++) { T_tmp[j] = T[j]; + } err = err || KF->Propagate(T_tmp, nullptr, ZReference, QP0); - if (err) return err; + if (err) { + return err; + } Double_t cz1, ct1, cz2, ct2; { Double_t iz1, iz2, iZ1, iZ2; Bool_t err1 = CbmKFMath::intersectCone(z1, z2, r1, r2, T_tmp, &iz1, &iz2); Bool_t err2 = CbmKFMath::intersectCone(z1, z2, R1, R2, T_tmp, &iZ1, &iZ2); - if (err1 || iz1 < z1 || iz1 > z2) iz1 = -200; - if (err1 || iz2 < z1 || iz2 > z2) iz2 = -200; - if (err2 || iZ1 < z1 || iZ1 > z2) iZ1 = -200; - if (err2 || iZ2 < z1 || iZ2 > z2) iZ2 = -200; + if (err1 || iz1 < z1 || iz1 > z2) { + iz1 = -200; + } + if (err1 || iz2 < z1 || iz2 > z2) { + iz2 = -200; + } + if (err2 || iZ1 < z1 || iZ1 > z2) { + iZ1 = -200; + } + if (err2 || iZ2 < z1 || iZ2 > z2) { + iZ2 = -200; + } if (iz1 > -200 && iZ1 > -200) { cz1 = (iz1 + iZ1) / 2; @@ -177,7 +194,9 @@ Int_t CbmKFCone::Pass(Double_t ZCross, Double_t /*ZThick*/, CbmKFTrackInterface& err = err || CbmKFMath::GetNoise(zthick / RadLength, F, Fe, T[2], T[3], QP0, track.GetMass(), track.IsElectron(), downstream, &Q5, &Q8, &Q9, &Ecor); - if (err) return err; + if (err) { + return err; + } C[5] += Q5; C[8] += Q8; C[9] += Q9; diff --git a/reco/KF/CbmKFMath.cxx b/reco/KF/CbmKFMath.cxx index a331d174c67e2ae35347aaea0161bb2a3ddb6bfa..cd8d24ae678f540a67778926c24e726e79f8f07d 100644 --- a/reco/KF/CbmKFMath.cxx +++ b/reco/KF/CbmKFMath.cxx @@ -29,9 +29,13 @@ ClassImp(CbmKFMath) Double_t b = 2 * (tx * x0 + ty * y0 + t * A); Double_t c = x0 * x0 + y0 * y0 - A * A; - if (fabs(a) < 1.E-4) return 1; + if (fabs(a) < 1.E-4) { + return 1; + } Double_t D = b * b - 4 * a * c; - if (D < 0.) return 1; + if (D < 0.) { + return 1; + } D = sqrt(D); *z1 = (-b + D) / (2 * a); *z2 = (-b - D) / (2 * a); @@ -43,23 +47,21 @@ void CbmKFMath::multQSQt(Int_t N, const Double_t Q[], const Double_t S[], Double { Double_t A[N * N]; - { - for (Int_t i = 0, n = 0; i < N; i++) { - for (Int_t j = 0; j < N; j++, ++n) { - A[n] = 0; - for (Int_t k = 0; k < N; ++k) - A[n] += S[indexS(i, k)] * Q[j * N + k]; + for (Int_t i = 0, n = 0; i < N; i++) { + for (Int_t j = 0; j < N; j++, ++n) { + A[n] = 0; + for (Int_t k = 0; k < N; ++k) { + A[n] += S[indexS(i, k)] * Q[j * N + k]; } } } - { - for (Int_t i = 0; i < N; i++) { - for (Int_t j = 0; j <= i; j++) { - Int_t n = indexS(i, j); - S_out[n] = 0; - for (Int_t k = 0; k < N; k++) - S_out[n] += Q[i * N + k] * A[k * N + j]; + for (Int_t i = 0; i < N; i++) { + for (Int_t j = 0; j <= i; j++) { + Int_t n = indexS(i, j); + S_out[n] = 0; + for (Int_t k = 0; k < N; k++) { + S_out[n] += Q[i * N + k] * A[k * N + j]; } } } @@ -72,8 +74,9 @@ void CbmKFMath::multQtSQ(Int_t N, const Double_t Q[], const Double_t S[], Double for (Int_t i = 0, n = 0; i < N; i++) { for (Int_t j = 0; j < N; j++, ++n) { A[n] = 0; - for (Int_t k = 0; k < N; ++k) + for (Int_t k = 0; k < N; ++k) { A[n] += S[indexS(i, k)] * Q[k * N + j]; + } } } @@ -81,8 +84,9 @@ void CbmKFMath::multQtSQ(Int_t N, const Double_t Q[], const Double_t S[], Double for (Int_t j = 0; j <= i; j++) { Int_t n = indexS(i, j); S_out[n] = 0; - for (Int_t k = 0; k < N; k++) + for (Int_t k = 0; k < N; k++) { S_out[n] += Q[k * N + i] * A[k * N + j]; + } } } } @@ -93,8 +97,9 @@ void CbmKFMath::multSSQ(const Double_t* A, const Double_t* B, Double_t* C, Int_t for (Int_t j = 0; j < n; ++j) { Int_t ind = i * n + j; C[ind] = 0; - for (Int_t k = 0; k < n; ++k) + for (Int_t k = 0; k < n; ++k) { C[ind] += A[indexS(i, k)] * B[indexS(k, j)]; + } } } } @@ -108,29 +113,35 @@ void CbmKFMath::four_dim_inv(Double_t a[4][4]) Double_t b[4][4]; // Set b to I - for (i = 0; i < 4; i++) - for (j = 0; j < 4; j++) - if (i == j) - b[i][j] = 1.0; - else - b[i][j] = 0.0; + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + b[i][j] = (i == j) ? 1.0 : 0.0; + } + } for (i = 0; i < 4; i++) { - for (j = i + 1; j < 4; j++) + for (j = i + 1; j < 4; j++) { if (fabs(a[i][i]) < fabs(a[j][i])) { - for (l = 0; l < 4; l++) + for (l = 0; l < 4; l++) { temp[l] = a[i][l]; - for (l = 0; l < 4; l++) + } + for (l = 0; l < 4; l++) { a[i][l] = a[j][l]; - for (l = 0; l < 4; l++) + } + for (l = 0; l < 4; l++) { a[j][l] = temp[l]; - for (l = 0; l < 4; l++) + } + for (l = 0; l < 4; l++) { temp[l] = b[i][l]; - for (l = 0; l < 4; l++) + } + for (l = 0; l < 4; l++) { b[i][l] = b[j][l]; - for (l = 0; l < 4; l++) + } + for (l = 0; l < 4; l++) { b[j][l] = temp[l]; + } } + } factor = a[i][i]; for (j = 4 - 1; j > -1; j--) { b[i][j] /= factor; @@ -155,9 +166,11 @@ void CbmKFMath::four_dim_inv(Double_t a[4][4]) } // copy b to a and return - for (i = 0; i < 4; i++) - for (j = 0; j < 4; j++) + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { a[i][j] = b[i][j]; + } + } } void CbmKFMath::five_dim_inv(Double_t a[5][5]) @@ -171,30 +184,34 @@ void CbmKFMath::five_dim_inv(Double_t a[5][5]) // Set b to I for (i = 0; i < 5; i++) { for (j = 0; j < 5; j++) { - if (i == j) - b[i][j] = 1.0; - else - b[i][j] = 0.0; + b[i][j] = (i == j) ? 1. : 0.; } } for (i = 0; i < 5; i++) { - for (j = i + 1; j < 5; j++) + for (j = i + 1; j < 5; j++) { if (fabs(a[i][i]) < fabs(a[j][i])) { - for (l = 0; l < 5; l++) + for (l = 0; l < 5; l++) { temp[l] = a[i][l]; - for (l = 0; l < 5; l++) + } + for (l = 0; l < 5; l++) { a[i][l] = a[j][l]; - for (l = 0; l < 5; l++) + } + for (l = 0; l < 5; l++) { a[j][l] = temp[l]; - for (l = 0; l < 5; l++) + } + for (l = 0; l < 5; l++) { temp[l] = b[i][l]; - for (l = 0; l < 5; l++) + } + for (l = 0; l < 5; l++) { b[i][l] = b[j][l]; - for (l = 0; l < 5; l++) + } + for (l = 0; l < 5; l++) { b[j][l] = temp[l]; + } } + } factor = a[i][i]; // cout<<"Highest element "<<a[i][i]<<endl; for (j = 5 - 1; j > -1; j--) { @@ -220,9 +237,11 @@ void CbmKFMath::five_dim_inv(Double_t a[5][5]) } // copy b to a and return - for (i = 0; i < 5; i++) - for (j = 0; j < 5; j++) + for (i = 0; i < 5; i++) { + for (j = 0; j < 5; j++) { a[i][j] = b[i][j]; + } + } } Bool_t CbmKFMath::invS(Double_t A[], Int_t N) @@ -256,15 +275,17 @@ Bool_t CbmKFMath::invS(Double_t A[], Int_t N) x = 1 / x; for (Int_t step = 1; step <= N - j; ik += ++step) { // ik==Ai1 Double_t sum = 0; - for (Double_t* jk = j1; jk != jj; sum += *(jk++) * *(ik++)) + for (Double_t* jk = j1; jk != jj; sum += *(jk++) * *(ik++)) { ; + } *ik = (*ik - sum) * x; // ik == Aij } } else { Double_t* ji = jj; - for (Int_t i = j; i < N; i++) + for (Int_t i = j; i < N; i++) { *(ji += i) = 0.; + } ret = 1; } } @@ -281,16 +302,14 @@ Bool_t CbmKFMath::invS(Double_t A[], Int_t N) Double_t *ii = A, *ij = A; for (Int_t i = 1; i <= N; ij = ii + 1, ii += ++i) { if (*ii > ZERO) { - Double_t x = -(*ii = 1. / *ii); - { - Double_t* jj = A; - for (Int_t j = 1; j < i; jj += ++j, ij++) { - Double_t *ik = ij, *kj = jj, sum = 0.; - for (Int_t k = j; ik != ii; kj += k++, ik++) { - sum += *ik * *kj; - } - *kj = sum * x; + Double_t x = -(*ii = 1. / *ii); + Double_t* jj = A; + for (Int_t j = 1; j < i; jj += ++j, ij++) { + Double_t *ik = ij, *kj = jj, sum = 0.; + for (Int_t k = j; ik != ii; kj += k++, ik++) { + sum += *ik * *kj; } + *kj = sum * x; } } else { @@ -311,8 +330,9 @@ Bool_t CbmKFMath::invS(Double_t A[], Int_t N) for (Int_t i = 1; i <= N; ii += ++i) { do { Double_t *ki = ii, *kj = ij, sum = 0.; - for (Int_t k = i; k <= N; ki += k, kj += k++) + for (Int_t k = i; k <= N; ki += k, kj += k++) { sum += (*ki) * (*kj); + } *ij = sum; } while ((ij++) != ii); } @@ -336,7 +356,9 @@ Double_t CbmKFMath::getDeviation(Double_t x, Double_t y, Double_t C[], Double_t c[2] += Cv[2]; } Double_t d = c[0] * c[2] - c[1] * c[1]; - if (fabs(d) < 1.e-20) return 0; + if (fabs(d) < 1.e-20) { + return 0; + } return sqrt(fabs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d)); } @@ -375,7 +397,9 @@ Double_t CbmKFMath::AnalyticQP(const Double_t T[], // track parameters (x, Double_t B[3][3]; Int_t n = int(fabs(vz - z) / 5.); - if (n < 2) n = 2; + if (n < 2) { + n = 2; + } Double_t dz = (vz - z) / n; { @@ -399,13 +423,14 @@ Double_t CbmKFMath::AnalyticQP(const Double_t T[], // track parameters (x, Double_t c2[3][3] = {{5, -4, -1}, {44, 80, -4}, {11, 44, 5}}; // /=360. Double_t C2[3][3] = {{38, 8, -4}, {148, 208, -20}, {3, 36, 3}}; // /=2520. - for (Int_t k = 0; k < 3; k++) + for (Int_t k = 0; k < 3; k++) { for (Int_t m = 0; m < 3; m++) { syz += c2[k][m] * B[k][1] * B[m][2]; Syz += C2[k][m] * B[k][1] * B[m][2]; szy += c2[k][m] * B[k][2] * B[m][1]; Szy += C2[k][m] * B[k][2] * B[m][1]; } + } syz *= dz * dz * 4 / 360.; Syz *= dz * dz * dz * 8 / 2520.; @@ -486,12 +511,12 @@ Double_t CbmKFMath::AnalyticQP(const Double_t T[], // track parameters (x, Double_t c2[3][3] = {{5, -52, -13}, {-28, 320, 68}, {-37, 332, 125}}; // /=5760 Double_t C2[3][3] = {{13, -152, -29}, {-82, 1088, 170}, {-57, 576, 153}}; // /=80640 - { - for (Int_t k = 0; k < 3; k++) - for (Int_t m = 0; m < 3; m++) { - Syz_ += C2[k][m] * B[k][1] * B[m][2]; - Szy_ += C2[k][m] * B[k][2] * B[m][1]; - } + + for (Int_t k = 0; k < 3; k++) { + for (Int_t m = 0; m < 3; m++) { + Syz_ += C2[k][m] * B[k][1] * B[m][2]; + Szy_ += C2[k][m] * B[k][2] * B[m][1]; + } } Syz_ *= dz * dz * dz * 8 / 80640.; @@ -547,15 +572,15 @@ Double_t CbmKFMath::AnalyticQP(const Double_t T[], // track parameters (x, syz_ = Syz_ = szy_ = Szy_ = 0; - { - for (Int_t k = 0; k < 3; k++) - for (Int_t m = 0; m < 3; m++) { - syz_ += c2[k][m] * B[k][1] * B[m][2]; - Syz_ += C2[k][m] * B[k][1] * B[m][2]; - szy_ += c2[k][m] * B[k][2] * B[m][1]; - Szy_ += C2[k][m] * B[k][2] * B[m][1]; - } + for (Int_t k = 0; k < 3; k++) { + for (Int_t m = 0; m < 3; m++) { + syz_ += c2[k][m] * B[k][1] * B[m][2]; + Syz_ += C2[k][m] * B[k][1] * B[m][2]; + szy_ += c2[k][m] * B[k][2] * B[m][1]; + Szy_ += C2[k][m] * B[k][2] * B[m][1]; + } } + syz_ *= dz * dz * 4 / 5760.; Syz_ *= dz * dz * dz * 8 / 80640.; @@ -667,7 +692,9 @@ Double_t CbmKFMath::AnalyticQP(const Double_t T[], // track parameters (x, A = 3 * d + a; Double_t D = B * B - 4 * A * C; - if (D < 0) D = 0; + if (D < 0.) { + D = 0.; + } D = sqrt(D); Double_t dqp1 = (-B + D) / 2 / A; @@ -693,7 +720,9 @@ Bool_t CbmKFMath::GetThickness(Double_t z1, Double_t z2, Double_t mz, Double_t m } Double_t tmin = mz - mthick / 2, tmax = mz + mthick / 2; - if (tmin >= Z || z >= tmax) return 1; + if (tmin >= Z || z >= tmax) { + return 1; + } if (z <= tmin && tmax <= Z) // case z |**| Z { *mz_out = mz; @@ -727,10 +756,14 @@ Int_t CbmKFMath::GetNoise(Double_t Lrl, Double_t F, Double_t Fe, Double_t tx, Do *Q9 = 0; *Ecor = 1.; Double_t t = sqrt(1. + tx * tx + ty * ty); - if (!isfinite(t) || t > 1.e4) return 1; + if (!isfinite(t) || t > 1.e4) { + return 1; + } t = sqrt(t); Double_t l = t * Lrl; - if (l > 1.) l = 1.; // protect against l too big + if (l > 1.) { + l = 1.; // protect against l too big + } Double_t s0 = (l > exp(-1. / 0.038)) ? F * .0136 * (1 + 0.038 * log(l)) * qp : 0.; Double_t a = (1. + mass * mass * qp * qp) * s0 * s0 * t * t * l; @@ -752,10 +785,12 @@ Int_t CbmKFMath::GetNoise(Double_t Lrl, Double_t F, Double_t Fe, Double_t tx, Do Double_t m_energyLoss = Fe; // Double_t m_energyLoss = 0.02145;//0.060; Double_t corr = (1. - fabs(qp) * L * m_energyLoss); - if (corr > 0.001 * fabs(qp)) + if (corr > 0.001 * fabs(qp)) { *Ecor *= 1. / corr; - else + } + else { *Ecor = 1000. / fabs(qp); + } } return 0; } @@ -772,9 +807,11 @@ void CbmKFMath::CopyTC2TrackParam(FairTrackParam* par, Double_t T[], Double_t C[ par->SetQp(T[4]); } if (C) { - for (Int_t i = 0, iCov = 0; i < 5; i++) - for (Int_t j = 0; j <= i; j++, iCov++) + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { par->SetCovariance(i, j, C[iCov]); + } + } } } @@ -789,8 +826,10 @@ void CbmKFMath::CopyTrackParam2TC(const FairTrackParam* par, Double_t T[], Doubl T[5] = par->GetZ(); } if (C) { - for (Int_t i = 0, iCov = 0; i < 5; i++) - for (Int_t j = 0; j <= i; j++, iCov++) + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { C[iCov] = par->GetCovariance(i, j); + } + } } } diff --git a/reco/KF/CbmKFParticleFinder.cxx b/reco/KF/CbmKFParticleFinder.cxx index 64f847de8345ce14e0e289a8beaa0816f9971f19..defed1dfdd549d3390e1ae8e912594923bb5228a 100644 --- a/reco/KF/CbmKFParticleFinder.cxx +++ b/reco/KF/CbmKFParticleFinder.cxx @@ -60,7 +60,9 @@ CbmKFParticleFinder::CbmKFParticleFinder(const char* name, Int_t iVerbose) CbmKFParticleFinder::~CbmKFParticleFinder() { - if (fTopoReconstructor) delete fTopoReconstructor; + if (fTopoReconstructor) { + delete fTopoReconstructor; + } } InitStatus CbmKFParticleFinder::Init() @@ -79,14 +81,16 @@ InitStatus CbmKFParticleFinder::Init() fTimeSliceMode = 1; LOG(info) << GetName() << ": Running in the timeslice mode."; } - else + else { LOG(info) << GetName() << ": Running in the event by event mode."; + } // Get reconstructed events if (fTimeSliceMode) { fEvents = (TClonesArray*) ioman->GetObject("CbmEvent"); - if (fEvents == nullptr) + if (fEvents == nullptr) { Fatal("CbmKFParticleFinder::Init", "No events available. Running in the event-by-event mode."); + } } // Get input collection @@ -100,7 +104,9 @@ InitStatus CbmKFParticleFinder::Init() if (fPVFindMode == 0) { if (fTimeSliceMode) { CbmMCDataManager* mcManager = (CbmMCDataManager*) ioman->GetObject("MCDataManager"); - if (mcManager == nullptr) Error("CbmKFParticleFinder::Init", "MC Data Manager not found!"); + if (mcManager == nullptr) { + Error("CbmKFParticleFinder::Init", "MC Data Manager not found!"); + } fMCTrackArray = mcManager->InitBranch("MCTrack"); @@ -138,14 +144,18 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) fTopoReconstructor->Clear(); int nEvents = 1; - if (fTimeSliceMode) nEvents = fEvents->GetEntriesFast(); + if (fTimeSliceMode) { + nEvents = fEvents->GetEntriesFast(); + } vector<KFParticleTopoReconstructor> eventTopoReconstructor; eventTopoReconstructor.resize(nEvents); for (int iEvent = 0; iEvent < nEvents; iEvent++) { CbmEvent* event = nullptr; - if (fTimeSliceMode) event = (CbmEvent*) fEvents->At(iEvent); + if (fTimeSliceMode) { + event = (CbmEvent*) fEvents->At(iEvent); + } eventTopoReconstructor[iEvent].SetTarget(fTopoReconstructor->GetTargetPosition()); eventTopoReconstructor[iEvent].SetChi2PrimaryCut(InversedChi2Prob(0.0001, 2)); eventTopoReconstructor[iEvent].CopyCuts(fTopoReconstructor); @@ -153,10 +163,12 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) fTopoReconstructor->GetKFParticleFinder()->GetReconstructionList()); int nTracksEvent = 0; - if (fTimeSliceMode) + if (fTimeSliceMode) { nTracksEvent = event->GetNofStsTracks(); - else + } + else { nTracksEvent = fTrackArray->GetEntriesFast(); + } Int_t ntracks = 0; //fTrackArray->GetEntriesFast(); @@ -165,7 +177,9 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) if (fPID) { if ((int) fPID->GetPID().size() == nTracksEvent) { for (int iTr = 0; iTr < nTracksEvent; iTr++) { - if (TMath::Abs(fPID->GetPID()[iTr]) == 1000010029) nCandidatesDHe4++; + if (TMath::Abs(fPID->GetPID()[iTr]) == 1000010029) { + nCandidatesDHe4++; + } } } else { @@ -180,20 +194,26 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) for (int iTr = 0; iTr < nTracksEvent; iTr++) { int stsTrackIndex = 0; - if (fTimeSliceMode) + if (fTimeSliceMode) { stsTrackIndex = event->GetStsTrackIndex(iTr); - else + } + else { stsTrackIndex = iTr; + } CbmStsTrack* stsTrack = ((CbmStsTrack*) fTrackArray->At(stsTrackIndex)); const FairTrackParam* parameters = stsTrack->GetParamFirst(); Double_t V[15] = {0.f}; - for (Int_t i = 0, iCov = 0; i < 5; i++) - for (Int_t j = 0; j <= i; j++, iCov++) + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { V[iCov] = parameters->GetCovariance(i, j); + } + } - if (stsTrack->GetTotalNofHits() < 3) continue; + if (stsTrack->GetTotalNofHits() < 3) { + continue; + } bool ok = 1; ok = ok && std::isfinite(parameters->GetX()); @@ -203,15 +223,21 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) ok = ok && std::isfinite(parameters->GetTy()); ok = ok && std::isfinite(parameters->GetQp()); - for (unsigned short iC = 0; iC < 15; iC++) + for (unsigned short iC = 0; iC < 15; iC++) { ok = ok && std::isfinite(V[iC]); + } + ok = ok && (V[0] < 1. && V[0] > 0.) && (V[2] < 1. && V[2] > 0.) && (V[5] < 1. && V[5] > 0.) && (V[9] < 1. && V[9] > 0.) && (V[14] < 1. && V[14] > 0.); ok = ok && stsTrack->GetChiSq() < 10 * stsTrack->GetNDF(); - if (!ok) continue; + if (!ok) { + continue; + } if (fPID) { - if (fPID->GetPID()[stsTrackIndex] == -2) continue; + if (fPID->GetPID()[stsTrackIndex] == -2) { + continue; + } //not clear separation between d and He4 if (TMath::Abs(fPID->GetPID()[stsTrackIndex]) == 1000010029) { @@ -223,11 +249,13 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) pdg[ntracks] = sgn * 1000020040; } - else + else { pdg[ntracks] = fPID->GetPID()[stsTrackIndex]; + } } - else + else { pdg[ntracks] = -1; + } vRTracks[ntracks] = *stsTrack; trackId[ntracks] = stsTrackIndex; @@ -244,23 +272,29 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) CbmKFVertex kfVertex; if (fPVFindMode == 0) { int nMCEvents = 1; - if (fTimeSliceMode) nMCEvents = fEventList->GetNofEvents(); + if (fTimeSliceMode) { + nMCEvents = fEventList->GetNofEvents(); + } float mcpv[3] = {0.f, 0.f, 0.f}; bool isMCPVFound = false; for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) { int nMCTracks = 0; - if (fTimeSliceMode) + if (fTimeSliceMode) { nMCTracks = fMCTrackArray->Size(0, iMCEvent); - else + } + else { nMCTracks = fMCTrackArrayEvent->GetEntriesFast(); + } for (Int_t iMC = 0; iMC < nMCTracks; iMC++) { CbmMCTrack* cbmMCTrack; - if (fTimeSliceMode) + if (fTimeSliceMode) { cbmMCTrack = (CbmMCTrack*) fMCTrackArray->Get(0, iMCEvent, iMC); - else + } + else { cbmMCTrack = (CbmMCTrack*) fMCTrackArrayEvent->At(iMC); + } if (cbmMCTrack->GetMotherId() < 0) { mcpv[0] = cbmMCTrack->GetStartX(); @@ -270,7 +304,9 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) break; } } - if (isMCPVFound) break; + if (isMCPVFound) { + break; + } } kfVertex.GetRefX() = mcpv[0]; @@ -322,10 +358,12 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) KFVertex pv(primVtx_tmp); eventTopoReconstructor[iEvent].AddPV(pv, pvTrackIds); } - else if (fPVFindMode == 1) + else if (fPVFindMode == 1) { eventTopoReconstructor[iEvent].ReconstructPrimVertex(); - else if (fPVFindMode == 2) + } + else if (fPVFindMode == 2) { eventTopoReconstructor[iEvent].ReconstructPrimVertex(0); + } eventTopoReconstructor[iEvent].SortTracks(); eventTopoReconstructor[iEvent].ReconstructParticles(); @@ -338,8 +376,12 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) const FairTrackParam* parameters = vRTracks[iTr].GetParamFirst(); float a = parameters->GetTx(), b = parameters->GetTy(), qp = parameters->GetQp(); Int_t q = 0; - if (qp > 0.f) q = 1; - if (qp < 0.f) q = -1; + if (qp > 0.f) { + q = 1; + } + if (qp < 0.f) { + q = -1; + } float c2 = 1.f / (1.f + a * a + b * b); float pq = 1.f / qp * TMath::Abs(q); float p2 = pq * pq; @@ -351,12 +393,16 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) bool save = 0; if (vChiToPrimVtx[iTr] < 3) { - if ((fabs(pdg[iTr]) == 11 && pt > 0.2f) || (fabs(pdg[iTr]) == 13) || (fabs(pdg[iTr]) == 19)) save = 1; + if ((fabs(pdg[iTr]) == 11 && pt > 0.2f) || (fabs(pdg[iTr]) == 13) || (fabs(pdg[iTr]) == 19)) { + save = 1; + } } if (vChiToPrimVtx[iTr] > 3) { - if ((fabs(pdg[iTr]) == 211 || fabs(pdg[iTr]) == 321 || fabs(pdg[iTr]) == 2212 || pdg[iTr] == -1) && pt > 0.2f) + if ((fabs(pdg[iTr]) == 211 || fabs(pdg[iTr]) == 321 || fabs(pdg[iTr]) == 2212 || pdg[iTr] == -1) + && pt > 0.2f) { save = 1; + } } if (save) { @@ -379,8 +425,9 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) for (int iPV = 0; iPV < eventTR->NPrimaryVertices(); iPV++) { PVTrackIndexArray = eventTR->GetPVTrackIndexArray(iPV); - for (unsigned int iTr = 0; iTr < PVTrackIndexArray.size(); iTr++) + for (unsigned int iTr = 0; iTr < PVTrackIndexArray.size(); iTr++) { PVTrackIndexArray[iTr] = PVTrackIndexArray[iTr] + indexAdd; + } fTopoReconstructor->AddPV(eventTR->GetPrimKFVertex(iPV), PVTrackIndexArray); PVTrackIndexArray.clear(); @@ -393,8 +440,12 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/) int idP = particleEvent.Id() + indexAdd; int idDaughter = 0; for (int nD = 0; nD < particleEvent.NDaughters(); nD++) { - if (particleEvent.NDaughters() == 1) idDaughter = particleEvent.DaughterIds()[nD]; - if (particleEvent.NDaughters() > 1) idDaughter = particleEvent.DaughterIds()[nD] + indexAdd; + if (particleEvent.NDaughters() == 1) { + idDaughter = particleEvent.DaughterIds()[nD]; + } + if (particleEvent.NDaughters() > 1) { + idDaughter = particleEvent.DaughterIds()[nD] + indexAdd; + } particle.AddDaughterId(idDaughter); } particle.SetId(idP); @@ -445,19 +496,27 @@ void CbmKFParticleFinder::FillKFPTrackVector(KFPTrackVector* tracks, const vecto //fill vector with tracks for (Int_t iTr = 0; iTr < ntracks; iTr++) { const FairTrackParam* parameters; - if (atFirstPoint) + if (atFirstPoint) { parameters = vRTracks[iTr].GetParamFirst(); - else + } + else { parameters = vRTracks[iTr].GetParamLast(); + } double par[6] = {0.f}; double tx = parameters->GetTx(), ty = parameters->GetTy(), qp = parameters->GetQp(); Int_t q = 0; - if (qp > 0.f) q = 1; - if (qp < 0.f) q = -1; - if (TMath::Abs(pdg[iTr]) == 1000020030 || TMath::Abs(pdg[iTr]) == 1000020040) q *= 2; + if (qp > 0.f) { + q = 1; + } + if (qp < 0.f) { + q = -1; + } + if (TMath::Abs(pdg[iTr]) == 1000020030 || TMath::Abs(pdg[iTr]) == 1000020040) { + q *= 2; + } double c2 = 1.f / (1.f + tx * tx + ty * ty); @@ -492,42 +551,50 @@ void CbmKFParticleFinder::FillKFPTrackVector(KFPTrackVector* tracks, const vecto {0.f, 0.f, dpydtx, dpydty, dpydqp}, {0.f, 0.f, dpzdtx, dpzdty, dpzdqp}}; double VFT[5][6]; - for (int i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { for (int j = 0; j < 6; j++) { VFT[i][j] = 0; for (int k = 0; k < 5; k++) { VFT[i][j] += parameters->GetCovariance(i, k) * F[j][k]; } } + } double cov[21]; - for (int i = 0, l = 0; i < 6; i++) + for (int i = 0, l = 0; i < 6; i++) { for (int j = 0; j <= i; j++, l++) { cov[l] = 0; for (int k = 0; k < 5; k++) { cov[l] += F[i][k] * VFT[k][j]; } } + } - for (Int_t iP = 0; iP < 6; iP++) + for (Int_t iP = 0; iP < 6; iP++) { tracks->SetParameter(par[iP], iP, iTr); - for (Int_t iC = 0; iC < 21; iC++) + } + for (Int_t iC = 0; iC < 21; iC++) { tracks->SetCovariance(cov[iC], iC, iTr); - for (Int_t iF = 0; iF < 10; iF++) + } + for (Int_t iF = 0; iF < 10; iF++) { tracks->SetFieldCoefficient(vField[iTr].fField[iF], iF, iTr); + } tracks->SetId(trackId[iTr], iTr); tracks->SetPDG(pdg[iTr], iTr); tracks->SetQ(q, iTr); tracks->SetNPixelHits(vRTracks[iTr].GetNofMvdHits(), iTr); if (fPVFindMode == 0 || fPVFindMode == 3) { - if (vChiToPrimVtx[iTr] < 3) + if (vChiToPrimVtx[iTr] < 3) { tracks->SetPVIndex(0, iTr); - else + } + else { tracks->SetPVIndex(-1, iTr); + } } - else + else { tracks->SetPVIndex(-1, iTr); + } } } diff --git a/reco/KF/CbmKFParticleFinderPID.cxx b/reco/KF/CbmKFParticleFinderPID.cxx index 8376cbeb7ea72fd78ee1053174a394ba219dc3a3..178d024bc3c3ac6bbae99fecd3497a900866d1c5 100644 --- a/reco/KF/CbmKFParticleFinderPID.cxx +++ b/reco/KF/CbmKFParticleFinderPID.cxx @@ -33,16 +33,20 @@ #include <vector> using std::vector; +ClassImp(CbmKFParticleFinderPID); + double vecMedian(const vector<double>& vec) { double median = 0.; vector<double> vecCopy = vec; sort(vecCopy.begin(), vecCopy.end()); int size = vecCopy.size(); - if (size % 2 == 0) + if (size % 2 == 0) { median = (vecCopy[size / 2 - 1] + vecCopy[size / 2]) / 2; - else + } + else { median = vecCopy[size / 2]; + } return median; } @@ -109,9 +113,9 @@ InitStatus CbmKFParticleFinderPID::Init() fLegacyEventMode = 0; LOG(info) << GetName() << ": Running in the timeslice mode."; } - else + else { LOG(info) << GetName() << ": Running in the event by event mode."; - + } if (fPIDMode == 1) { FairRootManager* fManger = FairRootManager::Instance(); @@ -194,7 +198,9 @@ InitStatus CbmKFParticleFinderPID::Init() fDigiManager->Init(); // --- Check input array (StsDigis) - if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) LOG(fatal) << GetName() << ": No StsDigi branch in input!"; + if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) { + LOG(fatal) << GetName() << ": No StsDigi branch in input!"; + } // Get ToF hits fTofHitArray = (TClonesArray*) ioman->GetObject(fTofBranchName); @@ -243,8 +249,12 @@ void CbmKFParticleFinderPID::Exec(Option_t* /*opt*/) Int_t nTracks = fTrackArray->GetEntriesFast(); fPID.resize(nTracks, -1); - if (fPIDMode == 1) SetMCPID(); - if (fPIDMode == 2) SetRecoPID(); + if (fPIDMode == 1) { + SetMCPID(); + } + if (fPIDMode == 2) { + SetRecoPID(); + } } void CbmKFParticleFinderPID::Finish() {} @@ -253,13 +263,17 @@ void CbmKFParticleFinderPID::SetMCPID() { Int_t nTracks = fTrackArray->GetEntriesFast(); Int_t nMCTracks = 0; - if (fLegacyEventMode) nMCTracks = fMCTrackArray->GetEntriesFast(); + if (fLegacyEventMode) { + nMCTracks = fMCTrackArray->GetEntriesFast(); + } for (int iTr = 0; iTr < nTracks; iTr++) { fPID[iTr] = -2; CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fTrackMatchArray->At(iTr); - if (stsTrackMatch->GetNofLinks() == 0) continue; + if (stsTrackMatch->GetNofLinks() == 0) { + continue; + } Float_t bestWeight = 0.f; Float_t totalWeight = 0.f; Int_t mcTrackId = -1; @@ -277,9 +291,13 @@ void CbmKFParticleFinderPID::SetMCPID() } } } - if (bestWeight / totalWeight < 0.7) continue; + if (bestWeight / totalWeight < 0.7) { + continue; + } - if ((fLegacyEventMode) && (mcTrackId >= nMCTracks || mcTrackId < 0)) continue; + if ((fLegacyEventMode) && (mcTrackId >= nMCTracks || mcTrackId < 0)) { + continue; + } // if(mcTrackId >= nMCTracks || mcTrackId < 0) // { // LOG(info) << "Sts Matching is wrong! StsTrackId = " << mcTrackId << " N mc tracks = " << nMCTracks; @@ -302,10 +320,12 @@ void CbmKFParticleFinderPID::SetMCPID() || TMath::Abs(cbmMCTrack->GetPdgCode()) == 211 || TMath::Abs(cbmMCTrack->GetPdgCode()) == 321 || TMath::Abs(cbmMCTrack->GetPdgCode()) == 2212 || TMath::Abs(cbmMCTrack->GetPdgCode()) == 1000010020 || TMath::Abs(cbmMCTrack->GetPdgCode()) == 1000010030 || TMath::Abs(cbmMCTrack->GetPdgCode()) == 1000020030 - || TMath::Abs(cbmMCTrack->GetPdgCode()) == 1000020040)) + || TMath::Abs(cbmMCTrack->GetPdgCode()) == 1000020040)) { fPID[iTr] = -1; - else + } + else { fPID[iTr] = cbmMCTrack->GetPdgCode(); + } } } @@ -336,7 +356,9 @@ void CbmKFParticleFinderPID::SetRecoPID() const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTrackArray->At(igt)); Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); - if (stsTrackIndex < 0) continue; + if (stsTrackIndex < 0) { + continue; + } Bool_t isElectronTRD = 0; Bool_t isElectronRICH = 0; @@ -375,8 +397,9 @@ void CbmKFParticleFinderPID::SetRecoPID() // if(fElIdAnn->DoSelect(richRing, p) > -0.5) isElectronRICH = 1; if (p < 5.) { if (fabs(axisA - fMeanA) < fRmsCoeff * fRmsA && fabs(axisB - fMeanB) < fRmsCoeff * fRmsB - && dist < fDistCut) + && dist < fDistCut) { isElectronRICH = 1; + } } else { ///3 sigma @@ -388,8 +411,9 @@ void CbmKFParticleFinderPID::SetRecoPID() Double_t polBaxis = 5.41106 - 4.49902 / (p - 3.52450); //Double_t polRaxis = 5.66516 - 6.62229/(momentum - 2.25304); if (axisA < (fMeanA + fRmsCoeff * fRmsA) && axisA > polAaxis && axisB < (fMeanB + fRmsCoeff * fRmsB) - && axisB > polBaxis && dist < fDistCut) + && axisB > polBaxis && dist < fDistCut) { isElectronRICH = 1; + } } } } @@ -402,7 +426,9 @@ void CbmKFParticleFinderPID::SetRecoPID() if (trdIndex > -1) { CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTrackArray->At(trdIndex); if (trdTrack) { - if (trdTrack->GetPidWkn() > 0.635) isElectronTRD = 1; + if (trdTrack->GetPidWkn() > 0.635) { + isElectronTRD = 1; + } } } } @@ -414,7 +440,9 @@ void CbmKFParticleFinderPID::SetRecoPID() if (trdIndex > -1) { CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTrackArray->At(trdIndex); if (trdTrack) { - if (trdTrack->GetPidANN() > 0.9) isElectronTRD = 1; + if (trdTrack->GetPidANN() > 0.9) { + isElectronTRD = 1; + } } } } @@ -433,7 +461,9 @@ void CbmKFParticleFinderPID::SetRecoPID() CbmTrdHit* trdHit = (CbmTrdHit*) fTrdHitArray->At(TRDindex); eloss += trdHit->GetELoss(); } - if (trdTrack->GetNofHits() > 0.) dEdXTRD = 1e6 * pz / p * eloss / trdTrack->GetNofHits(); + if (trdTrack->GetNofHits() > 0.) { + dEdXTRD = 1e6 * pz / p * eloss / trdTrack->GetNofHits(); + } } } } @@ -470,20 +500,29 @@ void CbmKFParticleFinderPID::SetRecoPID() //check if at least one digi in a cluster has overflow --- charge is registered in the last ADC channel #31 for (int iDigi = 0; iDigi < frontCluster->GetNofDigis(); ++iDigi) { - if (31 == (fDigiManager->Get<CbmStsDigi>(frontCluster->GetDigi(iDigi))->GetCharge())) frontVeto = kTRUE; + if (31 == (fDigiManager->Get<CbmStsDigi>(frontCluster->GetDigi(iDigi))->GetCharge())) { + frontVeto = kTRUE; + } } for (int iDigi = 0; iDigi < backCluster->GetNofDigis(); ++iDigi) { - if (31 == (fDigiManager->Get<CbmStsDigi>(backCluster->GetDigi(iDigi))->GetCharge())) backVeto = kTRUE; + if (31 == (fDigiManager->Get<CbmStsDigi>(backCluster->GetDigi(iDigi))->GetCharge())) { + backVeto = kTRUE; + } } - if (!frontVeto) dEdxAllveto.push_back((frontCluster->GetCharge()) / dr); - if (!backVeto) dEdxAllveto.push_back((backCluster->GetCharge()) / dr); - - if (0 == frontVeto) nClustersWveto--; - if (0 == backVeto) nClustersWveto--; + if (!frontVeto) { + dEdxAllveto.push_back((frontCluster->GetCharge()) / dr); + nClustersWveto--; + } + if (!backVeto) { + dEdxAllveto.push_back((backCluster->GetCharge()) / dr); + nClustersWveto--; + } } float dEdXSTS = 1.e6; - if (dEdxAllveto.size() != 0) dEdXSTS = vecMedian(dEdxAllveto); + if (dEdxAllveto.size() != 0) { + dEdXSTS = vecMedian(dEdxAllveto); + } int isMuon = 0; @@ -495,27 +534,42 @@ void CbmKFParticleFinderPID::SetRecoPID() if ((cbmStsTrack->GetChiSq() / cbmStsTrack->GetNDF()) < fMuchCutsFloat[0] && (muchTrack->GetChiSq() / muchTrack->GetNDF()) < fMuchCutsFloat[1] && cbmStsTrack->GetTotalNofHits() >= fMuchCutsInt[0]) { - if (muchTrack->GetNofHits() >= fMuchCutsInt[1]) isMuon = 2; - if (muchTrack->GetNofHits() >= fMuchCutsInt[2]) isMuon = 1; + if (muchTrack->GetNofHits() >= fMuchCutsInt[1]) { + isMuon = 2; + } + if (muchTrack->GetNofHits() >= fMuchCutsInt[2]) { + isMuon = 1; + } } } } } if (p > 1.5) { - if (isElectronRICH && isElectronTRD) isElectron = 1; - if (fRichPIDMode == 0 && isElectronTRD) isElectron = 1; - if (fTrdPIDMode == 0 && isElectronRICH) isElectron = 1; + if (isElectronRICH && isElectronTRD) { + isElectron = 1; + } + if (fRichPIDMode == 0 && isElectronTRD) { + isElectron = 1; + } + if (fTrdPIDMode == 0 && isElectronRICH) { + isElectron = 1; + } } - else if (isElectronRICH) + else if (isElectronRICH) { isElectron = 1; + } if (fTofHitArray && isMuon == 0) { Double_t l = globalTrack->GetLength(); // l is calculated by global tracking - if (!((l > fCuts.fTrackLengthMin) && (l < fCuts.fTrackLengthMax))) continue; + if (!((l > fCuts.fTrackLengthMin) && (l < fCuts.fTrackLengthMax))) { + continue; + } Int_t tofHitIndex = globalTrack->GetTofHitIndex(); - if (tofHitIndex < 0) continue; + if (tofHitIndex < 0) { + continue; + } const CbmTofHit* tofHit = static_cast<const CbmTofHit*>(fTofHitArray->At(tofHitIndex)); if (!tofHit) { @@ -525,7 +579,9 @@ void CbmKFParticleFinderPID::SetRecoPID() Double_t time = tofHit->GetTime() - eventTime; - if (!((time > fCuts.fTrackTofTimeMin) && (time < fCuts.fTrackTofTimeMax))) continue; + if (!((time > fCuts.fTrackTofTimeMin) && (time < fCuts.fTrackTofTimeMax))) { + continue; + } Double_t m2 = p * p * (1. / ((l / time / 29.9792458) * (l / time / 29.9792458)) - 1.); @@ -539,7 +595,9 @@ void CbmKFParticleFinderPID::SetRecoPID() } if (isElectron) { - if (dm2[3] > 3.) isElectron = 0; + if (dm2[3] > 3.) { + isElectron = 0; + } } int iPdg = 2; @@ -549,30 +607,44 @@ void CbmKFParticleFinderPID::SetRecoPID() // if(p>12.) continue; for (int jPDG = 0; jPDG < 7; jPDG++) { - if (jPDG == 3) continue; + if (jPDG == 3) { + continue; + } if (dm2[jPDG] < dm2min) { iPdg = jPDG; dm2min = dm2[jPDG]; } } - if (dm2min > 2) iPdg = -1; + if (dm2min > 2) { + iPdg = -1; + } Double_t dEdXTRDthresholdProton = 10.; if (iPdg == 6) { - if (fUseTRDdEdX && (dEdXTRD < dEdXTRDthresholdProton)) iPdg = 0; - if (fUseSTSdEdX && (dEdXSTS < 7.5e4)) iPdg = 0; + if (fUseTRDdEdX && (dEdXTRD < dEdXTRDthresholdProton)) { + iPdg = 0; + } + if (fUseSTSdEdX && (dEdXSTS < 7.5e4)) { + iPdg = 0; + } } - if (iPdg > -1) fPID[stsTrackIndex] = q * PdgHypo[iPdg]; + if (iPdg > -1) { + fPID[stsTrackIndex] = q * PdgHypo[iPdg]; + } } } - if (isElectron) fPID[stsTrackIndex] = q * PdgHypo[3]; + if (isElectron) { + fPID[stsTrackIndex] = q * PdgHypo[3]; + } - if (isMuon == 1) fPID[stsTrackIndex] = q * PdgHypo[7]; - if (isMuon == 2) fPID[stsTrackIndex] = q * PdgHypo[8]; + if (isMuon == 1) { + fPID[stsTrackIndex] = q * PdgHypo[7]; + } + if (isMuon == 2) { + fPID[stsTrackIndex] = q * PdgHypo[8]; + } } } - -ClassImp(CbmKFParticleFinderPID); diff --git a/reco/KF/CbmKFParticleFinderQa.cxx b/reco/KF/CbmKFParticleFinderQa.cxx index ee5dbd5f82ff588839d8048bf9e1cbcd6ec778a9..73f938e9b65e6954c5b318b593ba6a32b4789854 100644 --- a/reco/KF/CbmKFParticleFinderQa.cxx +++ b/reco/KF/CbmKFParticleFinderQa.cxx @@ -45,8 +45,9 @@ CbmKFParticleFinderQa::CbmKFParticleFinderQa(const char* name, Int_t iVerbose, c : FairTask(name, iVerbose) , fOutFileName(outFileName) { - for (Int_t i = 0; i < 5; i++) + for (Int_t i = 0; i < 5; i++) { fTime[i] = 0; + } fTopoPerformance = new KFTopoPerformance; fTopoPerformance->SetTopoReconstructor(tr); @@ -54,8 +55,9 @@ CbmKFParticleFinderQa::CbmKFParticleFinderQa(const char* name, Int_t iVerbose, c TFile* curFile = gFile; TDirectory* curDirectory = gDirectory; - if (!(fOutFileName == "")) + if (!(fOutFileName == "")) { fOutFile = new TFile(fOutFileName.Data(), "RECREATE"); + } else { fOutFile = gFile; } @@ -67,9 +69,13 @@ CbmKFParticleFinderQa::CbmKFParticleFinderQa(const char* name, Int_t iVerbose, c CbmKFParticleFinderQa::~CbmKFParticleFinderQa() { - if (fTopoPerformance) delete fTopoPerformance; + if (fTopoPerformance) { + delete fTopoPerformance; + } - if (fSaveParticles) fRecParticles->Delete(); + if (fSaveParticles) { + fRecParticles->Delete(); + } if (fSaveMCParticles) { fMCParticles->Delete(); fMatchParticles->Delete(); @@ -96,7 +102,9 @@ InitStatus CbmKFParticleFinderQa::Init() if (!fLegacyEventMode) { FairRootManager* fManger = FairRootManager::Instance(); CbmMCDataManager* mcManager = (CbmMCDataManager*) fManger->GetObject("MCDataManager"); - if (mcManager == nullptr) Error("CbmKFParticleFinderQa::Init", "MC Data Manager not found!"); + if (mcManager == nullptr) { + Error("CbmKFParticleFinderQa::Init", "MC Data Manager not found!"); + } fMCTrackArray = mcManager->InitBranch("MCTrack"); @@ -153,7 +161,9 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) return; } - if (fSaveParticles) fRecParticles->Delete(); + if (fSaveParticles) { + fRecParticles->Delete(); + } if (fSaveMCParticles) { fMCParticles->Delete(); fMatchParticles->Delete(); @@ -214,32 +224,45 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) Int_t pdg = cbmMCTrack->GetPdgCode(); Double_t q = 1; - if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) + if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; - else if (pdg == 1000010020) + } + else if (pdg == 1000010020) { q = 1; - else if (pdg == -1000010020) + } + else if (pdg == -1000010020) { q = -1; - else if (pdg == 1000010030) + } + else if (pdg == 1000010030) { q = 1; - else if (pdg == -1000010030) + } + else if (pdg == -1000010030) { q = -1; - else if (pdg == 1000020030) + } + else if (pdg == 1000020030) { q = 2; - else if (pdg == -1000020030) + } + else if (pdg == -1000020030) { q = -2; - else if (pdg == 1000020040) + } + else if (pdg == 1000020040) { q = 2; - else if (pdg == -1000020040) + } + else if (pdg == -1000020040) { q = -2; - else + } + else { q = 0; + } + Double_t p = cbmMCTrack->GetP(); - if (cbmMCTrack->GetMotherId() >= 0) + if (cbmMCTrack->GetMotherId() >= 0) { kfMCTrack.SetMotherId(cbmMCTrack->GetMotherId() + mcIndexOffset); - else + } + else { kfMCTrack.SetMotherId(-iMCEvent - 1); + } kfMCTrack.SetQP(q / p); kfMCTrack.SetPDG(pdg); kfMCTrack.SetNMCPoints(0); @@ -260,7 +283,9 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fTrackMatchArray->At(iTr); - if (stsTrackMatch->GetNofLinks() == 0) continue; + if (stsTrackMatch->GetNofLinks() == 0) { + continue; + } Float_t bestWeight = 0.f; Float_t totalWeight = 0.f; Int_t mcTrackId = -1; @@ -284,9 +309,13 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) } } - if (mcTrackId < 0) continue; + if (mcTrackId < 0) { + continue; + } - if (bestWeight / totalWeight < 0.7) continue; + if (bestWeight / totalWeight < 0.7) { + continue; + } // if(mcTrackId >= nMCTracks || mcTrackId < 0) // { // std::cout << "Sts Matching is wrong! StsTackId = " << mcTrackId << " N mc tracks = " << nMCTracks << std::endl; @@ -296,8 +325,9 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) if (TMath::Abs(mcTracks[mcTrackId].PDG()) > 4000 && !(TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010020 || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010030 || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020030 - || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020040)) + || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020040)) { continue; + } mcTracks[mcTrackId].SetReconstructed(); trackMatch[iTr] = mcTrackId; @@ -312,8 +342,9 @@ void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) fNEvents++; fTime[4] += fTopoPerformance->GetTopoReconstructor()->Time(); - for (int iT = 0; iT < 4; iT++) + for (int iT = 0; iT < 4; iT++) { fTime[iT] += fTopoPerformance->GetTopoReconstructor()->StatTime(iT); + } if (fNEvents % fPrintFrequency == 0) { std::cout << "Topo reconstruction time" << " Real = " << std::setw(10) << fTime[4] / fNEvents * 1.e3 << " ms" << std::endl; @@ -375,8 +406,9 @@ void CbmKFParticleFinderQa::Finish() fTopoPerformance->FillHistos(); fTime[4] += fTopoPerformance->GetTopoReconstructor()->Time(); - for (int iT = 0; iT < 4; iT++) + for (int iT = 0; iT < 4; iT++) { fTime[iT] += fTopoPerformance->GetTopoReconstructor()->StatTime(iT); + } std::cout << "Topo reconstruction time" << " Real = " << std::setw(10) << fTime[4] * 1.e3 << " ms" << std::endl; @@ -394,10 +426,12 @@ void CbmKFParticleFinderQa::Finish() WriteHistosCurFile(fTopoPerformance->GetHistosDirectory()); if (fCheckDecayQA && fDecayToAnalyse > -1) { - if (fDecayToAnalyse < 0) + if (fDecayToAnalyse < 0) { Error("CbmKFParticleFinderQa::Finish", "Please specify the decay to be analysed."); - else + } + else { CheckDecayQA(); + } } if (!(fOutFileName == "")) { @@ -414,8 +448,9 @@ void CbmKFParticleFinderQa::Finish() void CbmKFParticleFinderQa::WriteHistosCurFile(TObject* obj) { - if (!obj->IsFolder()) + if (!obj->IsFolder()) { obj->Write(); + } else { TDirectory* cur = gDirectory; TFile* currentFile = gFile; @@ -487,8 +522,11 @@ void CbmKFParticleFinderQa::FitDecayQAHistograms(float sigma[14], const bool sav fitCanvas.SaveAs(canvasFile.Data()); } - for (int iHisto = 0; iHisto < nParameters * 2; iHisto++) - if (fit[iHisto]) delete fit[iHisto]; + for (int iHisto = 0; iHisto < nParameters * 2; iHisto++) { + if (fit[iHisto]) { + delete fit[iHisto]; + } + } } void CbmKFParticleFinderQa::CheckDecayQA() @@ -519,11 +557,13 @@ void CbmKFParticleFinderQa::CheckDecayQA() "#sigma_{p_{y}}", "#sigma_{p_{z}}", "#sigma_{E}", "P_{x}", "P_{y}", "P_{z}", "P_{p_{x}}", "P_{p_{y}}", "P_{p_{z}}", "P_{E}", "#varepsilon_{4#pi}", "#varepsilon_{KFP}"}; - for (int iBin = 0; iBin < 16; iBin++) + for (int iBin = 0; iBin < 16; iBin++) { qaHisto->GetXaxis()->SetBinLabel(iBin + 1, binLabel[iBin].Data()); + } - for (int iSigma = 0; iSigma < 14; iSigma++) + for (int iSigma = 0; iSigma < 14; iSigma++) { qaHisto->SetBinContent(iSigma + 1, sigma[iSigma]); + } qaHisto->SetBinContent(15, fTopoPerformance->fParteff.GetTotal4piEfficiency(fDecayToAnalyse)); qaHisto->SetBinContent(16, fTopoPerformance->fParteff.GetTotalKFPEfficiency(fDecayToAnalyse)); @@ -536,18 +576,21 @@ void CbmKFParticleFinderQa::CheckDecayQA() TH1F* referenceHisto = referenceFile->Get<TH1F>(qaHistoName); if (referenceHisto) { fTestOk = true; - for (int iBin = 1; iBin <= 7; iBin++) + for (int iBin = 1; iBin <= 7; iBin++) { fTestOk &= fabs(referenceHisto->GetBinContent(iBin) - qaHisto->GetBinContent(iBin)) / referenceHisto->GetBinContent(iBin) < 0.25; - for (int iBin = 8; iBin <= 14; iBin++) + } + for (int iBin = 8; iBin <= 14; iBin++) { fTestOk &= fabs(referenceHisto->GetBinContent(iBin) - qaHisto->GetBinContent(iBin)) / referenceHisto->GetBinContent(iBin) < 0.25; - for (int iBin = 15; iBin <= 16; iBin++) + } + for (int iBin = 15; iBin <= 16; iBin++) { fTestOk &= fabs(referenceHisto->GetBinContent(iBin) - qaHisto->GetBinContent(iBin)) / referenceHisto->GetBinContent(iBin) < 0.1; + } } referenceFile->Close(); referenceFile->Delete(); diff --git a/reco/KF/CbmKFParticleInterface.cxx b/reco/KF/CbmKFParticleInterface.cxx index 850a0e9bd56847cf5a7695112926450cfb83ab69..51de3abbec8ee94096dd6daf0d17b5802a9b94da 100644 --- a/reco/KF/CbmKFParticleInterface.cxx +++ b/reco/KF/CbmKFParticleInterface.cxx @@ -28,6 +28,7 @@ #include <vector> using std::vector; +ClassImp(CbmKFParticleInterface); void CbmKFParticleInterface::SetKFParticleFromStsTrack(CbmStsTrack* track, KFParticle* particle, Int_t pdg, Bool_t firstPoint) @@ -48,25 +49,35 @@ void CbmKFParticleInterface::SetKFParticleFromStsTrack(CbmStsTrack* track, KFPar for (Int_t iTr = 0; iTr < 1; iTr++) { const FairTrackParam* parameters; - if (firstPoint) + if (firstPoint) { parameters = vRTracks[iTr].GetParamFirst(); - else + } + else { parameters = vRTracks[iTr].GetParamLast(); + } float par[6] = {0.f}; Double_t V[15] = {0.f}; - for (Int_t i = 0, iCov = 0; i < 5; i++) - for (Int_t j = 0; j <= i; j++, iCov++) + for (Int_t i = 0, iCov = 0; i < 5; i++) { + for (Int_t j = 0; j <= i; j++, iCov++) { V[iCov] = parameters->GetCovariance(i, j); + } + } float a = parameters->GetTx(), b = parameters->GetTy(), qp = parameters->GetQp(); Int_t q = 0; - if (qp > 0.f) q = 1; - if (qp < 0.f) q = -1; - if (TMath::Abs(pdg) == 1000020030 || TMath::Abs(pdg) == 1000020040) q *= 2; + if (qp > 0.f) { + q = 1; + } + if (qp < 0.f) { + q = -1; + } + if (TMath::Abs(pdg) == 1000020030 || TMath::Abs(pdg) == 1000020040) { + q *= 2; + } float c2 = 1.f / (1.f + a * a + b * b); float pq = 1.f / qp * TMath::Abs(q); @@ -150,7 +161,9 @@ void CbmKFParticleInterface::ExtrapolateTrackToPV(const CbmStsTrack* track, CbmV vector<float> vChiToPrimVtx; CbmKFVertex kfVertex; assert(pv); - if (pv) kfVertex = CbmKFVertex(*pv); + if (pv) { + kfVertex = CbmKFVertex(*pv); + } vector<CbmL1PFFitter::PFFieldRegion> vField; fitter.GetChiToVertex(vRTracks, vField, vChiToPrimVtx, kfVertex, 1000000.f); @@ -158,5 +171,3 @@ void CbmKFParticleInterface::ExtrapolateTrackToPV(const CbmStsTrack* track, CbmV chiPrim = vChiToPrimVtx[0]; *paramAtPV = *(vRTracks[0].GetParamFirst()); } - -ClassImp(CbmKFParticleInterface); diff --git a/reco/KF/CbmKFPixelMeasurement.cxx b/reco/KF/CbmKFPixelMeasurement.cxx index d3e150c02423ba6df51bc9fc35531c9d0db1a3b5..825b327d591eb4c686daf2a5b15d52b6c8e92c94 100644 --- a/reco/KF/CbmKFPixelMeasurement.cxx +++ b/reco/KF/CbmKFPixelMeasurement.cxx @@ -25,7 +25,9 @@ Int_t CbmKFPixelMeasurement::Filter(CbmKFTrackInterface& track) return err; } w = 1. / w; - if (!isfinite(w)) return 1; + if (!isfinite(w)) { + return 1; + } Double_t W[3] = {w * (C[2] + V[2]), -w * (C[1] + V[1]), w * (C[0] + V[0])}; Double_t zeta[2] = {T[0] - x, T[1] - y}; @@ -133,7 +135,9 @@ void CbmKFPixelMeasurement::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFP Ck = 2.0 * 3.1415926 * sqrt(Sk) * (1.0 - Pdetect * gateEff) / (gateArea * Pdetect * gateEff) * Lambda; double w = Sk; - if (w < 1.E-20) return; + if (w < 1.E-20) { + return; + } w = 1. / w; // 2. Kalman gain - assumed to be the same for all gathered hits in a window @@ -212,22 +216,24 @@ void CbmKFPixelMeasurement::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFP idx++; } double KD[5][2]; - for (i = 0; i < 5; i++) + for (i = 0; i < 5; i++) { for (j = 0; j < 2; j++) { KD[i][j] = 0.0; for (k = 0; k < 2; k++) KD[i][j] += K[i][k] * CR[k][j]; } + } // 7. Additional covariance term - reflects influence of the residual spread double Ga[5][5]; - for (i = 0; i < 5; i++) + for (i = 0; i < 5; i++) { for (j = 0; j < 5; j++) { Ga[i][j] = 0.0; for (k = 0; k < 2; k++) Ga[i][j] += KD[i][k] * K[j][k]; } + } // 8. Track parameters update - usual Kalman formalizm but with the merged // residual @@ -249,8 +255,9 @@ void CbmKFPixelMeasurement::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFP C[6] * K[4][0] + C[7] * K[4][1], C[10] * K[4][0] + C[11] * K[4][1]}; double Gs[15]; - for (i = 0; i < 15; i++) + for (i = 0; i < 15; i++) { Gs[i] = C[i]; + } // 9. Standard Kalman covariance @@ -274,11 +281,12 @@ void CbmKFPixelMeasurement::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFP // covariance produced by the Kalman filter + additional term idx = 0; - for (i = 0; i < 5; i++) + for (i = 0; i < 5; i++) { for (j = 0; j <= i; j++) { C[idx] = P0 * C[idx] + (1 - P0) * Gs[idx] + Ga[i][j]; idx++; } + } #ifdef _DEBUG_PDAF_ cout << "Updated params: x=" << T[0] << " y=" << T[1] << " tx=" << T[2] << " ty=" << T[3] << " Q=" << T[4] << endl; cout << "chi2 contrib.=" << chi2 << endl; diff --git a/reco/KF/CbmKFPrimaryVertexFinder.cxx b/reco/KF/CbmKFPrimaryVertexFinder.cxx index 981a8f15856a342282a0aec031e73750b2171748..8015e993c98243c1a66de98d2dce026eab47e2ac 100644 --- a/reco/KF/CbmKFPrimaryVertexFinder.cxx +++ b/reco/KF/CbmKFPrimaryVertexFinder.cxx @@ -42,8 +42,10 @@ void CbmKFPrimaryVertexFinder::Fit(CbmKFVertexInterface& vtx) //* Initialize the vertex - for (Int_t i = 0; i < 6; i++) + for (Int_t i = 0; i < 6; i++) { C[i] = 0; + } + if (CbmKF::Instance()->vTargets.empty()) { r[0] = r[1] = r[2] = 0.; C[0] = C[2] = 5.; @@ -65,10 +67,12 @@ void CbmKFPrimaryVertexFinder::Fit(CbmKFVertexInterface& vtx) Double_t r0[3], C0[6]; - for (Int_t i = 0; i < 3; i++) + for (Int_t i = 0; i < 3; i++) { r0[i] = r[i]; - for (Int_t i = 0; i < 6; i++) + } + for (Int_t i = 0; i < 6; i++) { C0[i] = C[i]; + } //* Initialize the vertex covariance, Chi^2 & NDF @@ -87,7 +91,9 @@ void CbmKFPrimaryVertexFinder::Fit(CbmKFVertexInterface& vtx) CbmKFTrack T(**itr); Bool_t err = T.Extrapolate(r0[2]); - if (err) continue; + if (err) { + continue; + } Double_t* m = T.GetTrack(); Double_t* V = T.GetCovMatrix(); @@ -159,8 +165,9 @@ void CbmKFPrimaryVertexFinder::Fit(CbmKFVertexInterface& vtx) //* New estimation of the vertex position r += K*zeta - for (Int_t i = 0; i < 3; ++i) + for (Int_t i = 0; i < 3; ++i) { r[i] += K[i][0] * zeta[0] + K[i][1] * zeta[1]; + } //* New covariance matrix C -= K*(CH')' diff --git a/reco/KF/CbmKFSecondaryVertexFinder.cxx b/reco/KF/CbmKFSecondaryVertexFinder.cxx index bcccb2ca0a23183f192ed02f21c91a5845653b04..10dc46bc25bb657b2751229dd34cbb81812d49e7 100644 --- a/reco/KF/CbmKFSecondaryVertexFinder.cxx +++ b/reco/KF/CbmKFSecondaryVertexFinder.cxx @@ -70,16 +70,18 @@ void CbmKFSecondaryVertexFinder::Fit() for (Int_t iteration = 0; iteration < MaxIter; ++iteration) { - for (int i = 0; i < 8; i++) + for (int i = 0; i < 8; i++) { r0[i] = r[i]; + } r[3] = 0; r[4] = 0; r[5] = 0; r[6] = 0; - for (Int_t i = 0; i < 36; ++i) + for (Int_t i = 0; i < 36; ++i) { C[i] = 0.0; + } C[0] = C[2] = C[5] = 100.0; C[35] = 1.E4; @@ -208,7 +210,9 @@ void CbmKFSecondaryVertexFinder::Fit() { Double_t s = S[0] * S[2] - S[1] * S[1]; - if (s < 1.E-20) continue; + if (s < 1.E-20) { + continue; + } s = 1. / s; Double_t S0 = S[0]; S[0] = s * S[2]; @@ -246,14 +250,16 @@ void CbmKFSecondaryVertexFinder::Fit() //* New estimation of the vertex position r += K*zeta - for (Int_t i = 0; i < 7; ++i) + for (Int_t i = 0; i < 7; ++i) { r[i] += K0[i] * zeta[0] + K1[i] * zeta[1]; + } //* New covariance matrix C -= K*(CH')' for (Int_t i = 0, k = 0; i < 7; ++i) { - for (Int_t j = 0; j <= i; ++j, ++k) + for (Int_t j = 0; j <= i; ++j, ++k) { C[k] -= K0[i] * CHt0[j] + K1[i] * CHt1[j]; + } } //* Calculate Chi^2 & NDF @@ -277,8 +283,9 @@ void CbmKFSecondaryVertexFinder::GetVertex(CbmKFVertexInterface& vtx) vtx.GetRefX() = r[0]; vtx.GetRefY() = r[1]; vtx.GetRefZ() = r[2]; - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { vtx.GetCovMatrix()[i] = C[i]; + } vtx.GetRefChi2() = Chi2; vtx.GetRefNDF() = NDF; } @@ -294,7 +301,9 @@ void CbmKFSecondaryVertexFinder::GetVertex(CbmVertex& vtx) void CbmKFSecondaryVertexFinder::GetMotherTrack(Double_t Par[], Double_t Cov[]) { - if (!Par) return; + if (!Par) { + return; + } Double_t px = r[3], py = r[4], pz = r[5]; @@ -308,14 +317,18 @@ void CbmKFSecondaryVertexFinder::GetMotherTrack(Double_t Par[], Double_t Cov[]) Par[4] = qp; Par[5] = r[2]; - if (!Cov) return; + if (!Cov) { + return; + } double qp3 = qp * qp * qp; Double_t J[5][6]; - for (Int_t i = 0; i < 5; i++) - for (Int_t j = 0; j < 6; ++j) + for (Int_t i = 0; i < 5; i++) { + for (Int_t j = 0; j < 6; ++j) { J[i][j] = 0.0; + } + } J[0][0] = 1; @@ -332,20 +345,24 @@ void CbmKFSecondaryVertexFinder::GetMotherTrack(Double_t Par[], Double_t Cov[]) Double_t JC[5][6]; - for (Int_t i = 0; i < 5; i++) + for (Int_t i = 0; i < 5; i++) { for (Int_t j = 0; j < 6; j++) { JC[i][j] = 0; - for (Int_t k = 0; k < 6; k++) + for (Int_t k = 0; k < 6; k++) { JC[i][j] += J[i][k] * Cij(k, j); + } } + } Int_t ii = 0; - for (Int_t i = 0; i < 5; i++) + for (Int_t i = 0; i < 5; i++) { for (Int_t j = i; j < 5; j++, ii++) { Cov[ii] = 0; - for (Int_t k = 0; k < 6; k++) + for (Int_t k = 0; k < 6; k++) { Cov[ii] += JC[i][k] * J[j][k]; + } } + } } @@ -359,14 +376,20 @@ void CbmKFSecondaryVertexFinder::GetMass(Double_t* M, Double_t* Error_) - r[6] * (r[3] * C[24] + r[4] * C[25] + r[5] * C[26]))); Double_t m2 = r[6] * r[6] - r[3] * r[3] - r[4] * r[4] - r[5] * r[5]; - if (M) *M = (m2 > 1.e-20) ? sqrt(m2) : 0.; - if (Error_) *Error_ = (S >= 0 && m2 > 1.e-20) ? sqrt(S / m2) : 1.e4; + if (M) { + *M = (m2 > 1.e-20) ? sqrt(m2) : 0.; + } + if (Error_) { + *Error_ = (S >= 0 && m2 > 1.e-20) ? sqrt(S / m2) : 1.e4; + } } void CbmKFSecondaryVertexFinder::AddMassConstraint() { - if (MassConstraint < 0) return; + if (MassConstraint < 0) { + return; + } double H[8]; H[0] = H[1] = H[2] = 0.; H[3] = -2 * r0[3]; @@ -377,27 +400,32 @@ void CbmKFSecondaryVertexFinder::AddMassConstraint() double m2 = r0[6] * r0[6] - r0[3] * r0[3] - r0[4] * r0[4] - r0[5] * r0[5]; double zeta = MassConstraint * MassConstraint - m2; - for (Int_t i = 0; i < 8; ++i) + for (Int_t i = 0; i < 8; ++i) { zeta -= H[i] * (r[i] - r0[i]); + } Double_t S = 0.; Double_t CHt[8]; for (Int_t i = 0; i < 8; ++i) { CHt[i] = 0.0; - for (Int_t j = 0; j < 8; ++j) + for (Int_t j = 0; j < 8; ++j) { CHt[i] += Cij(i, j) * H[j]; + } S += H[i] * CHt[i]; } - if (S < 1.e-20) return; + if (S < 1.e-20) { + return; + } S = 1. / S; Chi2 += zeta * zeta * S; NDF += 1; for (Int_t i = 0, ii = 0; i < 8; ++i) { Double_t Ki = CHt[i] * S; r[i] += Ki * zeta; - for (Int_t j = 0; j <= i; ++j) + for (Int_t j = 0; j <= i; ++j) { C[ii++] -= Ki * CHt[j]; + } } } @@ -446,7 +474,9 @@ void CbmKFSecondaryVertexFinder::Extrapolate(double T) void CbmKFSecondaryVertexFinder::AddTopoConstraint() { - if (!VParent) return; + if (!VParent) { + return; + } double m[3], *V; { diff --git a/reco/KF/CbmKFTrackInterface.cxx b/reco/KF/CbmKFTrackInterface.cxx index 6f8da98fbaa987867a6b7defad0071dc68fb1497..81eebadd0ac13cea0c8d284c2805a8c9fa2bc26d 100644 --- a/reco/KF/CbmKFTrackInterface.cxx +++ b/reco/KF/CbmKFTrackInterface.cxx @@ -61,7 +61,9 @@ Int_t CbmKFTrackInterface::Extrapolate(Double_t z_out, Double_t* QP0) vector<CbmKFMaterial*>::iterator i, ibeg, iend; ibeg = lower_bound(KF->vMaterial.begin(), KF->vMaterial.end(), z, CbmKFMaterial::compareP_z); iend = upper_bound(KF->vMaterial.begin(), KF->vMaterial.end(), Z, CbmKFMaterial::compareP_Z); - if (iend != KF->vMaterial.end() && (*iend)->ZReference - (*iend)->ZThickness / 2 < Z) iend++; + if (iend != KF->vMaterial.end() && (*iend)->ZReference - (*iend)->ZThickness / 2 < Z) { + iend++; + } if (downstream_direction) { } else { @@ -74,7 +76,9 @@ Int_t CbmKFTrackInterface::Extrapolate(Double_t z_out, Double_t* QP0) for (i = ibeg; i != iend; downstream_direction ? ++i : --i) { Double_t zthick = (*i)->ZThickness, zcross = (*i)->ZReference; - if (CbmKFMath::GetThickness(z, Z, zcross, zthick, &zcross, &zthick)) continue; + if (CbmKFMath::GetThickness(z, Z, zcross, zthick, &zcross, &zthick)) { + continue; + } double z0 = (downstream_direction) ? zcross - zthick / 2. : zcross + zthick / 2.; double d = (downstream_direction) ? 1 : -1; @@ -89,7 +93,9 @@ Int_t CbmKFTrackInterface::Extrapolate(Double_t z_out, Double_t* QP0) //(*i)->Pass( zcross, zthick, *this, downstream_direction, qp0 ); } err = err || Propagate(z_out, qp0); - if (QP0) *QP0 = qp0; + if (QP0) { + *QP0 = qp0; + } return err; } @@ -102,7 +108,9 @@ Int_t CbmKFTrackInterface::Fit(Bool_t downstream) Int_t NHits = GetNOfHits(); Bool_t err = 0; - if (NHits == 0) return 1; + if (NHits == 0) { + return 1; + } // use initial momentum // this fixed value will be used during fit instead of T[4] @@ -141,8 +149,9 @@ Int_t CbmKFTrackInterface::Fit(Bool_t downstream) for (Int_t i = 1; i < NHits; i++) { h = GetHit(i); Int_t ist = h->MaterialIndex; - for (Int_t j = istold + 1; j < ist; j++) + for (Int_t j = istold + 1; j < ist; j++) { err = err || KF->vMaterial[j]->Pass(*this, downstream, qp0); + } err = err || h->Filter(*this, downstream, qp0); istold = ist; } @@ -154,8 +163,9 @@ Int_t CbmKFTrackInterface::Fit(Bool_t downstream) for (Int_t i = NHits - 2; i >= 0; i--) { h = GetHit(i); Int_t ist = h->MaterialIndex; - for (Int_t j = istold - 1; j > ist; j--) + for (Int_t j = istold - 1; j > ist; j--) { err = err || KF->vMaterial[j]->Pass(*this, downstream, qp0); + } err = err || h->Filter(*this, downstream, qp0); istold = ist; } @@ -218,7 +228,9 @@ void CbmKFTrackInterface::Smooth(Double_t Z) Double_t* C = GetCovMatrix(); Int_t NHits = GetNOfHits(); - if (NHits == 0) return; + if (NHits == 0) { + return; + } // use initial momentum // this fixed value will be used during fit instead of T[4] @@ -259,10 +271,14 @@ void CbmKFTrackInterface::Smooth(Double_t Z) Int_t ist = h->MaterialIndex; Int_t j = istold + 1; for (; j < ist; j++) { - if (KF->vMaterial[j]->ZReference > Z) break; + if (KF->vMaterial[j]->ZReference > Z) { + break; + } KF->vMaterial[j]->Pass(*this, 1, qp0); } - if (j < ist || KF->vMaterial[h->MaterialIndex]->ZReference > Z) break; + if (j < ist || KF->vMaterial[h->MaterialIndex]->ZReference > Z) { + break; + } h->Filter(*this, 1, qp0); istold = ist; } @@ -271,10 +287,12 @@ void CbmKFTrackInterface::Smooth(Double_t Z) KF->Propagate(T, C, Z, qp0); double Ts[6], Cs[15]; - for (int k = 0; k < 6; k++) + for (int k = 0; k < 6; k++) { Ts[k] = T[k]; - for (int k = 0; k < 15; k++) + } + for (int k = 0; k < 15; k++) { Cs[k] = C[k]; + } C[0] = INF; C[1] = 0.; C[2] = INF; @@ -301,10 +319,14 @@ void CbmKFTrackInterface::Smooth(Double_t Z) Int_t ist = h->MaterialIndex; Int_t j = istold - 1; for (; j > ist; j--) { - if (KF->vMaterial[j]->ZReference <= Z) break; + if (KF->vMaterial[j]->ZReference <= Z) { + break; + } KF->vMaterial[j]->Pass(*this, 0, qp0); } - if (j > ist || KF->vMaterial[h->MaterialIndex]->ZReference <= Z) break; + if (j > ist || KF->vMaterial[h->MaterialIndex]->ZReference <= Z) { + break; + } h->Filter(*this, 0, qp0); istold = ist; } @@ -315,31 +337,39 @@ void CbmKFTrackInterface::Smooth(Double_t Z) { // smooth double I[15]; - for (int k = 0; k < 15; k++) + for (int k = 0; k < 15; k++) { I[k] = C[k] + Cs[k]; + } CbmKFMath::invS(I, 5); double K[25]; CbmKFMath::multSSQ(C, I, K, 5); double r[5]; - for (int k = 0; k < 5; k++) + for (int k = 0; k < 5; k++) { r[k] = T[k] - Ts[k]; - for (int k = 0; k < 5; k++) - for (int l = 0; l < 5; l++) + } + for (int k = 0; k < 5; k++) { + for (int l = 0; l < 5; l++) { T[k] -= K[k * 5 + l] * r[l]; + } + } double A[15]; - for (int l = 0; l < 5; l++) + for (int l = 0; l < 5; l++) { K[(5 + 1) * l] -= 1; - for (int l = 0; l < 5; ++l) + } + for (int l = 0; l < 5; ++l) { for (int j = 0; j <= l; ++j) { int ind = CbmKFMath::indexS(l, j); A[ind] = 0; - for (int k = 0; k < 5; ++k) + for (int k = 0; k < 5; ++k) { A[ind] -= K[l * 5 + k] * C[CbmKFMath::indexS(k, j)]; + } } - for (int l = 0; l < 15; l++) + } + for (int l = 0; l < 15; l++) { C[l] = A[l]; + } } // correct NDF value to number of fitted track parameters (straight line(4) or helix(5) ) @@ -381,7 +411,9 @@ void CbmKFTrackInterface::Fit2Vertex(CbmKFVertexInterface& vtx) //* Invert S { Double_t s = S[0] * S[2] - S[1] * S[1]; - if (s < 1.E-20) return; + if (s < 1.E-20) { + return; + } s = 1. / s; Double_t S0 = S[0]; S[0] = s * S[2]; diff --git a/reco/KF/CbmKFUMeasurement.cxx b/reco/KF/CbmKFUMeasurement.cxx index b7938711dbb152091858a59ecd02e904b09128e4..f212161d2497a61fd4fab845eaf3bcd6c8af99fc 100644 --- a/reco/KF/CbmKFUMeasurement.cxx +++ b/reco/KF/CbmKFUMeasurement.cxx @@ -30,9 +30,13 @@ Int_t CbmKFUMeasurement::Filter(CbmKFTrackInterface& track) Double_t* C = track.GetCovMatrix(); Double_t W = sigma2 + phi_cc * C[0] + phi_2sc * C[1] + phi_ss * C[2]; - if (!isfinite(W) || W < 1.e-10) return 1; + if (!isfinite(W) || W < 1.e-10) { + return 1; + } W = 1. / W; - if (!isfinite(W)) return 1; + if (!isfinite(W)) { + return 1; + } Double_t zeta = phi_c * T[0] + phi_s * T[1] - u; diff --git a/reco/KF/CbmKFVertexInterface.cxx b/reco/KF/CbmKFVertexInterface.cxx index 28dc02835f1e6c8e54954ddc2880aecfc7f2c0f6..65554d3935176c3575952f051d6c299d40ec3691 100644 --- a/reco/KF/CbmKFVertexInterface.cxx +++ b/reco/KF/CbmKFVertexInterface.cxx @@ -45,16 +45,20 @@ void CbmKFVertexInterface::SetVertex(CbmVertex& v) GetRefNTracks() = v.GetNTracks(); TMatrixFSym tmp(3); v.CovMatrix(tmp); - for (int i = 0, k = 0; i < 3; i++) - for (int j = 0; j <= i; j++, k++) + for (int i = 0, k = 0; i < 3; i++) { + for (int j = 0; j <= i; j++, k++) { GetCovMatrix()[k] = tmp(i, j); + } + } } void CbmKFVertexInterface::GetVertex(CbmVertex& v) { TMatrixFSym covMat(3); - for (int i = 0, k = 0; i < 3; i++) - for (int j = 0; j <= i; j++, k++) + for (int i = 0, k = 0; i < 3; i++) { + for (int j = 0; j <= i; j++, k++) { covMat(i, j) = GetCovMatrix()[k]; + } + } v.SetVertex(GetRefX(), GetRefY(), GetRefZ(), GetRefChi2(), GetRefNDF(), GetRefNTracks(), covMat); } diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx index 827f6cc87c7be8edfba9ec72fbe7ac6d086f08d3..770826f5fa3c47dbd9ff9e53f21ffc4d9213d7c1 100644 --- a/reco/KF/CbmKfTrackFitter.cxx +++ b/reco/KF/CbmKfTrackFitter.cxx @@ -67,7 +67,9 @@ CbmKfTrackFitter::~CbmKfTrackFitter() {} void CbmKfTrackFitter::Init() { - if (fIsInitialized) return; + if (fIsInitialized) { + return; + } if (!CbmL1::Instance() || !CbmL1::Instance()->fpAlgo) { LOG(fatal) << "CbmKfTrackFitter: no CbmL1 task initialized "; @@ -123,7 +125,9 @@ void CbmKfTrackFitter::SetElectronFlag(bool isElectron) { fIsElectron = isElectr bool CbmKfTrackFitter::CreateMvdStsTrack(CbmKfTrackFitter::Track& kfTrack, int stsTrackIndex) { Init(); - if (!fIsInitialized) return false; + if (!fIsInitialized) { + return false; + } if (!fInputStsTracks) { LOG(error) << "CbmKfTrackFitter: Sts track array not found!"; @@ -150,7 +154,9 @@ bool CbmKfTrackFitter::CreateMvdStsTrack(CbmKfTrackFitter::Track& kfTrack, int s bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, int globalTrackIndex) { Init(); - if (!fIsInitialized) return false; + if (!fIsInitialized) { + return false; + } if (!fInputGlobalTracks) { LOG(error) << "CbmKfTrackFitter: Global track array not found!"; @@ -176,7 +182,9 @@ bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, const { kfTrack = {}; Init(); - if (!fIsInitialized) return false; + if (!fIsInitialized) { + return false; + } const ca::Parameters& caPar = CbmL1::Instance()->fpAlgo->GetParameters(); @@ -207,7 +215,9 @@ bool CbmKfTrackFitter::CreateGlobalTrack(CbmKfTrackFitter::Track& kfTrack, const // lambda to set the node from the pixel hit auto setNode = [&](const CbmPixelHit& h, int stIdx, int hitIdx, ca::EDetectorID detId) { int ista = caPar.GetStationIndexActive(stIdx, detId); - if (ista < 0) return; + if (ista < 0) { + return; + } assert(ista < nStations); CbmKfTrackFitter::FitNode& n = t.fNodes[ista]; n.fZ = h.GetZ(); diff --git a/reco/KF/Interface/CbmEcalTrackExtrapolationKF.cxx b/reco/KF/Interface/CbmEcalTrackExtrapolationKF.cxx index 21ea5d58310cea3e99e7abf04c0beaf06d54c3aa..7b1c7a86ec41df7ec7a5322d4fc62d59e912ef93 100644 --- a/reco/KF/Interface/CbmEcalTrackExtrapolationKF.cxx +++ b/reco/KF/Interface/CbmEcalTrackExtrapolationKF.cxx @@ -54,11 +54,15 @@ void CbmEcalTrackExtrapolationKF::Init() CbmEcalTrackExtrapolation::Init(); //Get and check FairRootManager FairRootManager* ioman = FairRootManager::Instance(); - if (!ioman) Fatal("CbmEcalTrackExtrapolationKF::Init()", "Can't instantise Rootmanager"); + if (!ioman) { + Fatal("CbmEcalTrackExtrapolationKF::Init()", "Can't instantise Rootmanager"); + } //Get STS track array fSTSArray = (TClonesArray*) ioman->GetObject("StsTrack"); - if (!fSTSArray) cout << "-W- CbmEcalTrackExtrapolationKF::Init: No sts track array!" << endl; + if (!fSTSArray) { + cout << "-W- CbmEcalTrackExtrapolationKF::Init: No sts track array!" << endl; + } } // ------------------------------------------------------------------------- @@ -86,8 +90,12 @@ Int_t CbmEcalTrackExtrapolationKF::DoExtrapolate(TClonesArray* gTrackArray, TClo for (; i < n; i++) { Map()[nTr] = -1111; tr = (CbmGlobalTrack*) gTrackArray->At(i); - if (!tr) continue; - if (tr->GetStsTrackIndex() < 0 || tr->GetTrdTrackIndex() < 0) continue; + if (!tr) { + continue; + } + if (tr->GetStsTrackIndex() < 0 || tr->GetTrdTrackIndex() < 0) { + continue; + } kfTr.SetTrackParam(*(const_cast<FairTrackParam*>(tr->GetParamLast()))); kfTr.Extrapolate(Str()->GetEcalInf()->GetZPos()); kfTr.GetTrackParam(trpar); diff --git a/reco/KF/Interface/CbmGlobalTrackFitterKF.cxx b/reco/KF/Interface/CbmGlobalTrackFitterKF.cxx index 5a6f48ecdc938ed64830f5d7afe2e327c16fda84..2c1efd86a26b665cc75794569a7bc20ffa74203d 100644 --- a/reco/KF/Interface/CbmGlobalTrackFitterKF.cxx +++ b/reco/KF/Interface/CbmGlobalTrackFitterKF.cxx @@ -126,8 +126,9 @@ void CbmGlobalTrackFitterKF::DoFit(CbmGlobalTrack* glbTrack) { // Implementation of the fitting algorithm if (nullptr == glbTrack || nullptr == fArrayStsTrack || nullptr == fArrayTrdTrack || nullptr == fArrayStsHit - || nullptr == fArrayTrdHit || nullptr == fPrimVertex) + || nullptr == fArrayTrdHit || nullptr == fPrimVertex) { return; + } Double_t x_old; diff --git a/reco/KF/Interface/CbmKFStsHit.cxx b/reco/KF/Interface/CbmKFStsHit.cxx index b76ac6c3d0f4d3893d3b554c1cc54fedd2938394..ab2b2461d2d5e29b4880b04c8656a57a23d422a3 100644 --- a/reco/KF/Interface/CbmKFStsHit.cxx +++ b/reco/KF/Interface/CbmKFStsHit.cxx @@ -27,8 +27,9 @@ void CbmKFStsHit::Create(CbmStsHit* h) MaterialIndex = KF->GetMaterialIndex(id); - if (MaterialIndex >= 0) + if (MaterialIndex >= 0) { tube = (CbmKFTube*) KF->vMaterial[MaterialIndex]; + } else { st_tube.z = st_tube.dz = st_tube.r = st_tube.R = st_tube.rr = st_tube.RR = 0; tube = &st_tube; @@ -61,8 +62,9 @@ void CbmKFStsHit::Create(CbmMvdHit* h) MaterialIndex = KF->GetMaterialIndex(id); // cout << " and material index = " << MaterialIndex << endl; - if (MaterialIndex >= 0) + if (MaterialIndex >= 0) { tube = (CbmKFTube*) KF->vMaterial[MaterialIndex]; + } else { st_tube.z = st_tube.dz = st_tube.r = st_tube.R = st_tube.rr = st_tube.RR = 0; tube = &st_tube; @@ -121,7 +123,9 @@ void CbmKFStsHit::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFStsHit*>& v { best_hit_idx = 0; - if (vpHits.empty()) return; + if (vpHits.empty()) { + return; + } double qp0 = (QP0) ? *QP0 : track.GetTrack()[4]; @@ -171,5 +175,7 @@ void CbmKFStsHit::FilterPDAF(CbmKFTrackInterface& track, vector<CbmKFStsHit*>& v tube->Pass(zlst, zthick, track, downstream, qp0); track.Propagate(zend, qp0); - if (QP0) *QP0 = qp0; + if (QP0) { + *QP0 = qp0; + } } diff --git a/reco/KF/Interface/CbmKFTrack.cxx b/reco/KF/Interface/CbmKFTrack.cxx index cf42e24e430d8bdb213f61a5bf44edfaf84220b2..8bcba6dabf12847ab34abf15bb9ff6834cdf4bda 100644 --- a/reco/KF/Interface/CbmKFTrack.cxx +++ b/reco/KF/Interface/CbmKFTrack.cxx @@ -20,19 +20,23 @@ ClassImp(CbmKFTrack) , fNDF(0) , fHits() { - for (Int_t i = 0; i < 6; i++) + for (Int_t i = 0; i < 6; i++) { fT[i] = 0.; - for (Int_t i = 0; i < 15; i++) + } + for (Int_t i = 0; i < 15; i++) { fC[i] = 0.; + } } void CbmKFTrack::SetTrack(CbmKFTrackInterface& track) { - for (Int_t i = 0; i < 6; i++) + for (Int_t i = 0; i < 6; i++) { fT[i] = track.GetTrack()[i]; - for (Int_t i = 0; i < 15; i++) + } + for (Int_t i = 0; i < 15; i++) { fC[i] = track.GetCovMatrix()[i]; + } fMass = track.GetMass(); fIsElectron = track.IsElectron(); fChi2 = track.GetRefChi2(); diff --git a/reco/KF/Interface/CbmPVFinderKF.cxx b/reco/KF/Interface/CbmPVFinderKF.cxx index 48a75b5b80467fc87bad1f21792578f7f84a20ff..958c8f7b71ae3b741da23f666e8f4dfdfe366e2c 100644 --- a/reco/KF/Interface/CbmPVFinderKF.cxx +++ b/reco/KF/Interface/CbmPVFinderKF.cxx @@ -25,12 +25,20 @@ ClassImp(CbmPVFinderKF) for (Int_t i = 0; i < NTracks; i++) { CbmStsTrack* st = (CbmStsTrack*) tracks->At(i); Int_t NHits = st->GetTotalNofHits(); - if (NHits < 4) continue; - if (st->GetFlag()) continue; - if (st->GetChiSq() < 0. || st->GetChiSq() > 3.5 * 3.5 * st->GetNDF()) continue; + if (NHits < 4) { + continue; + } + if (st->GetFlag()) { + continue; + } + if (st->GetChiSq() < 0. || st->GetChiSq() > 3.5 * 3.5 * st->GetNDF()) { + continue; + } CbmKFTrack& T = CloneArray[i]; T.SetStsTrack(*st); - if (!isfinite(T.GetTrack()[0]) || !isfinite(T.GetCovMatrix()[0])) continue; + if (!isfinite(T.GetTrack()[0]) || !isfinite(T.GetCovMatrix()[0])) { + continue; + } Finder.AddTrack(&T); } CbmKFVertex v; @@ -53,21 +61,29 @@ Int_t CbmPVFinderKF::FindEventVertex(CbmEvent* event, TClonesArray* tracks) // Copy input tracks to KF tracks Int_t nTracks = event->GetNofData(ECbmDataType::kStsTrack); - if (nTracks <= 0) return 0; + if (nTracks <= 0) { + return 0; + } CbmKFTrack* trackArray = new CbmKFTrack[nTracks]; for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { Int_t trackIndex = event->GetIndex(ECbmDataType::kStsTrack, iTrack); CbmStsTrack* track = (CbmStsTrack*) tracks->At(trackIndex); Int_t nHits = track->GetTotalNofHits(); - if (nHits < 4) continue; // use only tracks with at least 4 hits - if (track->GetFlag()) continue; // do not use suspicious tracks - if (track->GetChiSq() < 0. // use only good-quality tracks - || track->GetChiSq() > 12.25 * track->GetNDF()) + if (nHits < 4) { // use only tracks with at least 4 hits + continue; + } + if (track->GetFlag()) { // do not use suspicious tracks continue; + } + if (track->GetChiSq() < 0. || track->GetChiSq() > 12.25 * track->GetNDF()) { // use only good-quality tracks + continue; + } CbmKFTrack& kTrack = trackArray[iTrack]; kTrack.SetStsTrack(*track); - if (!isfinite(kTrack.GetTrack()[0]) || !isfinite(kTrack.GetCovMatrix()[0])) continue; + if (!isfinite(kTrack.GetTrack()[0]) || !isfinite(kTrack.GetCovMatrix()[0])) { + continue; + } vFinder.AddTrack(&kTrack); } diff --git a/reco/KF/Interface/CbmStsFitPerformanceTask.cxx b/reco/KF/Interface/CbmStsFitPerformanceTask.cxx index 1e9b07df2b4ae53871f8c432b1f2258b563fe378..03f494c8a461613aed44bd84eb402b548a3daadd 100644 --- a/reco/KF/Interface/CbmStsFitPerformanceTask.cxx +++ b/reco/KF/Interface/CbmStsFitPerformanceTask.cxx @@ -43,20 +43,24 @@ using std::endl; using std::fabs; using std::vector; +ClassImp(CbmStsFitPerformanceTask); + // ------------------------------------------------------------------------- void writedir2current(TObject* obj) { - if (!obj->IsFolder()) + if (!obj->IsFolder()) { obj->Write(); + } else { TDirectory* cur = gDirectory; TDirectory* sub = cur->mkdir(obj->GetName()); sub->cd(); TList* listSub = ((TDirectory*) obj)->GetList(); TIter it(listSub); - while (TObject* obj1 = it()) + while (TObject* obj1 = it()) { writedir2current(obj1); + } cur->cd(); } } @@ -395,7 +399,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) Int_t nTracks = fRecStsTrackArray->GetEntriesFast(); cout << " nTracks: " << nTracks << endl; - if (!fSTSTrackMatch) return; + if (!fSTSTrackMatch) { + return; + } Double_t mc[6]; @@ -404,19 +410,35 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { Quality[iTrack] = 0; CbmStsTrack* track = (CbmStsTrack*) fRecStsTrackArray->At(iTrack); - if (!track) continue; + if (!track) { + continue; + } CbmTrackMatch* m = (CbmTrackMatch*) fSTSTrackMatch->At(iTrack); - if (!m) continue; - if (m->GetNofTrueHits() < 0.7 * (m->GetNofTrueHits() + m->GetNofWrongHits() + m->GetNofFakeHits())) continue; + if (!m) { + continue; + } + if (m->GetNofTrueHits() < 0.7 * (m->GetNofTrueHits() + m->GetNofWrongHits() + m->GetNofFakeHits())) { + continue; + } Int_t mcTrackID = m->GetMCTrackId(); - if (mcTrackID < 0) continue; - if (!IsLong(track)) continue; + if (mcTrackID < 0) { + continue; + } + if (!IsLong(track)) { + continue; + } CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTrackArray->At(mcTrackID); - if (!mcTrack) continue; - if (fabs((mcTrack->GetStartZ())) > 1.) continue; + if (!mcTrack) { + continue; + } + if (fabs((mcTrack->GetStartZ())) > 1.) { + continue; + } - if (mcTrack->GetMotherId() != -1) continue; + if (mcTrack->GetMotherId() != -1) { + continue; + } Quality[iTrack] = 1; } @@ -426,9 +448,13 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) int nMvdHits = fMvdHitArray->GetEntriesFast(); cout << "Mvd hit density..." << endl; for (int ih = 0; ih < nMvdHits; ih++) { - if (ih % 100 == 0) cout << ih << endl; + if (ih % 100 == 0) { + cout << ih << endl; + } CbmMvdHit* h1 = (CbmMvdHit*) fMvdHitArray->At(ih); - if (!h1) continue; + if (!h1) { + continue; + } double V[3] = {2 * h1->GetDx() * h1->GetDx(), 2 * 0, 2 * h1->GetDy() * h1->GetDy()}; Double_t v = 1. / (V[0] * V[2] - V[1] * V[1]); V[0] *= v; @@ -436,22 +462,34 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) V[2] *= v; double D2 = 1.e10; for (int jh = 0; jh < nMvdHits; jh++) { - if (ih == jh) continue; + if (ih == jh) { + continue; + } CbmMvdHit* h2 = (CbmMvdHit*) fMvdHitArray->At(jh); - if (!h2) continue; - if (h1->GetStationNr() != h2->GetStationNr()) continue; + if (!h2) { + continue; + } + if (h1->GetStationNr() != h2->GetStationNr()) { + continue; + } Double_t dx = h1->GetX() - h2->GetX(); Double_t dy = h1->GetY() - h2->GetY(); Double_t d2 = fabs(dx * dx * V[0] - 2 * dx * dy * V[1] + dy * dy * V[2]); - if (d2 < D2) D2 = d2; + if (d2 < D2) { + D2 = d2; + } } fhHitDensity[h1->GetStationNr() - 1]->Fill(sqrt(D2 / 2.)); } cout << "Sts hit density..." << endl; for (int ih = 0; ih < nStsHits; ih++) { - if (ih % 1000 == 0) cout << ih << endl; + if (ih % 1000 == 0) { + cout << ih << endl; + } CbmStsHit* h1 = (CbmStsHit*) fStsHitArray->At(ih); - if (!h1) continue; + if (!h1) { + continue; + } double V[3] = {2 * h1->GetDx() * h1->GetDx(), 2 * h1->GetDxy(), 2 * h1->GetDy() * h1->GetDy()}; Double_t v = 1. / (V[0] * V[2] - V[1] * V[1]); V[0] *= v; @@ -459,12 +497,17 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) V[2] *= v; double D2 = 1.e10; for (int jh = 0; jh < nStsHits; jh++) { - if (ih == jh) continue; + if (ih == jh) { + continue; + } CbmStsHit* h2 = (CbmStsHit*) fStsHitArray->At(jh); - if (!h2) continue; + if (!h2) { + continue; + } if (CbmStsSetup::Instance()->GetStationNumber(h1->GetAddress()) - != CbmStsSetup::Instance()->GetStationNumber(h2->GetAddress())) + != CbmStsSetup::Instance()->GetStationNumber(h2->GetAddress())) { continue; + } // if( h1->GetBackDigiId()>=0 ){ // if( h1->GetFrontDigiId()!=h2->GetFrontDigiId() && // h1->GetBackDigiId()!=h2->GetBackDigiId() ) continue; @@ -472,7 +515,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) Double_t dx = h1->GetX() - h2->GetX(); Double_t dy = h1->GetY() - h2->GetY(); Double_t d2 = fabs(dx * dx * V[0] - 2 * dx * dy * V[1] + dy * dy * V[2]); - if (d2 < D2) D2 = d2; + if (d2 < D2) { + D2 = d2; + } } fhHitDensity[CbmKF::Instance()->MvdStationIDMap.size() + CbmStsSetup::Instance()->GetStationNumber(h1->GetAddress()) - 1] @@ -489,7 +534,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) flag[iTrack] = 0; CbmStsTrack* trackI = (CbmStsTrack*) fRecStsTrackArray->At(iTrack); - if (!IsLong(trackI)) continue; + if (!IsLong(trackI)) { + continue; + } flag[iTrack] = 1; double z[8] = {0.2, 1, 2, 3, 4, 5, 10, 100}; for (int iz = 0; iz < 8; iz++) { @@ -504,45 +551,62 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) } for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - if (iTrack % 10 == 0) cout << iTrack << endl; - if (!flag[iTrack]) continue; + if (iTrack % 10 == 0) { + cout << iTrack << endl; + } + if (!flag[iTrack]) { + continue; + } double Chi2[8]; - for (int iz = 0; iz < 8; iz++) + for (int iz = 0; iz < 8; iz++) { Chi2[iz] = 1.e10; + } double Chi2L = 1.e10; for (Int_t jTrack = 0; jTrack < nTracks; jTrack++) { - if (jTrack == iTrack) continue; - if (!flag[jTrack]) continue; + if (jTrack == iTrack) { + continue; + } + if (!flag[jTrack]) { + continue; + } for (int iz = 0; iz < 8; iz++) { Double_t C[15], d[5]; - for (int i = 0; i < 15; i++) + for (int i = 0; i < 15; i++) { C[i] = sC[iTrack][iz][i] + sC[jTrack][iz][i]; - for (int i = 0; i < 5; i++) + } + for (int i = 0; i < 5; i++) { d[i] = sT[iTrack][iz][i] - sT[jTrack][iz][i]; + } CbmKFMath::invS(C, 5); double chi2 = 0; for (int i = 0; i < 5; i++) { double dCi = 0; - for (int j = 0; j < 5; j++) + for (int j = 0; j < 5; j++) { dCi += d[j] * C[CbmKFMath::indexS(i, j)]; + } chi2 += dCi * d[i]; } chi2 = fabs(chi2); - if (chi2 < Chi2[iz]) Chi2[iz] = chi2; + if (chi2 < Chi2[iz]) { + Chi2[iz] = chi2; + } if (iz == 0) { chi2 = 0; for (int i = 0; i < 4; i++) { double dCi = 0; - for (int j = 0; j < 4; j++) + for (int j = 0; j < 4; j++) { dCi += d[j] * C[CbmKFMath::indexS(i, j)]; + } chi2 += dCi * d[i]; } chi2 = fabs(chi2); - if (chi2 < Chi2L) Chi2L = chi2; + if (chi2 < Chi2L) { + Chi2L = chi2; + } } } } @@ -554,42 +618,64 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) } if (fTrackAnalysis) { for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - if (Quality[iTrack] < 1) continue; + if (Quality[iTrack] < 1) { + continue; + } CbmStsTrack* track = (CbmStsTrack*) fRecStsTrackArray->At(iTrack); CbmTrackMatch* m = (CbmTrackMatch*) fSTSTrackMatch->At(iTrack); CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTrackArray->At(m->GetMCTrackId()); - if (mcTrack->GetMotherId() != -1) continue; // not primary track //IK + if (mcTrack->GetMotherId() != -1) { + continue; // not primary track //IK + } // Get MC points; vector<CbmStsPoint*> vPoints; for (Int_t i = 0; i < track->GetNofStsHits(); i++) { Int_t hitID = track->GetStsHitIndex(i); - if (hitID < 0) continue; + if (hitID < 0) { + continue; + } CbmStsHit* hit = (CbmStsHit*) fStsHitArray->At(hitID); - if (!hit) continue; + if (!hit) { + continue; + } Int_t pointID = hit->GetRefId(); - if (pointID < 0) continue; + if (pointID < 0) { + continue; + } CbmStsPoint* point = (CbmStsPoint*) fStsPointArray->At(pointID); - if (!point) continue; + if (!point) { + continue; + } vPoints.push_back(point); } vector<CbmMvdPoint*> vMPoints; for (Int_t i = 0; i < track->GetNofMvdHits(); i++) { Int_t hitID = track->GetMvdHitIndex(i); - if (hitID < 0) continue; + if (hitID < 0) { + continue; + } CbmMvdHit* hit = (CbmMvdHit*) fMvdHitArray->At(hitID); - if (!hit) continue; + if (!hit) { + continue; + } Int_t pointID = hit->GetRefIndex(); - if (pointID < 0) continue; + if (pointID < 0) { + continue; + } CbmMvdPoint* point = (CbmMvdPoint*) fMvdPointArray->At(pointID); - if (!point) continue; + if (!point) { + continue; + } vMPoints.push_back(point); } // impact mc track Double_t mci[6]; Double_t pzi = mcTrack->GetPz(); - if (fabs(pzi) > 1.e-4) pzi = 1. / pzi; + if (fabs(pzi) > 1.e-4) { + pzi = 1. / pzi; + } mci[0] = mcTrack->GetStartX(); mci[1] = mcTrack->GetStartY(); mci[2] = mcTrack->GetPx() * pzi; @@ -610,10 +696,12 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) fhDsP2->Fill(p1, dp / s); } // first point - if (!vMPoints.empty()) + if (!vMPoints.empty()) { FillMCStateVectors(vMPoints.front(), mc); - else + } + else { FillMCStateVectors(vPoints.front(), mc); + } FillTrackHisto(mc, track, fhFrst); fhZMCf->Fill(mc[5]); fhZRecof->Fill(track->GetParamFirst()->GetZ()); @@ -626,18 +714,24 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) // impact point - if (TMath::Abs(mcTrack->GetPdgCode()) == 211) + if (TMath::Abs(mcTrack->GetPdgCode()) == 211) { FillTrackHisto(mci, track, fhImpPi); - else + } + else { FillTrackHisto(mci, track, fhImp); + } // 4th station (check smoother) for (vector<CbmStsPoint*>::iterator i = vPoints.begin(); i != vPoints.end(); i++) { int id = (*i)->GetDetectorID(); CbmKF& KF = *CbmKF::Instance(); - if (KF.StsStationIDMap.find(id) == KF.StsStationIDMap.end()) continue; - if (KF.StsStationIDMap[id] != 3) continue; + if (KF.StsStationIDMap.find(id) == KF.StsStationIDMap.end()) { + continue; + } + if (KF.StsStationIDMap[id] != 3) { + continue; + } FillMCStateVectors(*i, mc, 1); //FillTrackHisto( mc, track, fhMid ); break; @@ -659,8 +753,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) double z = mci[5]; FairTrackParam param; fFitter.Extrapolate(track, z, ¶m); - if (fabs(mci[4]) > 1.e-4 && fabs(param.GetQp()) > 1.e-4) + if (fabs(mci[4]) > 1.e-4 && fabs(param.GetQp()) > 1.e-4) { fhPq->Fill(fabs(1. / mci[4]), (mci[4] / param.GetQp() - 1.)); + } } if (track->GetNDF() > 0) { @@ -682,7 +777,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) Int_t nMCTracks = fMCTrackArray->GetEntriesFast(); for (Int_t iTrack = 0; iTrack < nMCTracks; iTrack++) { CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTrackArray->At(iTrack); - if (mcTrack->GetMotherId() != -1) continue; // not primary track + if (mcTrack->GetMotherId() != -1) { + continue; // not primary track + } double z = mcTrack->GetStartZ(); if (!Is_MC_V || fabs(z - MC_Vcurr.Z()) > 1.e-7) { // new vertex Is_MC_V = 1; @@ -693,10 +790,13 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) mcTrack->GetStartVertex(MC_Vcurr); nvtrackscurr = 1; } - else + else { nvtrackscurr++; + } + } + if (nvtrackscurr > nvtracks) { + MC_V = MC_Vcurr; } - if (nvtrackscurr > nvtracks) MC_V = MC_Vcurr; } @@ -712,13 +812,19 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) TClonesArray TracksToFit("CbmStsTrack"); Int_t N = 0; for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { - if (Quality[iTrack] < 1) continue; + if (Quality[iTrack] < 1) { + continue; + } new (TracksToFit[N]) CbmStsTrack(); *(CbmStsTrack*) TracksToFit.At(N) = *(CbmStsTrack*) fRecStsTrackArray->At(iTrack); N++; - if (N % 50 != 0) continue; + if (N % 50 != 0) { + continue; + } Int_t i = N / 50; - if (i >= 13) continue; + if (i >= 13) { + continue; + } CbmVertex V; Finder.FindPrimaryVertex(&TracksToFit, &V); FillVertexHisto(MC_V, &V, fhV[i]); @@ -733,27 +839,45 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) if (fD0Analysis && fabs(MC_V.Z()) < 1.e-5) { // search for Kaon for (Int_t iK = 0; iK < nTracks; iK++) { - if (Quality[iK] < 1) continue; + if (Quality[iK] < 1) { + continue; + } CbmStsTrack* tK = (CbmStsTrack*) fRecStsTrackArray->At(iK); CbmTrackMatch* m = (CbmTrackMatch*) fSTSTrackMatch->At(iK); CbmMCTrack* mcK = (CbmMCTrack*) fMCTrackArray->At(m->GetMCTrackId()); - if (mcK->GetMotherId() != -1) continue; // not primary or D0 track - if (abs(mcK->GetPdgCode()) != 321) continue; + if (mcK->GetMotherId() != -1) { + continue; // not primary or D0 track + } + if (abs(mcK->GetPdgCode()) != 321) { + continue; + } double zK = mcK->GetPz(); - if (zK - MC_V.Z() < 1.e-5) continue; // primary + if (zK - MC_V.Z() < 1.e-5) { + continue; // primary + } double MCPK = mcK->GetP(); fFitter.DoFit(tK, 321); // refit // search for Pion for (Int_t iP = 0; iP < nTracks; iP++) { - if (Quality[iP] < 1) continue; + if (Quality[iP] < 1) { + continue; + } CbmStsTrack* tP = (CbmStsTrack*) fRecStsTrackArray->At(iP); m = (CbmTrackMatch*) fSTSTrackMatch->At(iP); CbmMCTrack* mcP = (CbmMCTrack*) fMCTrackArray->At(m->GetMCTrackId()); - if (mcP->GetMotherId() != -1) continue; // not primary or D0 track - if (abs(mcP->GetPdgCode()) != 211) continue; - if (mcK->GetPdgCode() * mcP->GetPdgCode() >= 0) continue; //same charge + if (mcP->GetMotherId() != -1) { + continue; // not primary or D0 track + } + if (abs(mcP->GetPdgCode()) != 211) { + continue; + } + if (mcK->GetPdgCode() * mcP->GetPdgCode() >= 0) { + continue; //same charge + } double zP = mcP->GetStartZ(); - if (fabs(zP - zK) > 1.e-8) continue; // different vertex + if (fabs(zP - zK) > 1.e-8) { + continue; // different vertex + } double MCPP = mcP->GetP(); const double D0 = 1.8645; @@ -794,7 +918,9 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) //SVFinder.GetFittedTrack(0, &KFitted ); //SVFinder.GetFittedTrack(1, &PFitted ); SVFinder.GetMass(&mass, &mass_err); - if (sv.GetNDF() <= 0) continue; + if (sv.GetNDF() <= 0) { + continue; + } Double_t dx = sv.GetX() - mc_.X(); Double_t dy = sv.GetY() - mc_.Y(); Double_t dz = sv.GetZ() - mc_.Z(); @@ -804,21 +930,31 @@ void CbmStsFitPerformanceTask::Exec(Option_t* /*option*/) fhD0[iconst][0]->Fill(1.E4 * dx); fhD0[iconst][1]->Fill(1.E4 * dy); fhD0[iconst][2]->Fill(1.E4 * dz); - if (sx > 1.e-10) fhD0[iconst][3]->Fill(dx / sqrt(sx)); - if (sy > 1.e-10) fhD0[iconst][4]->Fill(dy / sqrt(sy)); - if (sz > 1.e-10) fhD0[iconst][5]->Fill(dz / sqrt(sz)); + if (sx > 1.e-10) { + fhD0[iconst][3]->Fill(dx / sqrt(sx)); + } + if (sy > 1.e-10) { + fhD0[iconst][4]->Fill(dy / sqrt(sy)); + } + if (sz > 1.e-10) { + fhD0[iconst][5]->Fill(dz / sqrt(sz)); + } if (sv.GetNDF() > 0) { fhD0[iconst][6]->Fill(sv.GetChi2() / sv.GetNDF()); fhD0[iconst][7]->Fill(TMath::Prob(sv.GetChi2(), sv.GetNDF())); } fhD0[iconst][8]->Fill((mass - D0)); - if (mass_err > 1.e-10) fhD0[iconst][9]->Fill((mass - D0) / mass_err); + if (mass_err > 1.e-10) { + fhD0[iconst][9]->Fill((mass - D0) / mass_err); + } //fhD0[iconst][10]->Fill( ( fabs(1./KFitted.GetQp()) - MCPK)/MCPK ); //fhD0[iconst][11]->Fill( (fabs(1./PFitted.GetQp()) - MCPP)/MCPP ); - if (fabs(tK->GetParamFirst()->GetQp()) > 1.e-4 && MCPK > 1.e-4) + if (fabs(tK->GetParamFirst()->GetQp()) > 1.e-4 && MCPK > 1.e-4) { fhD0[iconst][12]->Fill((fabs(1. / tK->GetParamFirst()->GetQp()) - MCPK) / MCPK); - if (fabs(tP->GetParamFirst()->GetQp()) > 1.e-4 && MCPP > 1.e-4) + } + if (fabs(tP->GetParamFirst()->GetQp()) > 1.e-4 && MCPP > 1.e-4) { fhD0[iconst][13]->Fill((fabs(1. / tP->GetParamFirst()->GetQp()) - MCPP) / MCPP); + } } } } @@ -851,28 +987,44 @@ void CbmStsFitPerformanceTask::FillTrackHisto(const Double_t mc[6], CbmStsTrack* double t[6], c[15]; CbmKFMath::CopyTrackParam2TC(¶m, t, c); Bool_t ok = 1; - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { ok = ok && isfinite(t[i]); - for (int i = 0; i < 15; i++) + } + for (int i = 0; i < 15; i++) { ok = ok && isfinite(c[i]); - if (!ok) return; + } + if (!ok) { + return; + } //fhExtraTracks2ndMVD->Fill(t[0], t[1]); hist[0]->Fill(1.E4 * (t[0] - mc[0])); hist[1]->Fill(1.E4 * (t[1] - mc[1])); hist[2]->Fill(1.E3 * (t[2] - mc[2])); hist[3]->Fill(1.E3 * (t[3] - mc[3])); - if (fabs(t[4]) > 1.e-10) hist[4]->Fill((mc[4] / t[4] - 1.)); + if (fabs(t[4]) > 1.e-10) { + hist[4]->Fill((mc[4] / t[4] - 1.)); + } //if (z < 7.0) // first hit //if ( fabs( t[4] )>1.e-10 ) fhRes_vs_Mom_f->Fill(fabs(1.0/t[4]),1.E3*(t[2] - mc[2])); //if (z > 80.0) // last hit //if ( fabs( t[4] )>1.e-10 ) fhRes_vs_Mom_l->Fill(fabs(1.0/t[4]),1.E3*(t[2] - mc[2])); - if (c[0] > 1.e-10) hist[5]->Fill((t[0] - mc[0]) / sqrt(c[0])); - if (c[2] > 1.e-10) hist[6]->Fill((t[1] - mc[1]) / sqrt(c[2])); - if (c[5] > 1.e-10) hist[7]->Fill((t[2] - mc[2]) / sqrt(c[5])); - if (c[9] > 1.e-10) hist[8]->Fill((t[3] - mc[3]) / sqrt(c[9])); - if (c[14] > 1.e-10) hist[9]->Fill((t[4] - mc[4]) / sqrt(c[14])); + if (c[0] > 1.e-10) { + hist[5]->Fill((t[0] - mc[0]) / sqrt(c[0])); + } + if (c[2] > 1.e-10) { + hist[6]->Fill((t[1] - mc[1]) / sqrt(c[2])); + } + if (c[5] > 1.e-10) { + hist[7]->Fill((t[2] - mc[2]) / sqrt(c[5])); + } + if (c[9] > 1.e-10) { + hist[8]->Fill((t[3] - mc[3]) / sqrt(c[9])); + } + if (c[14] > 1.e-10) { + hist[9]->Fill((t[4] - mc[4]) / sqrt(c[14])); + } } // ------------------------------------------------------------------------- @@ -888,9 +1040,15 @@ void CbmStsFitPerformanceTask::FillVertexHisto(TVector3& mc, CbmVertex* V, TH1D* hist[0]->Fill(1.E4 * dx); hist[1]->Fill(1.E4 * dy); hist[2]->Fill(1.E4 * dz); - if (s2x > 1.e-10) hist[3]->Fill(dx / sqrt(s2x)); - if (s2y > 1.e-10) hist[4]->Fill(dy / sqrt(s2y)); - if (s2z > 1.e-10) hist[5]->Fill(dz / sqrt(s2z)); + if (s2x > 1.e-10) { + hist[3]->Fill(dx / sqrt(s2x)); + } + if (s2y > 1.e-10) { + hist[4]->Fill(dy / sqrt(s2y)); + } + if (s2z > 1.e-10) { + hist[5]->Fill(dz / sqrt(s2z)); + } hist[6]->Fill(V->GetChi2() / V->GetNDF()); hist[7]->Fill(TMath::Prob(V->GetChi2(), V->GetNDF())); hist[8]->Fill(V->GetNTracks()); @@ -899,10 +1057,14 @@ void CbmStsFitPerformanceTask::FillVertexHisto(TVector3& mc, CbmVertex* V, TH1D* // ------------------------------------------------------------------------- void CbmStsFitPerformanceTask::FillMCStateVectors(CbmStsPoint* point, Double_t mc[], Bool_t out) { - if (!point) return; + if (!point) { + return; + } Int_t mcTrackID = point->GetTrackID(); CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTrackArray->At(mcTrackID); - if (!mcTrack) return; + if (!mcTrack) { + return; + } Double_t q = GetCharge(mcTrack); // Get MC state vector @@ -916,7 +1078,9 @@ void CbmStsFitPerformanceTask::FillMCStateVectors(CbmStsPoint* point, Double_t m point->PositionOut(r); } double pzi = p.z(); - if (fabs(pzi) > 1.e-4) pzi = 1. / pzi; + if (fabs(pzi) > 1.e-4) { + pzi = 1. / pzi; + } mc[2] = p.x() * pzi; mc[3] = p.y() * pzi; mc[4] = p.Mag() > 1.e-4 ? q / p.Mag() : 0; @@ -927,10 +1091,14 @@ void CbmStsFitPerformanceTask::FillMCStateVectors(CbmStsPoint* point, Double_t m void CbmStsFitPerformanceTask::FillMCStateVectors(CbmMvdPoint* point, Double_t mc[], Bool_t out) { - if (!point) return; + if (!point) { + return; + } Int_t mcTrackID = point->GetTrackID(); CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTrackArray->At(mcTrackID); - if (!mcTrack) return; + if (!mcTrack) { + return; + } Double_t q = GetCharge(mcTrack); // Get MC state vector @@ -944,7 +1112,9 @@ void CbmStsFitPerformanceTask::FillMCStateVectors(CbmMvdPoint* point, Double_t m point->PositionOut(r); } double pzi = p.z(); - if (fabs(pzi) > 1.e-4) pzi = 1. / pzi; + if (fabs(pzi) > 1.e-4) { + pzi = 1. / pzi; + } mc[2] = p.x() * pzi; mc[3] = p.y() * pzi; mc[4] = p.Mag() > 1.e-4 ? q / p.Mag() : 0; @@ -960,10 +1130,12 @@ Double_t CbmStsFitPerformanceTask::GetCharge(CbmMCTrack* mcTrack) Double_t q; Int_t pdgCode = mcTrack->GetPdgCode(); TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(pdgCode); - if (particlePDG) + if (particlePDG) { q = particlePDG->Charge() / 3.; - else + } + else { q = 0; + } return q; } @@ -972,17 +1144,23 @@ Double_t CbmStsFitPerformanceTask::GetCharge(CbmMCTrack* mcTrack) Bool_t CbmStsFitPerformanceTask::IsLong(CbmStsTrack* track) { Int_t nHits = track->GetNofStsHits(); - if (nHits < 4) return 0; + if (nHits < 4) { + return 0; + } Int_t stmin = 1000, stmax = -1000; Int_t st, iHit, hitID; for (iHit = 0; iHit < nHits; iHit++) { hitID = track->GetStsHitIndex(iHit); st = ((CbmStsHit*) fStsHitArray->At(hitID))->GetAddress(); - if (st < stmin) stmin = st; - if (st > stmax) stmax = st; + if (st < stmin) { + stmin = st; + } + if (st > stmax) { + stmax = st; + } + } + if (stmax - stmin + 1 < 4) { + return 0; } - if (stmax - stmin + 1 < 4) return 0; return 1; } - -ClassImp(CbmStsFitPerformanceTask) diff --git a/reco/KF/Interface/CbmStsKFSecondaryVertexFinder.cxx b/reco/KF/Interface/CbmStsKFSecondaryVertexFinder.cxx index f599b9f6dc20806075e251ae19886075a3b1157e..e19e61a97c2a91cdb0511b341478224ad4ef8f61 100644 --- a/reco/KF/Interface/CbmStsKFSecondaryVertexFinder.cxx +++ b/reco/KF/Interface/CbmStsKFSecondaryVertexFinder.cxx @@ -28,7 +28,9 @@ ClassImp(CbmStsKFSecondaryVertexFinder) void CbmStsKFSecondaryVertexFinder::AddTrack(CbmStsTrack* Track) { - if (!Track) return; + if (!Track) { + return; + } vStsTracks.push_back(Track); vKFTracks.push_back(CbmKFTrack(*Track)); } @@ -85,7 +87,9 @@ void CbmStsKFSecondaryVertexFinder::GetFittedTrack( Int_t itrack, FairTrackParam */ void CbmStsKFSecondaryVertexFinder::GetMotherTrack(CbmStsTrack* MotherTrack) { - if (!MotherTrack) return; + if (!MotherTrack) { + return; + } double T[6], C[15]; Finder.GetMotherTrack(T, C); FairTrackParam parFirst(*MotherTrack->GetParamFirst()), parLast(*MotherTrack->GetParamLast()); diff --git a/reco/KF/Interface/CbmStsKFTrackFitter.cxx b/reco/KF/Interface/CbmStsKFTrackFitter.cxx index 4373310d6aabc24458b5efe17321e8d4ec3f00a4..8d545e5ce7495cea88b4b4eebcef6adab7f6bb35 100644 --- a/reco/KF/Interface/CbmStsKFTrackFitter.cxx +++ b/reco/KF/Interface/CbmStsKFTrackFitter.cxx @@ -58,7 +58,9 @@ void CbmStsKFTrackFitter::SetKFHits(CbmKFTrack& T, CbmStsTrack* track) T.fHits.clear(); - if (!fIsInitialised) Init(); + if (!fIsInitialised) { + Init(); + } Int_t NStsHits = (fStsHitsArray) ? track->GetNofStsHits() : 0; Int_t NMvdHits = (fMvdHitsArray) ? track->GetNofMvdHits() : 0; @@ -87,9 +89,10 @@ Int_t CbmStsKFTrackFitter::DoFit(CbmStsTrack* track, Int_t pidHypo) CbmKFTrack T; T.SetPID(pidHypo); SetKFHits(T, track); - for (Int_t i = 0; i < 6; i++) + for (Int_t i = 0; i < 6; i++) { T.GetTrack()[i] = 0.; // no guess taken - T.Fit(1); // fit downstream + } + T.Fit(1); // fit downstream CheckTrack(T); T.Fit(0); // fit upstream CheckTrack(T); @@ -102,15 +105,19 @@ Int_t CbmStsKFTrackFitter::DoFit(CbmStsTrack* track, Int_t pidHypo) track->SetParamLast(&par); err = T.Fit(0); // fit upstream ok = ok && (!err) && CheckTrack(T); - if (ok) T.GetStsTrack(*track, 1); // store fitted track & cov.matrix & chi2 & NDF + if (ok) { + T.GetStsTrack(*track, 1); // store fitted track & cov.matrix & chi2 & NDF + } } if (!ok) { Double_t* t = T.GetTrack(); Double_t* c = T.GetCovMatrix(); - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { t[i] = 0; - for (int i = 0; i < 15; i++) + } + for (int i = 0; i < 15; i++) { c[i] = 0; + } c[0] = c[2] = c[5] = c[9] = c[14] = 100.; T.GetRefChi2() = 100.; T.GetRefNDF() = 0; @@ -127,17 +134,23 @@ Int_t CbmStsKFTrackFitter::DoFit(CbmStsTrack* track, Int_t pidHypo) void CbmStsKFTrackFitter::Extrapolate(FairTrackParam* track, Double_t z, FairTrackParam* e_track) { - if (!track) return; + if (!track) { + return; + } CbmKFTrack T; T.SetTrackParam(*track); T.Extrapolate(z); - if (e_track) T.GetTrackParam(*e_track); + if (e_track) { + T.GetTrackParam(*e_track); + } } void CbmStsKFTrackFitter::Extrapolate(CbmStsTrack* track, Double_t z, FairTrackParam* e_track) { - if (!track) return; + if (!track) { + return; + } CbmKFTrack T; T.SetPID(track->GetPidHypo()); const FairTrackParam *fpar = track->GetParamFirst(), *lpar = track->GetParamLast(); @@ -163,7 +176,9 @@ void CbmStsKFTrackFitter::Extrapolate(CbmStsTrack* track, Double_t z, FairTrackP T.SetTrackParam(*fpar); T.Smooth(z); } - if (e_track) T.GetTrackParam(*e_track); + if (e_track) { + T.GetTrackParam(*e_track); + } } @@ -199,14 +214,18 @@ Double_t CbmStsKFTrackFitter::GetChiToVertex(CbmStsTrack* track, CbmVertex* vtx) Double_t CbmStsKFTrackFitter::FitToVertex(CbmStsTrack* track, CbmVertex* vtx, FairTrackParam* v_track) { Double_t ret = 100.; - if (!track || !vtx || !v_track) return ret; + if (!track || !vtx || !v_track) { + return ret; + } CbmKFTrack T(*track); CbmKFVertex V(*vtx); T.Fit2Vertex(V); T.GetTrackParam(*v_track); if (T.GetRefNDF() > 0 && T.GetRefChi2() >= 0) { ret = T.GetRefChi2() / T.GetRefNDF(); - if (isfinite(ret)) ret = sqrt(ret); + if (isfinite(ret)) { + ret = sqrt(ret); + } } return ret; } @@ -216,10 +235,12 @@ Bool_t CbmStsKFTrackFitter::CheckTrack(CbmKFTrack& T) Bool_t ok = 1; Double_t* t = T.GetTrack(); Double_t* c = T.GetCovMatrix(); - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { ok = ok && isfinite(t[i]) && TMath::Abs(t[i]) < 1.e5; - for (int i = 0; i < 15; i++) + } + for (int i = 0; i < 15; i++) { ok = ok && isfinite(c[i]); + } ok = ok && isfinite(T.GetMass()) && isfinite(T.GetRefChi2()); if (ok) { ok = ok && (c[0] > 0); @@ -235,41 +256,83 @@ Bool_t CbmStsKFTrackFitter::CheckTrack(CbmKFTrack& T) Double_t c33 = TMath::Sqrt(c[9]); Double_t c44 = TMath::Sqrt(c[14]); Double_t a = c11 * c00; - if (c[1] > a) c[1] = a; - if (c[1] < -a) c[1] = -a; + if (c[1] > a) { + c[1] = a; + } + if (c[1] < -a) { + c[1] = -a; + } a = c22 * c00; - if (c[3] > a) c[3] = a; - if (c[3] < -a) c[3] = -a; + if (c[3] > a) { + c[3] = a; + } + if (c[3] < -a) { + c[3] = -a; + } a = c22 * c11; - if (c[4] > a) c[4] = a; - if (c[4] < -a) c[4] = -a; + if (c[4] > a) { + c[4] = a; + } + if (c[4] < -a) { + c[4] = -a; + } a = c33 * c00; - if (c[6] > a) c[6] = a; - if (c[6] < -a) c[6] = -a; + if (c[6] > a) { + c[6] = a; + } + if (c[6] < -a) { + c[6] = -a; + } a = c33 * c11; - if (c[7] > a) c[7] = a; - if (c[7] < -a) c[7] = -a; + if (c[7] > a) { + c[7] = a; + } + if (c[7] < -a) { + c[7] = -a; + } a = c33 * c22; - if (c[8] > a) c[8] = a; - if (c[8] < -a) c[8] = -a; + if (c[8] > a) { + c[8] = a; + } + if (c[8] < -a) { + c[8] = -a; + } a = c44 * c00; - if (c[10] > a) c[10] = a; - if (c[10] < -a) c[10] = -a; + if (c[10] > a) { + c[10] = a; + } + if (c[10] < -a) { + c[10] = -a; + } a = c44 * c11; - if (c[11] > a) c[11] = a; - if (c[11] < -a) c[11] = -a; + if (c[11] > a) { + c[11] = a; + } + if (c[11] < -a) { + c[11] = -a; + } a = c44 * c22; - if (c[12] > a) c[12] = a; - if (c[12] < -a) c[12] = -a; + if (c[12] > a) { + c[12] = a; + } + if (c[12] < -a) { + c[12] = -a; + } a = c44 * c33; - if (c[13] > a) c[13] = a; - if (c[13] < -a) c[13] = -a; + if (c[13] > a) { + c[13] = a; + } + if (c[13] < -a) { + c[13] = -a; + } } if (!ok) { - for (int i = 0; i < 6; i++) + for (int i = 0; i < 6; i++) { t[i] = 0; - for (int i = 0; i < 15; i++) + } + for (int i = 0; i < 15; i++) { c[i] = 0; + } c[0] = c[2] = c[5] = c[9] = c[14] = 100.; T.GetRefChi2() = 100.; T.GetRefNDF() = 0; diff --git a/reco/KF/Interface/CbmTofTrackFitterKF.cxx b/reco/KF/Interface/CbmTofTrackFitterKF.cxx index 06bcadee40328a6beae29e07ea12c93ae88fca7f..38a42e69b6129fc7b466d6b186519aea9319e389 100644 --- a/reco/KF/Interface/CbmTofTrackFitterKF.cxx +++ b/reco/KF/Interface/CbmTofTrackFitterKF.cxx @@ -93,7 +93,9 @@ Int_t CbmTofTrackFitterKF::DoFit(CbmTofTracklet* pTrack) { LOG(debug1) << "CbmTofTrackFitterKF::DoFit starting "; // Implementation of the fitting algorithm - if (nullptr == fArrayTofHit) this->Init(); + if (nullptr == fArrayTofHit) { + this->Init(); + } if (nullptr == fArrayTofHit) { LOG(error) << "CbmTofTrackFitterKF::DoFit No fArrayTofHit "; @@ -177,11 +179,15 @@ Int_t CbmTofTrackFitterKF::DoFit(CbmTofTrack* /*pTrack*/) { return 0; } void CbmTofTrackFitterKF::Extrapolate(const FairTrackParam* track, Double_t z, FairTrackParam* e_track) { - if (!track) return; + if (!track) { + return; + } CbmKFTrack T; T.SetTrackParam(*track); T.Extrapolate(z); - if (e_track) T.GetTrackParam(*e_track); + if (e_track) { + T.GetTrackParam(*e_track); + } } /* diff --git a/reco/KF/Interface/CbmTrdTrackFitterKF.cxx b/reco/KF/Interface/CbmTrdTrackFitterKF.cxx index 082e84326d0f7313b077fd7a705060be5d0f91c8..8295aea50ecc27a33a84deed27f526f0914ecd4d 100644 --- a/reco/KF/Interface/CbmTrdTrackFitterKF.cxx +++ b/reco/KF/Interface/CbmTrdTrackFitterKF.cxx @@ -91,7 +91,9 @@ void CbmTrdTrackFitterKF::Init() Int_t CbmTrdTrackFitterKF::DoFit(CbmTrdTrack* pTrack) { // Implementation of the fitting algorithm - if (nullptr == fArrayTrdHit) return 1; + if (nullptr == fArrayTrdHit) { + return 1; + } // Declare variables outside the loop CbmTrdHit* pHit = nullptr; diff --git a/reco/KF/KFQA/CbmKFTrErrMCPoints.cxx b/reco/KF/KFQA/CbmKFTrErrMCPoints.cxx index 3f30e4fa882f79efc618e26ea97baa143054139c..3b32fa6d352eabd7a0ab9bc422364086ebc5a603 100644 --- a/reco/KF/KFQA/CbmKFTrErrMCPoints.cxx +++ b/reco/KF/KFQA/CbmKFTrErrMCPoints.cxx @@ -20,20 +20,15 @@ using std::vector; -ClassImp(CbmKFTrErrMCPoints) - - CbmKFTrErrMCPoints::CbmKFTrErrMCPoints() - : StsArray() - , MvdArray() - , TofArray() - , StsHitsArray() - , MvdHitsArray() -{ -} +ClassImp(CbmKFTrErrMCPoints); + +CbmKFTrErrMCPoints::CbmKFTrErrMCPoints() : StsArray(), MvdArray(), TofArray(), StsHitsArray(), MvdHitsArray() {} int CbmKFTrErrMCPoints::GetNConsMCStations() { - if ((GetNMvdPoints() + GetNStsPoints()) < 1) return 0; + if ((GetNMvdPoints() + GetNStsPoints()) < 1) { + return 0; + } // TODO get station number of the point using methods of the point class! float zStation[8] = {30., 40., 50., 60., 70., 80., 90., 100.}; vector<int> iStations; @@ -43,16 +38,23 @@ int CbmKFTrErrMCPoints::GetNConsMCStations() } for (int iSts = 0; iSts < GetNStsPoints(); ++iSts) { int stNumber = -1; - for (int iSt = 0; iSt < 8; iSt++) - if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5) stNumber = iSt; - if (stNumber >= 0) iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations()); + for (int iSt = 0; iSt < 8; iSt++) { + if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5) { + stNumber = iSt; + } + } + if (stNumber >= 0) { + iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations()); + } } std::sort(iStations.begin(), iStations.end()); int nMaxConsStations = 1; int nConsStations = 1; - if (iStations.size() == 0) return 0; + if (iStations.size() == 0) { + return 0; + } int iPrevSt = iStations[0]; for (unsigned int iP = 1; iP < iStations.size(); iP++) { if ((iStations[iP] - iPrevSt) == 1) { @@ -65,23 +67,28 @@ int CbmKFTrErrMCPoints::GetNConsMCStations() iPrevSt = iStations[iP]; } } - if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations; + if (nConsStations > nMaxConsStations) { + nMaxConsStations = nConsStations; + } return nMaxConsStations; } int CbmKFTrErrMCPoints::GetNConsHitStations() { - if ((GetNMvdHits() + GetNStsHits()) < 1) return 0; + if ((GetNMvdHits() + GetNStsHits()) < 1) { + return 0; + } // TODO get station number of the point using methods of the point class! vector<int> iStations; for (int iMvd = 0; iMvd < GetNMvdHits(); ++iMvd) { iStations.push_back(GetMvdHit(iMvd)->GetStationNr() - 1); // std::cout << GetMvdHit(iMvd)->GetStationNr() << " " << GetMvdHit(iMvd)->GetZ() << std::endl; } - for (int iSts = 0; iSts < GetNStsHits(); ++iSts) + for (int iSts = 0; iSts < GetNStsHits(); ++iSts) { iStations.push_back(CbmStsSetup::Instance()->GetStationNumber(GetStsHit(iSts)->GetAddress()) - 1 + CbmKF::Instance()->GetNMvdStations()); + } std::sort(iStations.begin(), iStations.end()); @@ -94,29 +101,35 @@ int CbmKFTrErrMCPoints::GetNConsHitStations() iPrevSt = iStations[iP]; } else if ((iStations[iP] - iPrevSt) > 1) { - if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations; + if (nConsStations > nMaxConsStations) { + nMaxConsStations = nConsStations; + } nConsStations = 1; iPrevSt = iStations[iP]; } } - if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations; + if (nConsStations > nMaxConsStations) { + nMaxConsStations = nConsStations; + } return nMaxConsStations; } int CbmKFTrErrMCPoints::GetNHitStations() { - if ((GetNMvdHits() + GetNStsHits()) < 1) return 0; + if ((GetNMvdHits() + GetNStsHits()) < 1) { + return 0; + } // TODO get station number of the point using methods of the point class! vector<int> iStations; for (int iMvd = 0; iMvd < GetNMvdHits(); ++iMvd) { iStations.push_back(GetMvdHit(iMvd)->GetStationNr() - 1); // std::cout << GetMvdHit(iMvd)->GetStationNr() << " " << GetMvdHit(iMvd)->GetZ() << std::endl; } - for (int iSts = 0; iSts < GetNStsHits(); ++iSts) + for (int iSts = 0; iSts < GetNStsHits(); ++iSts) { iStations.push_back(CbmStsSetup::Instance()->GetStationNumber(GetStsHit(iSts)->GetAddress()) - 1 + CbmKF::Instance()->GetNMvdStations()); - + } std::sort(iStations.begin(), iStations.end()); int nStations = 1; @@ -133,7 +146,9 @@ int CbmKFTrErrMCPoints::GetNHitStations() int CbmKFTrErrMCPoints::GetNMaxMCPointsOnStation() { - if ((GetNMvdPoints() + GetNStsPoints()) < 1) return 0; + if ((GetNMvdPoints() + GetNStsPoints()) < 1) { + return 0; + } // TODO get station number of the point using methods of the point class! float zStation[8] = {30., 40., 50., 60., 70., 80., 90., 100.}; vector<int> iStations; @@ -143,16 +158,23 @@ int CbmKFTrErrMCPoints::GetNMaxMCPointsOnStation() } for (int iSts = 0; iSts < GetNStsPoints(); ++iSts) { int stNumber = -1; - for (int iSt = 0; iSt < 8; iSt++) - if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5) stNumber = iSt; - if (stNumber >= 0) iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations()); + for (int iSt = 0; iSt < 8; iSt++) { + if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5) { + stNumber = iSt; + } + } + if (stNumber >= 0) { + iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations()); + } } std::sort(iStations.begin(), iStations.end()); int nMaxMCPointsOnStation = 1; int nMCPointsOnStation = 1; - if (iStations.size() == 0) return 0; + if (iStations.size() == 0) { + return 0; + } int iPrevSt = iStations[0]; for (unsigned int iP = 1; iP < iStations.size(); iP++) { if ((iStations[iP] - iPrevSt) == 0) { @@ -160,7 +182,9 @@ int CbmKFTrErrMCPoints::GetNMaxMCPointsOnStation() iPrevSt = iStations[iP]; } else { - if (nMCPointsOnStation > nMaxMCPointsOnStation) nMaxMCPointsOnStation = nMCPointsOnStation; + if (nMCPointsOnStation > nMaxMCPointsOnStation) { + nMaxMCPointsOnStation = nMCPointsOnStation; + } nMCPointsOnStation = 1; iPrevSt = iStations[iP]; } @@ -176,9 +200,15 @@ Bool_t CbmKFTrErrMCPoints::IsReconstructable(CbmMCTrack* mcTr, int MinNStations, // reject very slow tracks from analysis f &= (mcTr->GetP() > MinRecoMom); // detected at least in 4 stations - if (PerformanceMode == 3) f &= (GetNConsMCStations() >= MinNStations); // L1-MC - if (PerformanceMode == 2) f &= (GetNHitStations() >= MinNStations); // QA definition - if (PerformanceMode == 1) f &= (GetNConsHitStations() >= MinNStations); // L1 definition + if (PerformanceMode == 3) { + f &= (GetNConsMCStations() >= MinNStations); // L1-MC + } + if (PerformanceMode == 2) { + f &= (GetNHitStations() >= MinNStations); // QA definition + } + if (PerformanceMode == 1) { + f &= (GetNConsHitStations() >= MinNStations); // L1 definition + } // maximul 4 layers for a station. f &= (GetNMaxMCPointsOnStation() <= 4); diff --git a/reco/KF/KFQA/CbmKFTrackFitQa.cxx b/reco/KF/KFQA/CbmKFTrackFitQa.cxx index 6b12b2b6b8c640fb60353065cbe855b70a1f48d3..b0088acb5c5b0c0f40aa2c90f30ec50e11ae4992 100644 --- a/reco/KF/KFQA/CbmKFTrackFitQa.cxx +++ b/reco/KF/KFQA/CbmKFTrackFitQa.cxx @@ -172,39 +172,39 @@ ClassImp(CbmKFTrackFitQa) CbmKFTrackFitQa::~CbmKFTrackFitQa() { - if (res_STShit_x) delete res_STShit_x; - if (res_STShit_y) delete res_STShit_y; - if (pull_STShit_x) delete pull_STShit_x; - if (pull_STShit_y) delete pull_STShit_y; - - if (res_MVDhit_x) delete res_MVDhit_x; - if (res_MVDhit_y) delete res_MVDhit_y; - if (pull_MVDhit_x) delete pull_MVDhit_x; - if (pull_MVDhit_y) delete pull_MVDhit_y; - - if (res_AtPV_x) delete res_AtPV_x; - if (res_AtPV_y) delete res_AtPV_y; - if (res_AtPV_tx) delete res_AtPV_tx; - if (res_AtPV_ty) delete res_AtPV_ty; - if (res_AtPV_qp) delete res_AtPV_qp; - - if (pull_AtPV_x) delete pull_AtPV_x; - if (pull_AtPV_y) delete pull_AtPV_y; - if (pull_AtPV_tx) delete pull_AtPV_tx; - if (pull_AtPV_ty) delete pull_AtPV_ty; - if (pull_AtPV_qp) delete pull_AtPV_qp; - - if (res_AtFP_x) delete res_AtFP_x; - if (res_AtFP_y) delete res_AtFP_y; - if (res_AtFP_tx) delete res_AtFP_tx; - if (res_AtFP_ty) delete res_AtFP_ty; - if (res_AtFP_qp) delete res_AtFP_qp; - - if (pull_AtFP_x) delete pull_AtFP_x; - if (pull_AtFP_y) delete pull_AtFP_y; - if (pull_AtFP_tx) delete pull_AtFP_tx; - if (pull_AtFP_ty) delete pull_AtFP_ty; - if (pull_AtFP_qp) delete pull_AtFP_qp; + delete res_STShit_x; + delete res_STShit_y; + delete pull_STShit_x; + delete pull_STShit_y; + + delete res_MVDhit_x; + delete res_MVDhit_y; + delete pull_MVDhit_x; + delete pull_MVDhit_y; + + delete res_AtPV_x; + delete res_AtPV_y; + delete res_AtPV_tx; + delete res_AtPV_ty; + delete res_AtPV_qp; + + delete pull_AtPV_x; + delete pull_AtPV_y; + delete pull_AtPV_tx; + delete pull_AtPV_ty; + delete pull_AtPV_qp; + + delete res_AtFP_x; + delete res_AtFP_y; + delete res_AtFP_tx; + delete res_AtFP_ty; + delete res_AtFP_qp; + + delete pull_AtFP_x; + delete pull_AtFP_y; + delete pull_AtFP_tx; + delete pull_AtFP_ty; + delete pull_AtFP_qp; } InitStatus CbmKFTrackFitQa::ReInit() { return Init(); } @@ -241,12 +241,13 @@ void CbmKFTrackFitQa::Exec(Option_t* /*option*/) FillHitHistos(); CbmKFTrErrMCPoints* MCTrackSortedArray = new CbmKFTrErrMCPoints[listMCTracks->GetEntriesFast() + 1]; - if (CbmKF::Instance()->vMvdMaterial.size() != 0) + if (CbmKF::Instance()->vMvdMaterial.size() != 0) { for (Int_t iMvd = 0; iMvd < listMvdPts->GetEntriesFast(); iMvd++) { CbmMvdPoint* MvdPoint = (CbmMvdPoint*) listMvdPts->At(iMvd); //MCTrackSortedArray[MvdPoint->GetTrackID()].SetMvdPoint(MvdPoint); MCTrackSortedArray[MvdPoint->GetTrackID()].MvdArray.push_back(MvdPoint); } + } for (Int_t iSts = 0; iSts < listStsPts->GetEntriesFast(); iSts++) { CbmStsPoint* StsPoint = (CbmStsPoint*) listStsPts->At(iSts); //MCTrackSortedArray[StsPoint->GetTrackID()].SetStsPoint(StsPoint); @@ -255,7 +256,9 @@ void CbmKFTrackFitQa::Exec(Option_t* /*option*/) for (Int_t itrack = 0; itrack < listStsTracks->GetEntriesFast(); itrack++) { CbmStsTrack* StsTrack = (CbmStsTrack*) listStsTracks->At(itrack); CbmTrackMatch* StsTrackMatch = (CbmTrackMatch*) listStsTracksMatch->At(itrack); - if (StsTrackMatch->GetNofMCTracks() == 0) continue; + if (StsTrackMatch->GetNofMCTracks() == 0) { + continue; + } CbmMCTrack* MCTrack = (CbmMCTrack*) listMCTracks->At(StsTrackMatch->GetMCTrackId()); CbmKFTrack KFTrack(*StsTrack); FillHistoAtFirstPoint(&MCTrackSortedArray[StsTrackMatch->GetMCTrackId()], MCTrack, &KFTrack); @@ -282,10 +285,12 @@ void CbmKFTrackFitQa::FillHistoAtParticleVertex(CbmMCTrack* track_mc, CbmKFTrack // getting the q (charge) of the mc track using the pdg information Double_t qtrack = 0; - if (track_mc->GetPdgCode() < 9999999) + if (track_mc->GetPdgCode() < 9999999) { qtrack = TDatabasePDG::Instance()->GetParticle(track_mc->GetPdgCode())->Charge() / 3.0; - else + } + else { qtrack = 0; + } Double_t PointPx = track_mc->GetPx(); Double_t PointPy = track_mc->GetPy(); @@ -307,21 +312,37 @@ void CbmKFTrackFitQa::FillHistoAtParticleVertex(CbmMCTrack* track_mc, CbmKFTrack res_AtPV_y->Fill(ddy); res_AtPV_tx->Fill(ddtx); res_AtPV_ty->Fill(ddty); - if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) res_AtPV_qp->Fill(ddqp); + if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) { + res_AtPV_qp->Fill(ddqp); + } //pulls of the parameters - if (isfinite(fC[0]) && fC[0] > 0) pull_AtPV_x->Fill(ddx / sqrt(fC[0])); - if (isfinite(fC[2]) && fC[2] > 0) pull_AtPV_y->Fill(ddy / sqrt(fC[2])); - if (isfinite(fC[5]) && fC[5] > 0) pull_AtPV_tx->Fill(ddtx / sqrt(fC[5])); - if (isfinite(fC[9]) && fC[9] > 0) pull_AtPV_ty->Fill(ddty / sqrt(fC[9])); - if (isfinite(fC[14]) && fC[14] > 0) pull_AtPV_qp->Fill(ddqp_p / sqrt(fC[14])); + if (isfinite(fC[0]) && fC[0] > 0) { + pull_AtPV_x->Fill(ddx / sqrt(fC[0])); + } + if (isfinite(fC[2]) && fC[2] > 0) { + pull_AtPV_y->Fill(ddy / sqrt(fC[2])); + } + if (isfinite(fC[5]) && fC[5] > 0) { + pull_AtPV_tx->Fill(ddtx / sqrt(fC[5])); + } + if (isfinite(fC[9]) && fC[9] > 0) { + pull_AtPV_ty->Fill(ddty / sqrt(fC[9])); + } + if (isfinite(fC[14]) && fC[14] > 0) { + pull_AtPV_qp->Fill(ddqp_p / sqrt(fC[14])); + } if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) { - if (qtrack == (fabs(fT[4]) / fT[4])) + if (qtrack == (fabs(fT[4]) / fT[4])) { q_QA->Fill(P_mc, 100.0); - else + } + else { q_QA->Fill(P_mc, 0.0); + } - if (isfinite(fC[14]) && fC[14] > 0) dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1); + if (isfinite(fC[14]) && fC[14] > 0) { + dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1); + } } } @@ -339,13 +360,19 @@ void CbmKFTrackFitQa::FillHistoAtFirstPoint(CbmKFTrErrMCPoints* mc_points, CbmMC // if (fabs(track_mc->GetPz()) < 0.0001) return; //if (fT[4]<0.) return; - if (!mc_points) return; - if ((mc_points->MvdArray.size() == 0) && (mc_points->StsArray.size() == 0)) return; + if (!mc_points) { + return; + } + if ((mc_points->MvdArray.size() == 0) && (mc_points->StsArray.size() == 0)) { + return; + } - if (mc_points->MvdArray.size() > 0) + if (mc_points->MvdArray.size() > 0) { MCFirstPoint = *mc_points->MvdArray.begin(); - else + } + else { MCFirstPoint = *mc_points->StsArray.begin(); + } for (vector<CbmMvdPoint*>::iterator imvd = mc_points->MvdArray.begin(); imvd != mc_points->MvdArray.end(); ++imvd) { MCFirstPoint = *imvd; @@ -372,10 +399,12 @@ void CbmKFTrackFitQa::FillHistoAtFirstPoint(CbmKFTrErrMCPoints* mc_points, CbmMC // getting of the q (charge) of the mc track using the pdg information { Double_t qtrack = 0; - if (track_mc->GetPdgCode() < 9999999) + if (track_mc->GetPdgCode() < 9999999) { qtrack = TDatabasePDG::Instance()->GetParticle(track_mc->GetPdgCode())->Charge() / 3.0; - else + } + else { qtrack = 0; + } Double_t PointPx = MCFirstPoint->GetPx(); Double_t PointPy = MCFirstPoint->GetPy(); @@ -383,6 +412,7 @@ void CbmKFTrackFitQa::FillHistoAtFirstPoint(CbmKFTrErrMCPoints* mc_points, CbmMC Double_t P_mc = sqrt(PointPx * PointPx + PointPy * PointPy + PointPz * PointPz); //differences of the KF track and MC track parameters calculation + Double_t ddx = fT[0] - MCFirstPoint->GetX(); // The difference between x coordinates of the reconstructed and MC tracks Double_t ddy = @@ -397,23 +427,39 @@ void CbmKFTrackFitQa::FillHistoAtFirstPoint(CbmKFTrErrMCPoints* mc_points, CbmMC res_AtFP_y->Fill(ddy * 10000.); res_AtFP_tx->Fill(ddtx); res_AtFP_ty->Fill(ddty); - if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) res_AtFP_qp->Fill(ddqp); + if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) { + res_AtFP_qp->Fill(ddqp); + } //pulls of the parameters - if (isfinite(fC[0]) && fC[0] > 0) pull_AtFP_x->Fill(ddx / sqrt(fC[0])); - if (isfinite(fC[2]) && fC[2] > 0) pull_AtFP_y->Fill(ddy / sqrt(fC[2])); - if (isfinite(fC[5]) && fC[5] > 0) pull_AtFP_tx->Fill(ddtx / sqrt(fC[5])); - if (isfinite(fC[9]) && fC[9] > 0) pull_AtFP_ty->Fill(ddty / sqrt(fC[9])); - if (isfinite(fC[14]) && fC[14] > 0) pull_AtFP_qp->Fill(ddqp_p / sqrt(fC[14])); + if (isfinite(fC[0]) && fC[0] > 0) { + pull_AtFP_x->Fill(ddx / sqrt(fC[0])); + } + if (isfinite(fC[2]) && fC[2] > 0) { + pull_AtFP_y->Fill(ddy / sqrt(fC[2])); + } + if (isfinite(fC[5]) && fC[5] > 0) { + pull_AtFP_tx->Fill(ddtx / sqrt(fC[5])); + } + if (isfinite(fC[9]) && fC[9] > 0) { + pull_AtFP_ty->Fill(ddty / sqrt(fC[9])); + } + if (isfinite(fC[14]) && fC[14] > 0) { + pull_AtFP_qp->Fill(ddqp_p / sqrt(fC[14])); + } if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) { - if (qtrack == (fabs(fT[4]) / fT[4])) + if (qtrack == (fabs(fT[4]) / fT[4])) { q_QA->Fill(P_mc, 100.0); - else + } + else { q_QA->Fill(P_mc, 0.0); + } - if (isfinite(fC[14]) && fC[14] > 0) dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1); + if (isfinite(fC[14]) && fC[14] > 0) { + dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1); + } + //cout << ddx <<" "<< ddx/sqrt(fC[0]) << endl; } - //cout << ddx <<" "<< ddx/sqrt(fC[0]) << endl; } } @@ -601,8 +647,9 @@ void CbmKFTrackFitQa::StsHitMatch() CbmLink link3 = stsDigiMatch->GetLink(iL3); int stsPointId = link3.GetIndex(); - if (!find(&(stsPointIds[0]), &(stsPointIds[stsPointIds.size()]), stsPointId)) + if (!find(&(stsPointIds[0]), &(stsPointIds[stsPointIds.size()]), stsPointId)) { continue; // check if first cluster matched with same mc-point + } // CbmStsPoint *pt = dynamic_cast<CbmStsPoint*>( listStsPts->At(stsPointId) ); // if(!pt) continue; // CbmMCTrack *MCTrack = dynamic_cast<CbmMCTrack*>( listMCTracks->At( pt->GetTrackID() ) ); diff --git a/reco/KF/KFQA/CbmKFTrackQa.cxx b/reco/KF/KFQA/CbmKFTrackQa.cxx index 94445636d20de8b01dbb60007602e40ee1c27ba4..7120c05a2be4b8fa126e30af4887fdb5e4e84d93 100644 --- a/reco/KF/KFQA/CbmKFTrackQa.cxx +++ b/reco/KF/KFQA/CbmKFTrackQa.cxx @@ -92,10 +92,12 @@ CbmKFTrackQa::CbmKFTrackQa(const char* name, Int_t iVerbose, TString outFileName TFile* curFile = gFile; TDirectory* curDirectory = gDirectory; - if (!(fOutFileName == "")) + if (!(fOutFileName == "")) { fOutFile = new TFile(fOutFileName.Data(), "RECREATE"); - else + } + else { fOutFile = gFile; + } fOutFile->cd(); @@ -297,21 +299,31 @@ InitStatus CbmKFTrackQa::Init() // Get global tracks fGlobalTrackArray = (TClonesArray*) ioman->GetObject(fGlobalTrackBranchName); - if (fGlobalTrackArray == nullptr) Warning("CbmKFTrackQa::Init", "global track array not found!"); + if (fGlobalTrackArray == nullptr) { + Warning("CbmKFTrackQa::Init", "global track array not found!"); + } // Get ToF hits fTofHitArray = (TClonesArray*) ioman->GetObject(fTofBranchName); - if (fTofHitArray == nullptr) Warning("CbmKFTrackQa::Init", "TOF hit-array not found!"); + if (fTofHitArray == nullptr) { + Warning("CbmKFTrackQa::Init", "TOF hit-array not found!"); + } // TRD fTrdTrackArray = (TClonesArray*) ioman->GetObject(fTrdBranchName); - if (fTrdTrackArray == nullptr) Warning("CbmKFTrackQa::Init", "TRD track-array not found!"); + if (fTrdTrackArray == nullptr) { + Warning("CbmKFTrackQa::Init", "TRD track-array not found!"); + } fTrdHitArray = (TClonesArray*) ioman->GetObject(fTrdHitBranchName); - if (fTrdHitArray == nullptr) Warning("CbmKFTrackQa::Init", "TRD hit-array not found!"); + if (fTrdHitArray == nullptr) { + Warning("CbmKFTrackQa::Init", "TRD hit-array not found!"); + } fRichRingArray = (TClonesArray*) ioman->GetObject(fRichBranchName); - if (fRichRingArray == nullptr) Warning("CbmKFTrackQa::Init", "Rich ring array not found!"); + if (fRichRingArray == nullptr) { + Warning("CbmKFTrackQa::Init", "Rich ring array not found!"); + } fMCTrackArray = (TClonesArray*) ioman->GetObject(fMCTracksBranchName); if (fMCTrackArray == nullptr) { @@ -328,15 +340,21 @@ InitStatus CbmKFTrackQa::Init() //Ring match fRichRingMatchArray = (TClonesArray*) ioman->GetObject(fRichRingMatchBranchName); - if (fRichRingMatchArray == nullptr) Warning("CbmKFTrackQa::Init", "RichRing match array not found!"); + if (fRichRingMatchArray == nullptr) { + Warning("CbmKFTrackQa::Init", "RichRing match array not found!"); + } //Tof match fTofHitMatchArray = (TClonesArray*) ioman->GetObject(fTofHitMatchBranchName); - if (fTofHitMatchArray == nullptr) Warning("CbmKFTrackQa::Init", "TofHit match array not found!"); + if (fTofHitMatchArray == nullptr) { + Warning("CbmKFTrackQa::Init", "TofHit match array not found!"); + } //TRD match fTrdTrackMatchArray = (TClonesArray*) ioman->GetObject(fTrdTrackMatchBranchName); - if (fTrdTrackMatchArray == nullptr) Warning("CbmKFTrackQa::Init", "TrdTrack match array not found!"); + if (fTrdTrackMatchArray == nullptr) { + Warning("CbmKFTrackQa::Init", "TrdTrack match array not found!"); + } //Much track match fMuchTrackMatchArray = (TClonesArray*) ioman->GetObject(fMuchTrackMatchBranchName); @@ -383,26 +401,36 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) Int_t pdg = cbmMCTrack->GetPdgCode(); Double_t q = 1; - if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) + if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; - else if (pdg == 1000010020) + } + else if (pdg == 1000010020) { q = 1; - else if (pdg == -1000010020) + } + else if (pdg == -1000010020) { q = -1; - else if (pdg == 1000010030) + } + else if (pdg == 1000010030) { q = 1; - else if (pdg == -1000010030) + } + else if (pdg == -1000010030) { q = -1; - else if (pdg == 1000020030) + } + else if (pdg == 1000020030) { q = 2; - else if (pdg == -1000020030) + } + else if (pdg == -1000020030) { q = -2; - else if (pdg == 1000020040) + } + else if (pdg == 1000020040) { q = 2; - else if (pdg == -1000020040) + } + else if (pdg == -1000020040) { q = -2; - else + } + else { q = 0; + } Double_t p = cbmMCTrack->GetP(); mcTracks[iMC].SetMotherId(cbmMCTrack->GetMotherId()); @@ -416,7 +444,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) for (int iTr = 0; iTr < ntrackMatches; iTr++) { CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fStsTrackMatchArray->At(iTr); - if (stsTrackMatch->GetNofLinks() == 0) continue; + if (stsTrackMatch->GetNofLinks() == 0) { + continue; + } Float_t bestWeight = 0.f; Float_t totalWeight = 0.f; Int_t mcTrackId = -1; @@ -427,7 +457,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) mcTrackId = stsTrackMatch->GetLink(iLink).GetIndex(); } } - if (bestWeight / totalWeight < 0.7) continue; + if (bestWeight / totalWeight < 0.7) { + continue; + } if (mcTrackId >= nMCTracks || mcTrackId < 0) { std::cout << "Sts Matching is wrong! StsTackId = " << mcTrackId << " N mc tracks = " << nMCTracks << std::endl; continue; @@ -443,8 +475,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) for (int iTr = 0; iTr < fStsTrackArray->GetEntriesFast(); iTr++) { CbmStsTrack* stsTrack = ((CbmStsTrack*) fStsTrackArray->At(iTr)); vRTracks[iTr] = *stsTrack; - - if (trackMatch[iTr] > -1) pdg[iTr] = mcTracks[trackMatch[iTr]].PDG(); + if (trackMatch[iTr] > -1) { + pdg[iTr] = mcTracks[trackMatch[iTr]].PDG(); + } } CbmL1PFFitter fitter; @@ -456,10 +489,14 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) // fitter.GetChiToVertex(vRTracks, vField, vChiToPrimVtx, kfVertex, 3000000); for (unsigned int iTr = 0; iTr < vRTracks.size(); iTr++) { - if (trackMatch[iTr] < 0) continue; + if (trackMatch[iTr] < 0) { + continue; + } const KFMCTrack& mcTrack = mcTracks[trackMatch[iTr]]; - if (mcTrack.MotherId() > -1) continue; + if (mcTrack.MotherId() > -1) { + continue; + } // if ( vRTracks[iTr].GetNofHits() < 11 ) continue; const FairTrackParam* parameters = vRTracks[iTr].GetParamFirst(); @@ -507,7 +544,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) for (int iTr = 0; iTr < nMuchTrackMatches; iTr++) { CbmTrackMatchNew* muchTrackMatch = (CbmTrackMatchNew*) fMuchTrackMatchArray->At(iTr); - if (muchTrackMatch->GetNofLinks() == 0) continue; + if (muchTrackMatch->GetNofLinks() == 0) { + continue; + } Float_t bestWeight = 0.f; Float_t totalWeight = 0.f; Int_t mcTrackId = -1; @@ -518,7 +557,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) mcTrackId = muchTrackMatch->GetLink(iLink).GetIndex(); } } - if (bestWeight / totalWeight < 0.7) continue; + if (bestWeight / totalWeight < 0.7) { + continue; + } if (mcTrackId >= nMCTracks || mcTrackId < 0) { std::cout << "Much Matching is wrong! MuchTackId = " << mcTrackId << " N mc tracks = " << nMCTracks << std::endl; @@ -529,8 +570,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) } } - if (fGlobalTrackArray == nullptr) + if (fGlobalTrackArray == nullptr) { Warning("KF Track QA", "No GlobalTrack array!"); + } else { for (Int_t igt = 0; igt < fGlobalTrackArray->GetEntriesFast(); igt++) { const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTrackArray->At(igt)); @@ -575,16 +617,22 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) muchHistoData[3] = muchTrack->GetChiSq() / muchTrack->GetNDF(); //Chi2/NDF muchHistoData[4] = TMath::Prob(muchTrack->GetChiSq(), muchTrack->GetNDF()); //prob - if (stsTrackMCIndex < 0 || stsTrackMCIndex != muchTrackMCIndex) //ghost - for (int iH = 0; iH < NMuchHisto; iH++) + if (stsTrackMCIndex < 0 || stsTrackMCIndex != muchTrackMCIndex) { //ghost + for (int iH = 0; iH < NMuchHisto; iH++) { hMuchHisto[2][iH]->Fill(muchHistoData[iH]); + } + } else { - if (TMath::Abs(mcTracks[stsTrackMCIndex].PDG()) == 13) //muon - for (int iH = 0; iH < NMuchHisto; iH++) + if (TMath::Abs(mcTracks[stsTrackMCIndex].PDG()) == 13) { //muon + for (int iH = 0; iH < NMuchHisto; iH++) { hMuchHisto[0][iH]->Fill(muchHistoData[iH]); - else //BG - for (int iH = 0; iH < NMuchHisto; iH++) + } + } + else { //BG + for (int iH = 0; iH < NMuchHisto; iH++) { hMuchHisto[1][iH]->Fill(muchHistoData[iH]); + } + } } delete[] muchHistoData; } @@ -617,7 +665,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) bestMCTrackId = richRingMatch->GetLink(iLink).GetIndex(); } } - if (bestWeight / totalWeight >= 0.7) richTrackMCIndex = bestMCTrackId; + if (bestWeight / totalWeight >= 0.7) { + richTrackMCIndex = bestMCTrackId; + } } } @@ -630,14 +680,18 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) hRichRingHisto2D[0][2]->Fill(p, axisB); int iTrackCategory = -1; - if (stsTrackMCIndex < 0) + if (stsTrackMCIndex < 0) { iTrackCategory = 8; // ghost sts track + any ring - else if (richTrackMCIndex < 0) + } + else if (richTrackMCIndex < 0) { iTrackCategory = 9; // normal sts track + ghost ring - else if (stsTrackMCIndex != richTrackMCIndex) + } + else if (stsTrackMCIndex != richTrackMCIndex) { iTrackCategory = 7; // mismatched sts track and ring - else + } + else { iTrackCategory = GetHistoIndex(pdg[stsTrackIndex]); + } if (iTrackCategory < 10) { hRichRingHisto2D[iTrackCategory][0]->Fill(p, r); @@ -664,7 +718,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) Int_t trdTrackMCIndex = -1; CbmTrackMatchNew* trdTrackMatch = (CbmTrackMatchNew*) fTrdTrackMatchArray->At(trdIndex); for (int iTr = 0; iTr < ntrackTrdMatches; iTr++) { - if (trdTrackMatch->GetNofLinks() == 0) continue; + if (trdTrackMatch->GetNofLinks() == 0) { + continue; + } Float_t bestWeight = 0.f; Float_t totalWeight = 0.f; for (int iLink = 0; iLink < trdTrackMatch->GetNofLinks(); iLink++) { @@ -674,7 +730,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) trdTrackMCIndex = trdTrackMatch->GetLink(iLink).GetIndex(); } } - if (bestWeight / totalWeight < 0.7) continue; + if (bestWeight / totalWeight < 0.7) { + continue; + } if (trdTrackMCIndex >= nMCTracks || trdTrackMCIndex < 0) { std::cout << "Trd Matching is wrong! TrdTackId = " << trdTrackMCIndex << " N mc tracks = " << nMCTracks << std::endl; @@ -695,17 +753,22 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) CbmTrdHit* trdHit = (CbmTrdHit*) fTrdHitArray->At(TRDindex); eloss += trdHit->GetELoss(); } - if (trdTrack->GetNofHits() > 0.) eloss = eloss / trdTrack->GetNofHits(); // average of dE/dx per station - + if (trdTrack->GetNofHits() > 0.) { + eloss = eloss / trdTrack->GetNofHits(); // average of dE/dx per station + } int iTrackCategory = -1; - if (stsTrackMCIndex < 0) + if (stsTrackMCIndex < 0) { iTrackCategory = 8; // ghost sts track + any trd track - else if (trdTrackMCIndex < 0) + } + else if (trdTrackMCIndex < 0) { iTrackCategory = 9; // normal sts track + ghost trd track - else if (stsTrackMCIndex != trdTrackMCIndex) + } + else if (stsTrackMCIndex != trdTrackMCIndex) { iTrackCategory = 7; // mismatched sts track and trd track - else + } + else { iTrackCategory = GetHistoIndex(pdg[stsTrackIndex]); + } for (int iH = 0; iH < NTrdHisto; iH++) { hTrdHisto[0][iH]->Fill(trdHistoData[iH]); @@ -745,7 +808,9 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) for (Int_t ih = 0; ih < fTofHitMatchArray->GetEntriesFast(); ih++) { CbmMatch* tofHitMatch = (CbmMatch*) fTofHitMatchArray->At(ih); CbmTofPoint* tofPoint = (CbmTofPoint*) fTofPoints->Get(tofHitMatch->GetMatchedLink()); - if (!tofPoint) continue; + if (!tofPoint) { + continue; + } Int_t tofMCTrackId = tofPoint->GetTrackID(); if (tofMCTrackId == stsTrackMCIndex) { isReconstructible = 1; @@ -828,14 +893,18 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/) // std::cout << " mcl: " << tofMCl << " mc time: " << tofMCTime << " mc beta: " << tofMCl/tofMCTime/29.9792458 << std::endl; int iHitCategory = -1; - if (stsTrackMCIndex < 0) + if (stsTrackMCIndex < 0) { iHitCategory = 8; // ghost sts track + any tof hit - else if (tofMCTrackId < 0) + } + else if (tofMCTrackId < 0) { iHitCategory = 9; // normal sts track + ghost tof hit - else if (stsTrackMCIndex != tofMCTrackId) + } + else if (stsTrackMCIndex != tofMCTrackId) { iHitCategory = 7; // mismatched sts track and tof hit - else + } + else { iHitCategory = GetHistoIndex(pdg[stsTrackIndex]); + } hTofHisto2D[iHitCategory][0]->Fill(p * q, m2); hTofHisto2D[iHitCategory][1]->Fill(eloss * 1e6, m2); @@ -865,8 +934,9 @@ void CbmKFTrackQa::WriteHistosCurFile(TObject* obj) { - if (!obj->IsFolder()) + if (!obj->IsFolder()) { obj->Write(); + } else { TDirectory* cur = gDirectory; TFile* currentFile = gFile; @@ -887,32 +957,70 @@ int CbmKFTrackQa::GetHistoIndex(int pdg) { map<int, int>::iterator it; it = fPDGtoIndexMap.find(TMath::Abs(pdg)); - if (it != fPDGtoIndexMap.end()) + if (it != fPDGtoIndexMap.end()) { return it->second; - else + } + else { return 6; + } } Int_t CbmKFTrackQa::GetZtoNStation(Double_t getZ) { - if (TMath::Abs(getZ - 145) <= 2.0) return 1; - if (TMath::Abs(getZ - 155) <= 2.0) return 2; - if (TMath::Abs(getZ - 165) <= 2.0) return 3; - if (TMath::Abs(getZ - 195) <= 2.0) return 4; - if (TMath::Abs(getZ - 205) <= 2.0) return 5; - if (TMath::Abs(getZ - 215) <= 2.0) return 6; - if (TMath::Abs(getZ - 245) <= 2.0) return 7; - if (TMath::Abs(getZ - 255) <= 2.0) return 8; - if (TMath::Abs(getZ - 265) <= 2.0) return 9; - if (TMath::Abs(getZ - 305) <= 2.0) return 10; - if (TMath::Abs(getZ - 315) <= 2.0) return 11; - if (TMath::Abs(getZ - 325) <= 2.0) return 12; - if (TMath::Abs(getZ - 370) <= 2.0) return 13; - if (TMath::Abs(getZ - 380) <= 2.0) return 14; - if (TMath::Abs(getZ - 390) <= 2.0) return 15; - if (TMath::Abs(getZ - 500) <= 2.0) return 16; - if (TMath::Abs(getZ - 510) <= 2.0) return 17; - if (TMath::Abs(getZ - 520) <= 2.0) return 18; + if (TMath::Abs(getZ - 145) <= 2.0) { + return 1; + } + if (TMath::Abs(getZ - 155) <= 2.0) { + return 2; + } + if (TMath::Abs(getZ - 165) <= 2.0) { + return 3; + } + if (TMath::Abs(getZ - 195) <= 2.0) { + return 4; + } + if (TMath::Abs(getZ - 205) <= 2.0) { + return 5; + } + if (TMath::Abs(getZ - 215) <= 2.0) { + return 6; + } + if (TMath::Abs(getZ - 245) <= 2.0) { + return 7; + } + if (TMath::Abs(getZ - 255) <= 2.0) { + return 8; + } + if (TMath::Abs(getZ - 265) <= 2.0) { + return 9; + } + if (TMath::Abs(getZ - 305) <= 2.0) { + return 10; + } + if (TMath::Abs(getZ - 315) <= 2.0) { + return 11; + } + if (TMath::Abs(getZ - 370) <= 2.0) { + return 13; + } + if (TMath::Abs(getZ - 325) <= 2.0) { + return 12; + } + if (TMath::Abs(getZ - 380) <= 2.0) { + return 14; + } + if (TMath::Abs(getZ - 390) <= 2.0) { + return 15; + } + if (TMath::Abs(getZ - 500) <= 2.0) { + return 16; + } + if (TMath::Abs(getZ - 510) <= 2.0) { + return 17; + } + if (TMath::Abs(getZ - 520) <= 2.0) { + return 18; + } return -1; } diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx index 1d2f7ee5f845994c2c90629f74bdc53b217ab8c7..d29c56df12530c51783e6f74e55ccc882864b068 100644 --- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx +++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx @@ -99,7 +99,9 @@ CbmL1PFFitter::~CbmL1PFFitter() {} inline void CbmL1PFFitter::Initialize() { - if (fIsInitialised) return; + if (fIsInitialised) { + return; + } fNmvdStationsActive = 0; fNstsStationsActive = 0; @@ -214,7 +216,9 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM T.Vi() = 0.; T.ResetErrors(1.e2, 1.e2, 1., 1., 1., 1.e6, 1.e2); - if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack; + if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) { + nTracks_SIMD = N_vTracks - itrack; + } for (int iVec = 0; iVec < nTracks_SIMD; iVec++) { tr[iVec] = &Tracks[itrack + iVec]; // current track T.X()[iVec] = tr[iVec]->GetParamFirst()->GetX(); @@ -233,7 +237,9 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM } int pid = pidHypo[itrack + iVec]; - if (pid == -1) pid = 211; + if (pid == -1) { + pid = 211; + } // mass[i] = TDatabasePDG::Instance()->GetParticle(pid)->Mass(); mass[iVec] = KFParticleDatabase::Instance()->GetMass(pid); } @@ -271,14 +277,18 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM int hitIndex = tr[iVec]->GetMvdHitIndex(i); const CbmMvdHit* mvdHit = &(vMvdHits[hitIndex]); ista = GetMvdStationIndex(mvdHit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } hit = mvdHit; } else { int hitIndex = tr[iVec]->GetStsHitIndex(i - nHitsTrackMvd); const CbmStsHit* stsHit = &(vStsHits[hitIndex]); ista = GetStsStationIndex(stsHit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } hit = stsHit; } @@ -499,7 +509,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe int nStations = fNmvdStationsActive + fNstsStationsActive; const ca::StationV* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin(); - fvec* zSta = new fvec[nStations]; + fvec* zSta = new fvec[nStations]; for (int iSta = 0; iSta < nStations; iSta++) { zSta[iSta] = sta[iSta].fZ; } @@ -542,7 +552,9 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe float posx = 0.f, posy = 0.f, posz = 0.f; if (iH < nHitsTrackMvd) { - if (!fMvdHitArray) continue; + if (!fMvdHitArray) { + continue; + } int hitIndex = tr[iVec]->GetMvdHitIndex(iH); const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex)); @@ -550,10 +562,14 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe posy = hit->GetY(); posz = hit->GetZ(); ista = GetMvdStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } else { - if (!fStsHitArray) continue; + if (!fStsHitArray) { + continue; + } int hitIndex = tr[iVec]->GetStsHitIndex(iH - nHitsTrackMvd); const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex)); @@ -561,7 +577,9 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe posy = hit->GetY(); posz = hit->GetZ(); ista = GetStsStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy); @@ -659,10 +677,13 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF unsigned short N_vTracks = Tracks.size(); for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) { - if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack; + if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) { + nTracks_SIMD = N_vTracks - itrack; + } - for (int i = 0; i < nTracks_SIMD; i++) + for (int i = 0; i < nTracks_SIMD; i++) { tr[i] = &Tracks[itrack + i]; // current track + } for (int iVec = 0; iVec < nTracks_SIMD; iVec++) { int nHitsTrackMvd = tr[iVec]->GetNofMvdHits(); @@ -670,7 +691,9 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF float posx = 0.f, posy = 0.f, posz = 0.f; if (iH < nHitsTrackMvd) { - if (!fMvdHitArray) continue; + if (!fMvdHitArray) { + continue; + } int hitIndex = tr[iVec]->GetMvdHitIndex(iH); const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex)); @@ -678,10 +701,14 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF posy = hit->GetY(); posz = hit->GetZ(); ista = GetMvdStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } else { - if (!fStsHitArray) continue; + if (!fStsHitArray) { + continue; + } int hitIndex = tr[iVec]->GetStsHitIndex(iH - nHitsTrackMvd); const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex)); @@ -689,7 +716,9 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF posy = hit->GetY(); posz = hit->GetZ(); ista = GetStsStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy); @@ -730,7 +759,9 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, unsigned short N_vTracks = Tracks.size(); for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) { - if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack; + if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) { + nTracks_SIMD = N_vTracks - itrack; + } for (int i = 0; i < nTracks_SIMD; i++) { tr[i] = &Tracks[itrack + i]; // current track @@ -744,7 +775,9 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, int hitNumber = nHits - iH - 1; if (hitNumber < nHitsTrackMvd) { - if (!fMvdHitArray) continue; + if (!fMvdHitArray) { + continue; + } int hitIndex = tr[iVec]->GetMvdHitIndex(hitNumber); const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex)); @@ -752,10 +785,14 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, posy = hit->GetY(); posz = hit->GetZ(); ista = GetMvdStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } else { - if (!fStsHitArray) continue; + if (!fStsHitArray) { + continue; + } int hitIndex = tr[iVec]->GetStsHitIndex(hitNumber - nHitsTrackMvd); const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex)); @@ -763,7 +800,9 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, posy = hit->GetY(); posz = hit->GetZ(); ista = GetStsStationIndex(hit); - if (ista < 0) continue; + if (ista < 0) { + continue; + } } fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);