File indexing completed on 2023-10-25 10:02:15
0001
0002
0003
0004
0005
0006
0007
0008
0009 #define Conv4HitsReco2_cxx
0010
0011
0012
0013 #include "Conv4HitsReco2.h"
0014 #include <ctime>
0015
0016 Conv4HitsReco2::Conv4HitsReco2(
0017 math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4) {
0018 Refresh(vPhotVertex, h1, h2, h3, h4);
0019 }
0020
0021 Conv4HitsReco2::~Conv4HitsReco2() {}
0022
0023
0024 void Conv4HitsReco2::Refresh(
0025 math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4) {
0026
0027 fPV = vPhotVertex;
0028 fPV.SetZ(0.);
0029 fHitv11 = h3;
0030 fHitv11.SetZ(0.);
0031 fHitv11 = fHitv11 - fPV;
0032 fHitv21 = h4;
0033 fHitv21.SetZ(0.);
0034 fHitv21 = fHitv21 - fPV;
0035 fHitv12 = h2;
0036 fHitv12.SetZ(0.);
0037 fHitv12 = fHitv12 - fPV;
0038 fHitv22 = h1;
0039 fHitv22.SetZ(0.);
0040 fHitv22 = fHitv22 - fPV;
0041
0042
0043 fMaxNumberOfIterations = 40;
0044 fFixedNumberOfIterations = 0;
0045 fRadiusECut = 10.0;
0046 fPhiECut = 0.03;
0047 fRECut = 0.5;
0048 fBField = 3.8;
0049
0050
0051 fCutSatisfied = 0;
0052 fSignSatisfied = 0;
0053 fSolved = 0;
0054
0055 fRecPhi = 0.;
0056 fRecR = 0.;
0057 fRecR1 = 0.;
0058 fRecR2 = 0.;
0059 fRadiusE = 0.;
0060 fRE = 0.;
0061 fPhiE = 0.;
0062 }
0063
0064 void Conv4HitsReco2::LocalTransformation(const math::XYZVector &v11,
0065 const math::XYZVector &v12,
0066 const math::XYZVector &v21,
0067 const math::XYZVector &v22,
0068 math::XYZVector &V11,
0069 math::XYZVector &V12,
0070 math::XYZVector &V21,
0071 math::XYZVector &V22,
0072 double NextPhi) {
0073 double x11, x12, x21, x22, y11, y12, y21, y22;
0074
0075 x11 = v11.X();
0076 y11 = v11.Y();
0077 x12 = v12.X();
0078 y12 = v12.Y();
0079 x21 = v21.X();
0080 y21 = v21.Y();
0081 x22 = v22.X();
0082 y22 = v22.Y();
0083
0084 double SINP = std::sin(NextPhi);
0085 double COSP = std::cos(NextPhi);
0086 double SignCOSP = 1.;
0087 if (COSP < 0.)
0088 SignCOSP = -1.;
0089 double AbsCOSP = std::fabs(COSP);
0090
0091 double X11 = -std::fabs(x11 * SINP * SignCOSP - y11 * AbsCOSP);
0092 double Y11 = std::fabs(y11 * SINP * SignCOSP + x11 * AbsCOSP);
0093
0094 double X21 = -std::fabs(x21 * SINP * SignCOSP - y21 * AbsCOSP);
0095 double Y21 = std::fabs(y21 * SINP * SignCOSP + x21 * AbsCOSP);
0096
0097 double X12 = std::fabs(x12 * SINP * SignCOSP - y12 * AbsCOSP);
0098 double Y12 = std::fabs(y12 * SINP * SignCOSP + x12 * AbsCOSP);
0099
0100 double X22 = std::fabs(x22 * SINP * SignCOSP - y22 * AbsCOSP);
0101 double Y22 = std::fabs(y22 * SINP * SignCOSP + x22 * AbsCOSP);
0102
0103 V11.SetXYZ(X11, Y11, 0.);
0104 V12.SetXYZ(X12, Y12, 0.);
0105 V21.SetXYZ(X21, Y21, 0.);
0106 V22.SetXYZ(X22, Y22, 0.);
0107 }
0108
0109 int Conv4HitsReco2::ConversionCandidate(math::XYZVector &vtx, double &ptplus, double &ptminus) {
0110 Reconstruct();
0111
0112 ptplus = fRecR2 * fBField * 0.01 * 0.3;
0113 ptminus = fRecR1 * fBField * 0.01 * 0.3;
0114 vtx = fRecV;
0115
0116 return fLoop;
0117 }
0118
0119 void Conv4HitsReco2::Reconstruct() {
0120 double x11, x12, x21, x22, y11, y12, y21, y22;
0121
0122 x11 = fHitv11.X();
0123 y11 = fHitv11.Y();
0124 x12 = fHitv12.X();
0125 y12 = fHitv12.Y();
0126 x21 = fHitv21.X();
0127 y21 = fHitv21.Y();
0128 x22 = fHitv22.X();
0129 y22 = fHitv22.Y();
0130
0131 if (fFixedNumberOfIterations == 0)
0132 fLoop = fMaxNumberOfIterations;
0133 else
0134 fLoop = fFixedNumberOfIterations;
0135
0136
0137 double tempr1 = std::sqrt(y11 * y11 + x11 * x11);
0138 double tempr2 = std::sqrt(y12 * y12 + x12 * x12);
0139
0140 double Phi1 = 2.0 * std::atan(y11 / (x11 + tempr1));
0141 double Phi2 = 2.0 * std::atan(y12 / (x12 + tempr2));
0142
0143 if (Phi1 < Phi2)
0144 Phi1 += 2.0 * 3.141592653;
0145
0146 fPhiE = std::fabs((Phi1 - Phi2)) / std::pow(2.0, fLoop + 1);
0147
0148 double NextPhi = (Phi1 + Phi2) / 2.0;
0149 double D1, D2 = 0.0;
0150 double prevR1 = 0;
0151 double prevR2 = 0;
0152 double R1 = 0;
0153 double R2 = 0;
0154
0155
0156 for (int i = 0; i < fLoop; i++) {
0157
0158 double SINP = std::sin(NextPhi);
0159 double COSP = std::cos(NextPhi);
0160 double SignCOSP = 1.;
0161 if (COSP < 0.)
0162 SignCOSP = -1.;
0163 double AbsCOSP = std::fabs(COSP);
0164
0165 double X11 = -std::fabs(x11 * SINP * SignCOSP - y11 * AbsCOSP);
0166 double Y11 = std::fabs(y11 * SINP * SignCOSP + x11 * AbsCOSP);
0167
0168 double X21 = -std::fabs(x21 * SINP * SignCOSP - y21 * AbsCOSP);
0169 double Y21 = std::fabs(y21 * SINP * SignCOSP + x21 * AbsCOSP);
0170
0171 double X12 = std::fabs(x12 * SINP * SignCOSP - y12 * AbsCOSP);
0172 double Y12 = std::fabs(y12 * SINP * SignCOSP + x12 * AbsCOSP);
0173
0174 double X22 = std::fabs(x22 * SINP * SignCOSP - y22 * AbsCOSP);
0175 double Y22 = std::fabs(y22 * SINP * SignCOSP + x22 * AbsCOSP);
0176
0177
0178
0179 double d1 = Y21 - Y11;
0180 double d2 = Y22 - Y12;
0181
0182 if (((X11 * X11 * d1 * d1 / (X21 - X11) / (X21 - X11) + X11 * X21 + X11 * d1 * d1 / (X21 - X11)) < 0) ||
0183 ((X12 * X12 * d2 * d2 / (X22 - X12) / (X22 - X12) + X12 * X22 + X12 * d2 * d2 / (X22 - X12)) < 0)) {
0184 fSolved = -1;
0185 fLoop = i;
0186 return;
0187 }
0188
0189 else {
0190 fSolved = 1;
0191 D1 = X11 * d1 / (X21 - X11);
0192 D1 = D1 + std::sqrt(X11 * X11 * d1 * d1 / (X21 - X11) / (X21 - X11) + X11 * X21 + X11 * d1 * d1 / (X21 - X11));
0193 D2 = X12 * d2 / (X22 - X12);
0194 D2 = D2 + std::sqrt(X12 * X12 * d2 * d2 / (X22 - X12) / (X22 - X12) + X12 * X22 + X12 * d2 * d2 / (X22 - X12));
0195
0196 R1 = std::fabs((X11 + X21) / 2.0 + (D1 + d1 / 2.0) * d1 / (X21 - X11));
0197 R2 = std::fabs((X12 + X22) / 2.0 + (D2 + d2 / 2.0) * d2 / (X22 - X12));
0198
0199 if ((Y11 - D1) >= (Y12 - D2)) {
0200 Phi1 = NextPhi;
0201
0202 NextPhi = (Phi1 + Phi2) / 2.0;
0203 } else if ((Y11 - D1) < (Y12 - D2)) {
0204
0205 Phi2 = NextPhi;
0206 NextPhi = (Phi1 + Phi2) / 2.0;
0207 }
0208
0209
0210 double tmpPhiE = std::fabs(Phi1 - Phi2);
0211 double tmpRE = std::fabs((Y11 - D1) - (Y12 - D2));
0212 double tmpRadiusE = (std::fabs(R1 - prevR1) + std::fabs(R2 - prevR2)) / 2.;
0213
0214
0215 if ((tmpPhiE <= fPhiECut) && (tmpRE <= fRECut) && (tmpRadiusE <= fRadiusECut) &&
0216 (fFixedNumberOfIterations == 0)) {
0217 fSolved = 1;
0218 fCutSatisfied = 1;
0219 fLoop = i + 1;
0220 fPhiE = tmpPhiE;
0221 fRE = tmpRE;
0222 fRadiusE = tmpRadiusE;
0223 fRecR1 = R1;
0224 fRecR2 = R2;
0225 fRecR = ((Y11 - D1) + (Y12 - D2)) / 2.0;
0226 fRecPhi = NextPhi;
0227 fRecV.SetX(fRecR * cos(fRecPhi));
0228 fRecV.SetY(fRecR * sin(fRecPhi));
0229 fRecC1.SetXYZ(fRecV.X() - fRecR1 * sin(fRecPhi), fRecV.Y() + fRecR1 * cos(fRecPhi), 0.);
0230 fRecC2.SetXYZ(fRecV.X() + fRecR2 * sin(fRecPhi), fRecV.Y() - fRecR2 * cos(fRecPhi), 0.);
0231 fRecV = fRecV + fPV;
0232 fRecC1 = fRecC1 + fPV;
0233 fRecC2 = fRecC2 + fPV;
0234 fCutSatisfied = 1;
0235 if ((R1 > 0) && (R2 > 0) && (D1 > 0) && (D2 > 0) && ((Y11 - D1) > 0) && ((Y12 - D2) > 0))
0236 fSignSatisfied = 1;
0237 else
0238 fSignSatisfied = 0;
0239 } else if (i == fLoop - 1) {
0240 fSolved = 1;
0241 fCutSatisfied = 1;
0242 fLoop = i + 1;
0243 fPhiE = tmpPhiE;
0244 fRE = tmpRE;
0245 fRadiusE = tmpRadiusE;
0246 fRecR1 = R1;
0247 fRecR2 = R2;
0248 fRecR = ((Y11 - D1) + (Y12 - D2)) / 2.0;
0249 fRecPhi = NextPhi;
0250 fRecV.SetX(fRecR * cos(fRecPhi));
0251 fRecV.SetY(fRecR * sin(fRecPhi));
0252 fRecC1.SetXYZ(fRecV.X() - fRecR1 * sin(fRecPhi), fRecV.Y() + fRecR1 * cos(fRecPhi), 0.);
0253 fRecC2.SetXYZ(fRecV.X() + fRecR2 * sin(fRecPhi), fRecV.Y() - fRecR2 * cos(fRecPhi), 0.);
0254 fRecV = fRecV + fPV;
0255 fRecC1 = fRecC1 + fPV;
0256 fRecC2 = fRecC2 + fPV;
0257 fCutSatisfied = 0;
0258 if ((R1 > 0) && (R2 > 0) && (D1 > 0) && (D2 > 0) && ((Y11 - D1) > 0) && ((Y12 - D2) > 0))
0259 fSignSatisfied = 1;
0260 else
0261 fSignSatisfied = 0;
0262 }
0263
0264 prevR1 = R1;
0265 prevR2 = R2;
0266 }
0267 }
0268 }
0269
0270 void Conv4HitsReco2::Dump() {
0271 std::cout << std::endl << "================================================" << std::endl;
0272 std::cout << " Nothing happend here.";
0273 if (fSolved == 1)
0274 std::cout << "Solved.";
0275 if (fCutSatisfied == 1)
0276 std::cout << "Cut good.";
0277 if (fSignSatisfied == 1)
0278 std::cout << "Sign good.";
0279 }
0280
0281 math::XYZVector Conv4HitsReco2::GetPlusCenter(double &plusR) {
0282 plusR = fRecR1;
0283 return fRecC1;
0284 }
0285 math::XYZVector Conv4HitsReco2::GetMinusCenter(double &minusR) {
0286 minusR = fRecR2;
0287 return fRecC2;
0288 }