Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:02:15 

0001 //
0002 //  Simple photon conversion seeding class (src)
0003 //
0004 //  Author: E Song
0005 //
0006 //  Version: 1;     6 Aug 2012
0007 //
0008 
0009 #define Conv4HitsReco2_cxx
0010 
0011 //#include "Conv4HitsReco2.h"
0012 //#include "FWCore/MessegeLogger/interface/MessegeLogger.h"
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 //Conv4HitsReco2::~Conv4HitsReco() { }
0023 
0024 void Conv4HitsReco2::Refresh(
0025     math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4) {
0026   // Fix 2D plane, make primary vertex the original point
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   // DEFAULT setup
0043   fMaxNumberOfIterations = 40;
0044   fFixedNumberOfIterations = 0;
0045   fRadiusECut = 10.0;  //cm
0046   fPhiECut = 0.03;     //rad
0047   fRECut = 0.5;        //cm
0048   fBField = 3.8;       //T
0049 
0050   // TRIVIAL initialization
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;   // 2 - positron
0113   ptminus = fRecR1 * fBField * 0.01 * 0.3;  // 1 - electron
0114   vtx = fRecV;
0115   //std::cout << ".";
0116   return fLoop;
0117 }
0118 
0119 void Conv4HitsReco2::Reconstruct() {
0120   double x11, x12, x21, x22, y11, y12, y21, y22;
0121   //double X11,X12,X21,X22,Y11,Y12,Y21,Y22;
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   // Setting Phi1, Phi2 initial guess range, and first guess
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;  // stupid Atan correction
0145 
0146   fPhiE = std::fabs((Phi1 - Phi2)) / std::pow(2.0, fLoop + 1);
0147 
0148   double NextPhi = (Phi1 + Phi2) / 2.0;  // first guess
0149   double D1, D2 = 0.0;
0150   double prevR1 = 0;
0151   double prevR2 = 0;
0152   double R1 = 0;
0153   double R2 = 0;
0154 
0155   // Iterations
0156   for (int i = 0; i < fLoop; i++) {
0157     // LOCAL TRANFORMATION & EXTRACTION
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     // I'm not using LocalTransform() function because this direct way turns out to be faster
0177 
0178     // SOLVING EQUATIONS
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     }  // No real root.  Break out.
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)) {  // Moving RIGHT
0200         Phi1 = NextPhi;
0201         // Phi2 remains the same
0202         NextPhi = (Phi1 + Phi2) / 2.0;
0203       } else if ((Y11 - D1) < (Y12 - D2)) {  // Moving LEFT
0204         // Phi1 remains the same
0205         Phi2 = NextPhi;
0206         NextPhi = (Phi1 + Phi2) / 2.0;
0207       }
0208 
0209       // CHECK STOP CONDITION
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       // A. Cut threshold satisfied - STOP - record
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       // B. Cut threshold NOT satisfied - prepare for next loop
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!