Conv4HitsReco2.cc

Go to the documentation of this file.

//
//  Simple photon conversion seeding class (src)
//
//  Author: E Song
//
//  Version: 1;     6 Aug 2012
//

#define Conv4HitsReco2_cxx

//#include "Conv4HitsReco2.h"
//#include "FWCore/MessegeLogger/interface/MessegeLogger.h"
#include "Conv4HitsReco2.h"
#include <ctime>

Conv4HitsReco2::Conv4HitsReco2(math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4)
{
    Refresh(vPhotVertex, h1, h2, h3, h4);
}

Conv4HitsReco2::~Conv4HitsReco2() { }
//Conv4HitsReco2::~Conv4HitsReco() { }

void Conv4HitsReco2::Refresh(math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4)
{
    // Fix 2D plane, make primary vertex the original point
    fPV = vPhotVertex;  fPV.SetZ(0.);
    fHitv11 = h3;       fHitv11.SetZ(0.);   fHitv11 = fHitv11 - fPV;
    fHitv21 = h4;       fHitv21.SetZ(0.);   fHitv21 = fHitv21 - fPV;
    fHitv12 = h2;       fHitv12.SetZ(0.);   fHitv12 = fHitv12 - fPV;
    fHitv22 = h1;       fHitv22.SetZ(0.);   fHitv22 = fHitv22 - fPV;

    // DEFAULT setup
    fMaxNumberOfIterations = 40;
    fFixedNumberOfIterations = 0;
    fRadiusECut = 10.0;//cm
    fPhiECut = 0.03;//rad
    fRECut = 0.5;//cm   
    fBField = 3.8;//T   

    // TRIVIAL initialization
    fCutSatisfied = 0;
    fSignSatisfied = 0;
    fSolved = 0;

    fRecPhi = 0.;
    fRecR = 0.;
    fRecR1 = 0.;
    fRecR2 = 0.;
    fRadiusE = 0.;
    fRE = 0.;
    fPhiE = 0.;
}

void Conv4HitsReco2::LocalTransformation(const math::XYZVector& v11, const math::XYZVector& v12, const math::XYZVector& v21, const math::XYZVector& v22,
                     math::XYZVector &V11, math::XYZVector &V12, math::XYZVector &V21, math::XYZVector &V22,
                     double NextPhi) 
{
    double x11,x12,x21,x22,y11,y12,y21,y22;

    x11 = v11.X(); y11 = v11.Y();
    x12 = v12.X(); y12 = v12.Y();
    x21 = v21.X(); y21 = v21.Y();
    x22 = v22.X(); y22 = v22.Y();

    double SINP = std::sin(NextPhi);
    double COSP = std::cos(NextPhi);
    double SignCOSP = 1.; if(COSP < 0.) SignCOSP = -1.;
    double AbsCOSP = std::fabs(COSP);

    double X11 = -std::fabs(x11*SINP*SignCOSP - y11*AbsCOSP);
    double Y11 =  std::fabs(y11*SINP*SignCOSP + x11*AbsCOSP);

    double X21 = -std::fabs(x21*SINP*SignCOSP - y21*AbsCOSP);
    double Y21 =  std::fabs(y21*SINP*SignCOSP + x21*AbsCOSP);

    double X12 =  std::fabs(x12*SINP*SignCOSP - y12*AbsCOSP);
    double Y12 =  std::fabs(y12*SINP*SignCOSP + x12*AbsCOSP);

    double X22 =  std::fabs(x22*SINP*SignCOSP - y22*AbsCOSP);
    double Y22 =  std::fabs(y22*SINP*SignCOSP + x22*AbsCOSP);

    V11.SetXYZ(X11,Y11,0.);
    V12.SetXYZ(X12,Y12,0.);
    V21.SetXYZ(X21,Y21,0.);
    V22.SetXYZ(X22,Y22,0.);
}


int Conv4HitsReco2::ConversionCandidate(math::XYZVector &vtx, double &ptplus, double &ptminus)
{
    Reconstruct();

    ptplus = fRecR2 * fBField * 0.01 * 0.3;     // 2 - positron
    ptminus = fRecR1 * fBField * 0.01 * 0.3;    // 1 - electron
    vtx = fRecV;
    //std::cout << ".";
    return fLoop;
}


