Reference Manual

Go to the documentation of this file.
00001 //
00002 //      Simple photon conversion seeding class (src)
00003 //
00004 //      Author: E Song
00005 //
00006 //      Version: 1;     6 Aug 2012
00007 //
00008 
00009 #define Conv4HitsReco2_cxx
00010 
00011 //#include "RecoTracker/ConversionSeedGenerators/interface/Conv4HitsReco2.h"
00012 //#include "FWCore/MessegeLogger/interface/MessegeLogger.h"
00013 #include "RecoTracker/ConversionSeedGenerators/interface/Conv4HitsReco2.h"
00014 #include <time.h>
00015 
00016 Conv4HitsReco2::Conv4HitsReco2(math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4)
00017 {
00018         Refresh(vPhotVertex, h1, h2, h3, h4);
00019 }
00020 
00021 Conv4HitsReco2::~Conv4HitsReco2() { }
00022 //Conv4HitsReco2::~Conv4HitsReco() { }
00023 
00024 void Conv4HitsReco2::Refresh(math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4)
00025 {
00026         // Fix 2D plane, make primary vertex the original point
00027         fPV = vPhotVertex;      fPV.SetZ(0.);
00028         fHitv11 = h3;           fHitv11.SetZ(0.);       fHitv11 = fHitv11 - fPV;
00029         fHitv21 = h4;           fHitv21.SetZ(0.);       fHitv21 = fHitv21 - fPV;
00030         fHitv12 = h2;           fHitv12.SetZ(0.);       fHitv12 = fHitv12 - fPV;
00031         fHitv22 = h1;           fHitv22.SetZ(0.);       fHitv22 = fHitv22 - fPV;
00032 
00033         // DEFAULT setup
00034         fMaxNumberOfIterations = 40;
00035         fFixedNumberOfIterations = 0;
00036         fRadiusECut = 10.0;//cm
00037         fPhiECut = 0.03;//rad
00038         fRECut = 0.5;//cm       
00039         fBField = 3.8;//T       
00040 
00041         // TRIVIAL initialization
00042         fCutSatisfied = 0;
00043         fSignSatisfied = 0;
00044         fSolved = 0;
00045 
00046         fRecPhi = 0.;
00047         fRecR = 0.;
00048         fRecR1 = 0.;
00049         fRecR2 = 0.;
00050         fRadiusE = 0.;
00051         fRE = 0.;
00052         fPhiE = 0.;
00053 }
00054 
00055 void Conv4HitsReco2::LocalTransformation(math::XYZVector v11, math::XYZVector v12, math::XYZVector v21, math::XYZVector v22,
00056                                          math::XYZVector &V11, math::XYZVector &V12, math::XYZVector &V21, math::XYZVector &V22,
00057                                          double NextPhi) 
00058 {
00059         double x11,x12,x21,x22,y11,y12,y21,y22;
00060 
00061         x11 = v11.X(); y11 = v11.Y();
00062         x12 = v12.X(); y12 = v12.Y();
00063         x21 = v21.X(); y21 = v21.Y();
00064         x22 = v22.X(); y22 = v22.Y();
00065 
00066         double SINP = std::sin(NextPhi);
00067         double COSP = std::cos(NextPhi);
00068         double SignCOSP = 1.; if(COSP < 0.) SignCOSP = -1.;
00069         double AbsCOSP = std::fabs(COSP);
00070 
00071         double X11 = -std::fabs(x11*SINP*SignCOSP - y11*AbsCOSP);
00072         double Y11 =  std::fabs(y11*SINP*SignCOSP + x11*AbsCOSP);
00073 
00074         double X21 = -std::fabs(x21*SINP*SignCOSP - y21*AbsCOSP);
00075         double Y21 =  std::fabs(y21*SINP*SignCOSP + x21*AbsCOSP);
00076 
00077         double X12 =  std::fabs(x12*SINP*SignCOSP - y12*AbsCOSP);
00078         double Y12 =  std::fabs(y12*SINP*SignCOSP + x12*AbsCOSP);
00079 
00080         double X22 =  std::fabs(x22*SINP*SignCOSP - y22*AbsCOSP);
00081         double Y22 =  std::fabs(y22*SINP*SignCOSP + x22*AbsCOSP);
00082 
00083         V11.SetXYZ(X11,Y11,0.);
00084         V12.SetXYZ(X12,Y12,0.);
00085         V21.SetXYZ(X21,Y21,0.);
00086         V22.SetXYZ(X22,Y22,0.);
00087 }
00088 
00089 
00090 int Conv4HitsReco2::ConversionCandidate(math::XYZVector &vtx, double &ptplus, double &ptminus)
00091 {
00092         Reconstruct();
00093 
00094         ptplus = fRecR2 * fBField * 0.01 * 0.3;         // 2 - positron
00095         ptminus = fRecR1 * fBField * 0.01 * 0.3;        // 1 - electron
00096         vtx = fRecV;
00097         //std::cout << ".";
00098         return fLoop;
00099 }
00100 
00101 
00102 void Conv4HitsReco2::Reconstruct()
00103 {
00104         double x11,x12,x21,x22,y11,y12,y21,y22;
00105         //double X11,X12,X21,X22,Y11,Y12,Y21,Y22;
00106         x11 = fHitv11.X();      y11 = fHitv11.Y();
00107         x12 = fHitv12.X();      y12 = fHitv12.Y();
00108         x21 = fHitv21.X();      y21 = fHitv21.Y();
00109         x22 = fHitv22.X();      y22 = fHitv22.Y(); 
00110 
00111         if (fFixedNumberOfIterations==0)
00112         fLoop = fMaxNumberOfIterations;
00113         else fLoop = fFixedNumberOfIterations;
00114 
00115         // Setting Phi1, Phi2 initial guess range, and first guess
00116         double tempr1 = std::sqrt(y11*y11 + x11*x11);
00117         double tempr2 = std::sqrt(y12*y12 + x12*x12);
00118 
00119         double Phi1 = 2.0 * std::atan(y11 / (x11+tempr1));
00120         double Phi2 = 2.0 * std::atan(y12 / (x12+tempr2));
00121 
00122         if (Phi1<Phi2) Phi1 += 2.0 * 3.141592653;       // stupid Atan correction
00123 
00124         fPhiE = std::fabs((Phi1-Phi2)) / std::pow(2.0, fLoop + 1);
00125 
00126         double NextPhi = ( Phi1 + Phi2 ) / 2.0; // first guess
00127         double D1, D2 = 0.0;
00128         double prevR1 = 0; double prevR2 = 0;
00129         double R1 = 0; double R2 = 0;
00130 
00131         // Iterations
00132         for (int i=0; i<fLoop; i++) {
00133 
00134                 // LOCAL TRANFORMATION & EXTRACTION
00135                 double SINP = std::sin(NextPhi);
00136                 double COSP = std::cos(NextPhi);
00137                 double SignCOSP = 1.; if(COSP < 0.) SignCOSP = -1.;
00138                 double AbsCOSP = std::fabs(COSP);
00139 
00140                 double X11 = -std::fabs(x11*SINP*SignCOSP - y11*AbsCOSP);
00141                 double Y11 =  std::fabs(y11*SINP*SignCOSP + x11*AbsCOSP);
00142 
00143                 double X21 = -std::fabs(x21*SINP*SignCOSP - y21*AbsCOSP);
00144                 double Y21 =  std::fabs(y21*SINP*SignCOSP + x21*AbsCOSP);
00145 
00146                 double X12 =  std::fabs(x12*SINP*SignCOSP - y12*AbsCOSP);
00147                 double Y12 =  std::fabs(y12*SINP*SignCOSP + x12*AbsCOSP);
00148 
00149                 double X22 =  std::fabs(x22*SINP*SignCOSP - y22*AbsCOSP);
00150                 double Y22 =  std::fabs(y22*SINP*SignCOSP + x22*AbsCOSP);
00151                 // I'm not using LocalTransform() function because this direct way turns out to be faster
00152 
00153 
00154 
00155 
00156                 // SOLVING EQUATIONS
00157                 double d1 = Y21 - Y11;
00158                 double d2 = Y22 - Y12;
00159 
00160                 if ( ( (X11*X11*d1*d1/(X21-X11)/(X21-X11) + X11*X21 + X11*d1*d1/(X21-X11)) < 0 ) || 
00161                          ( (X12*X12*d2*d2/(X22-X12)/(X22-X12) + X12*X22 + X12*d2*d2/(X22-X12)) < 0 )    )       
00162                         { fSolved = -1; fLoop = i;      return; } // No real root.  Break out.
00163         
00164                 else {
00165                         fSolved = 1;
00166                         D1 = X11*d1/(X21-X11);
00167                         D1 = D1 + std::sqrt(X11*X11*d1*d1/(X21-X11)/(X21-X11) + X11*X21 + X11*d1*d1/(X21-X11));
00168                         D2 = X12*d2/(X22-X12);
00169                         D2 = D2 + std::sqrt(X12*X12*d2*d2/(X22-X12)/(X22-X12) + X12*X22 + X12*d2*d2/(X22-X12));
00170         
00171                         R1 = std::fabs((X11+X21)/2.0+(D1+d1/2.0)*d1/(X21-X11));
00172                         R2 = std::fabs((X12+X22)/2.0+(D2+d2/2.0)*d2/(X22-X12));
00173 
00174                         if ((Y11-D1)>=(Y12-D2)) {  // Moving RIGHT
00175                                 Phi1 = NextPhi;
00176                                 Phi2 = Phi2;
00177                                 NextPhi = (Phi1+Phi2)/2.0;
00178                         }
00179                         else if ((Y11-D1)<(Y12-D2)) {  // Moving LEFT
00180                                 Phi1 = Phi1;
00181                                 Phi2 = NextPhi;
00182                                 NextPhi = (Phi1+Phi2)/2.0; 
00183                         } 
00184                         
00185                         // CHECK STOP CONDITION
00186                         double tmpPhiE = std::fabs(Phi1-Phi2);
00187                         double tmpRE = std::fabs( (Y11 - D1) - (Y12 - D2) );
00188                         double tmpRadiusE = ( std::fabs(R1-prevR1) + std::fabs(R2-prevR2) ) / 2.;
00189                         
00190                         // A. Cut threshold satisfied - STOP - record
00191                         if (( tmpPhiE <= fPhiECut ) && ( tmpRE <= fRECut ) && ( tmpRadiusE <= fRadiusECut ) && ( fFixedNumberOfIterations ==0 ))
00192                         {
00193                                 fSolved = 1;  
00194                                 fCutSatisfied = 1; 
00195                                 fLoop = i+1;
00196                                 fPhiE = tmpPhiE; 
00197                                 fRE = tmpRE;    
00198                                 fRadiusE = tmpRadiusE;                          
00199                                 fRecR1 = R1; 
00200                                 fRecR2 = R2; 
00201                                 fRecR = ( (Y11 - D1) + (Y12 - D2) ) / 2.0; 
00202                                 fRecPhi = NextPhi;
00203                                 fRecV.SetX( fRecR * cos(fRecPhi) );
00204                                 fRecV.SetY( fRecR * sin(fRecPhi) );
00205                                 fRecC1.SetXYZ( fRecV.X()-fRecR1*sin(fRecPhi), fRecV.Y()+fRecR1*cos(fRecPhi), 0.);
00206                                 fRecC2.SetXYZ( fRecV.X()+fRecR2*sin(fRecPhi), fRecV.Y()-fRecR2*cos(fRecPhi), 0.);
00207                                 fRecV = fRecV + fPV;
00208                                 fRecC1 = fRecC1 + fPV;
00209                                 fRecC2 = fRecC2 + fPV;
00210                                 fCutSatisfied = 1;
00211                                 if ( (R1>0)&&(R2>0)&&(D1>0)&&(D2>0)&&((Y11-D1)>0)&&((Y12-D2)>0) ) fSignSatisfied = 1;
00212                                 else fSignSatisfied = 0;
00213                         }
00214                         else if (i==fLoop-1) {
00215                                 fSolved = 1;  
00216                                 fCutSatisfied = 1; 
00217                                 fLoop = i+1;
00218                                 fPhiE = tmpPhiE; 
00219                                 fRE = tmpRE;    
00220                                 fRadiusE = tmpRadiusE;                          
00221                                 fRecR1 = R1; 
00222                                 fRecR2 = R2; 
00223                                 fRecR = ( (Y11 - D1) + (Y12 - D2) ) / 2.0; 
00224                                 fRecPhi = NextPhi;
00225                                 fRecV.SetX( fRecR * cos(fRecPhi) );
00226                                 fRecV.SetY( fRecR * sin(fRecPhi) );
00227                                 fRecC1.SetXYZ( fRecV.X()-fRecR1*sin(fRecPhi), fRecV.Y()+fRecR1*cos(fRecPhi), 0.);
00228                                 fRecC2.SetXYZ( fRecV.X()+fRecR2*sin(fRecPhi), fRecV.Y()-fRecR2*cos(fRecPhi), 0.);
00229                                 fRecV = fRecV + fPV;
00230                                 fRecC1 = fRecC1 + fPV;
00231                                 fRecC2 = fRecC2 + fPV;
00232                                 fCutSatisfied = 0;
00233                                 if ( (R1>0)&&(R2>0)&&(D1>0)&&(D2>0)&&((Y11-D1)>0)&&((Y12-D2)>0) ) fSignSatisfied = 1;
00234                                 else fSignSatisfied = 0;
00235                         }
00236                         // B. Cut threshold NOT satisfied - prepare for next loop
00237                         prevR1 = R1;   prevR2 = R2; 
00238 
00239                 }
00240         }
00241 }
00242 
00243 
00244 void Conv4HitsReco2::Dump()
00245 {
00246         std::cout << std::endl<< "================================================" << std::endl;
00247         std::cout << "  Nothing happend here.";
00248         if (fSolved==1) std::cout << "Solved.";
00249         if (fCutSatisfied==1) std::cout << "Cut good.";
00250         if (fSignSatisfied==1) std::cout << "Sign good.";
00251         
00252 }
00253 
00254 math::XYZVector Conv4HitsReco2::GetPlusCenter(double &plusR)
00255 {
00256         plusR = fRecR1;
00257         return fRecC1;
00258         
00259 }
00260 math::XYZVector Conv4HitsReco2::GetMinusCenter(double &minusR)
00261 {
00262         minusR = fRecR2;
00263         return fRecC2;
00264 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoTracker/ConversionSeedGenerators/src/Conv4HitsReco2.cc