void Conv4HitsReco2::Reconstruct()
{
    double x11,x12,x21,x22,y11,y12,y21,y22;
    //double X11,X12,X21,X22,Y11,Y12,Y21,Y22;
    x11 = fHitv11.X();  y11 = fHitv11.Y();
    x12 = fHitv12.X();  y12 = fHitv12.Y();
    x21 = fHitv21.X();  y21 = fHitv21.Y();
    x22 = fHitv22.X();  y22 = fHitv22.Y(); 

    if (fFixedNumberOfIterations==0)
    fLoop = fMaxNumberOfIterations;
    else fLoop = fFixedNumberOfIterations;

    // Setting Phi1, Phi2 initial guess range, and first guess
    double tempr1 = std::sqrt(y11*y11 + x11*x11);
    double tempr2 = std::sqrt(y12*y12 + x12*x12);

    double Phi1 = 2.0 * std::atan(y11 / (x11+tempr1));
    double Phi2 = 2.0 * std::atan(y12 / (x12+tempr2));

    if (Phi1<Phi2) Phi1 += 2.0 * 3.141592653;   // stupid Atan correction

    fPhiE = std::fabs((Phi1-Phi2)) / std::pow(2.0, fLoop + 1);

    double NextPhi = ( Phi1 + Phi2 ) / 2.0; // first guess
    double D1, D2 = 0.0;
    double prevR1 = 0; double prevR2 = 0;
    double R1 = 0; double R2 = 0;

    // Iterations
    for (int i=0; i<fLoop; i++) {

        // LOCAL TRANFORMATION & EXTRACTION
        double SINP = std::sin(NextPhi);
        double COSP = std::cos(NextPhi);
        double SignCOSP = 1.; if(COSP < 0.) SignCOSP = -1.;
        double AbsCOSP = std::fabs(COSP);

        double X11 = -std::fabs(x11*SINP*SignCOSP - y11*AbsCOSP);
        double Y11 =  std::fabs(y11*SINP*SignCOSP + x11*AbsCOSP);

        double X21 = -std::fabs(x21*SINP*SignCOSP - y21*AbsCOSP);
        double Y21 =  std::fabs(y21*SINP*SignCOSP + x21*AbsCOSP);

        double X12 =  std::fabs(x12*SINP*SignCOSP - y12*AbsCOSP);
        double Y12 =  std::fabs(y12*SINP*SignCOSP + x12*AbsCOSP);

        double X22 =  std::fabs(x22*SINP*SignCOSP - y22*AbsCOSP);
        double Y22 =  std::fabs(y22*SINP*SignCOSP + x22*AbsCOSP);
        // I'm not using LocalTransform() function because this direct way turns out to be faster




        // SOLVING EQUATIONS
        double d1 = Y21 - Y11;
        double d2 = Y22 - Y12;

        if ( ( (X11*X11*d1*d1/(X21-X11)/(X21-X11) + X11*X21 + X11*d1*d1/(X21-X11)) < 0 ) || 
             ( (X12*X12*d2*d2/(X22-X12)/(X22-X12) + X12*X22 + X12*d2*d2/(X22-X12)) < 0 )    )   
            { fSolved = -1; fLoop = i;  return; } // No real root.  Break out.
    
        else {
            fSolved = 1;
            D1 = X11*d1/(X21-X11);
            D1 = D1 + std::sqrt(X11*X11*d1*d1/(X21-X11)/(X21-X11) + X11*X21 + X11*d1*d1/(X21-X11));
            D2 = X12*d2/(X22-X12);
            D2 = D2 + std::sqrt(X12*X12*d2*d2/(X22-X12)/(X22-X12) + X12*X22 + X12*d2*d2/(X22-X12));
    
            R1 = std::fabs((X11+X21)/2.0+(D1+d1/2.0)*d1/(X21-X11));
            R2 = std::fabs((X12+X22)/2.0+(D2+d2/2.0)*d2/(X22-X12));

            if ((Y11-D1)>=(Y12-D2)) {  // Moving RIGHT
                Phi1 = NextPhi;
                                // Phi2 remains the same
                NextPhi = (Phi1+Phi2)/2.0;
            }
            else if ((Y11-D1)<(Y12-D2)) {  // Moving LEFT
                // Phi1 remains the same
                Phi2 = NextPhi;
                NextPhi = (Phi1+Phi2)/2.0; 
            } 
            
            // CHECK STOP CONDITION
            double tmpPhiE = std::fabs(Phi1-Phi2);
            double tmpRE = std::fabs( (Y11 - D1) - (Y12 - D2) );
            double tmpRadiusE = ( std::fabs(R1-prevR1) + std::fabs(R2-prevR2) ) / 2.;
            
            // A. Cut threshold satisfied - STOP - record
            if (( tmpPhiE <= fPhiECut ) && ( tmpRE <= fRECut ) && ( tmpRadiusE <= fRadiusECut ) && ( fFixedNumberOfIterations ==0 ))
            {
                fSolved = 1;  
                fCutSatisfied = 1; 
                fLoop = i+1;
                fPhiE = tmpPhiE; 
                fRE = tmpRE;    
                fRadiusE = tmpRadiusE;              
                fRecR1 = R1; 
                fRecR2 = R2; 
                fRecR = ( (Y11 - D1) + (Y12 - D2) ) / 2.0; 
                fRecPhi = NextPhi;
                fRecV.SetX( fRecR * cos(fRecPhi) );
                fRecV.SetY( fRecR * sin(fRecPhi) );
                fRecC1.SetXYZ( fRecV.X()-fRecR1*sin(fRecPhi), fRecV.Y()+fRecR1*cos(fRecPhi), 0.);
                fRecC2.SetXYZ( fRecV.X()+fRecR2*sin(fRecPhi), fRecV.Y()-fRecR2*cos(fRecPhi), 0.);
                fRecV = fRecV + fPV;
                fRecC1 = fRecC1 + fPV;
                fRecC2 = fRecC2 + fPV;
                fCutSatisfied = 1;
                if ( (R1>0)&&(R2>0)&&(D1>0)&&(D2>0)&&((Y11-D1)>0)&&((Y12-D2)>0) ) fSignSatisfied = 1;
                else fSignSatisfied = 0;
            }
            else if (i==fLoop-1) {
                fSolved = 1;  
                fCutSatisfied = 1; 
                fLoop = i+1;
                fPhiE = tmpPhiE; 
                fRE = tmpRE;    
                fRadiusE = tmpRadiusE;              
                fRecR1 = R1; 
                fRecR2 = R2; 
                fRecR = ( (Y11 - D1) + (Y12 - D2) ) / 2.0; 
                fRecPhi = NextPhi;
                fRecV.SetX( fRecR * cos(fRecPhi) );
                fRecV.SetY( fRecR * sin(fRecPhi) );
                fRecC1.SetXYZ( fRecV.X()-fRecR1*sin(fRecPhi), fRecV.Y()+fRecR1*cos(fRecPhi), 0.);
                fRecC2.SetXYZ( fRecV.X()+fRecR2*sin(fRecPhi), fRecV.Y()-fRecR2*cos(fRecPhi), 0.);
                fRecV = fRecV + fPV;
                fRecC1 = fRecC1 + fPV;
                fRecC2 = fRecC2 + fPV;
                fCutSatisfied = 0;
                if ( (R1>0)&&(R2>0)&&(D1>0)&&(D2>0)&&((Y11-D1)>0)&&((Y12-D2)>0) ) fSignSatisfied = 1;
                else fSignSatisfied = 0;
            }
            // B. Cut threshold NOT satisfied - prepare for next loop
            prevR1 = R1;   prevR2 = R2; 

        }
    }
}


void Conv4HitsReco2::Dump()
{
    std::cout << std::endl<< "================================================" << std::endl;
    std::cout << "  Nothing happend here.";
    if (fSolved==1) std::cout << "Solved.";
    if (fCutSatisfied==1) std::cout << "Cut good.";
    if (fSignSatisfied==1) std::cout << "Sign good.";
    
}

math::XYZVector Conv4HitsReco2::GetPlusCenter(double &plusR)
{
    plusR = fRecR1;
    return fRecC1;
    
}
math::XYZVector Conv4HitsReco2::GetMinusCenter(double &minusR)
{
    minusR = fRecR2;
    return fRecC2;
}