---

ExtraHits.cxx

Go to the documentation of this file.

// $Id: ExtraHits.cxx,v 1.5 2000/12/11 16:38:45 burnett Exp $
// 
//  Original author: Sawyer Gillespie
//                   hgillesp@u.washington.edu
//

#include "reconstruction/ExtraHits.h"

#include "reconstruction/GlastFit.h"
#include "reconstruction/GlastRecon.h"
#include "reconstruction/LbldData.h"
#include "reconstruction/SiClusters.h"
#include "reconstruction/TrackerRecon.h"

//TODO: unhardwire this. It was SiTracker::convRadLen(ipln)
static double conv_rad_len = 0.25000;


// implementation of the ExtraHits class

ExtraHits::ExtraHits (GlastRecon* r)
: ReconAnalysis (r)
{
}

ExtraHits::~ExtraHits ()
{
}

void    ExtraHits::setUpData ()
{
    i_1st_count = data()->addElement("fst_Hit_Count");
    i_srp_hit_rat = data()->addElement("Surplus_Hit_Ratio");
    i_out_hit_rat = data()->addElement("OutSide_Hit_Ratio");
    i_ShowerHits1 = data()->addElement("showerHits1");
    i_ShowerHits2 = data()->addElement("showerHits2");
    i_No_SHits = data()->addElement("Sum_Hits");
}

void    ExtraHits::compute ()
{
    // get the TrackerRecon object which has reconstructed the
    // data within the silicon tracker
    TrackerRecon*   tr = recon()->getTracker();
    
    // declare local variables

    // get the gamma-ray fit from TrackerRecon
    GFgamma*    gamma = tr->gammaFit ();

    // get the instrument parameters from TrackerRecon
    const TrackerRecon::Parameters* parms = &tr->params ();

    // get parameters corresponding to the fit gamma-ray trajectory
    float t0x = tr->data.getData("Gamma_x_dir");
    float t0y = tr->data.getData("Gamma_y_dir");
    float t0z = tr->data.getData("Gamma_z_dir");
    Vector  t0 (t0x, t0y, t0z);
    Point   x0 = gamma->vertex();

    // get parameters corresponding to the track fit?
    float t1z = tr->data.getData("Fit_z_dir");

    // get the corrected CsI energy
    double   m_CsICorrEnergy = tr->getCsICorrEnergy();

    // get the silicon-strip cluster data
    SiClusters*   _mclsData = tr->getClusters();

    // native code - this code has been brought in from TrackerRecon
    // -------------------------------------------------------------
    double norma = 1.;
    double istHitCount = 0.;
    
    //  Find number of hits around first hit: 5 sigma out
    int firstLayer = gamma->firstLayer();
    Point firstHit = gamma->vertex();
    
    int diffXY = gamma->getXpair()->firstLayer()-gamma->getYpair()->firstLayer();
    float dz    =  gamma->getBest(SiData::X)->vertex().z() - 
        gamma->getBest(SiData::Y)->vertex().z();
    
    float hitRadFact = fabs(t0z) + sqrt(1. - t0z*t0z)/fabs(t0z);
    if(hitRadFact > 3.) hitRadFact = 3.;
    int ipln = parms->num_convs - firstLayer; 
    float radLenEff = (conv_rad_len + 2.*parms->si_thickness/9.36 +.003)/fabs(t1z);
    float thetaCone = .014/(m_CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
    float minSprd  = 5.* thetaCone * parms->conv_gap; // 5 sigma cut
    float dltSprd  = minSprd;
    float norm = sqrt(t0.x()*t0.x() + t0.y()*t0.y());
    float xFact = 1., yFact = 1.;
    if(norm > 10000.*FLT_MIN) {
        xFact = 1. + (hitRadFact - 1.)*fabs(t0.x())/norm;
        yFact = 1. + (hitRadFact - 1.)*fabs(t0.y())/norm;
    }
    float dltX = dltSprd*xFact;
    float dltY = dltSprd*yFact;
    
    // Compute the First_hit_count variable
    istHitCount = _mclsData->numberOfHitsNear(firstLayer, dltX, dltY, firstHit);
    // separate the special case of 3 into 3 or 2.5
    if(istHitCount == 3 && diffXY == 0) { // Check if extra hits on the far side of the first gap
        SiData::Axis far_axis = SiData::X;
        double dR = dltX;
        if (dz>0) {
            far_axis = SiData::Y;
            dR = dltY;
        }
        int ihit_second = _mclsData->numberOfHitsNear(far_axis, firstLayer, dltX, firstHit);
        if (ihit_second == 2) istHitCount -=.5;
    }
    
    //  Find the number of hits around the rest of the gamma trajectory
    float deltaS = parms->conv_gap/fabs(t0.z());
    double disp = deltaS;
    int lastLayer = (int)(4+2.*log(m_CsICorrEnergy/.01) + firstLayer);
    if(lastLayer > parms->num_convs) lastLayer = parms->num_convs;
    
    float sumHits = _mclsData->numberOfHitsNear(firstLayer, .5*dltX, .5*dltY, firstHit);
    float shwHits1 = sumHits;
    float shwHits2 = 0; 
    if(firstLayer > 11){
        shwHits1 = 0; 
        shwHits2 = sumHits;
    }
    
    if(sumHits < 2.) sumHits = 2.;
    if(lastLayer - firstLayer < 5 && lastLayer == parms->num_convs) sumHits +=.5;
    float outHits = 0;
    
    xFact *= sqrt(t0.z() *t0.z() + t0.x()*t0.x());
    yFact *= sqrt(t0.z() *t0.z() + t0.y()*t0.y());
    
    int i;
    for(i=firstLayer+1; i<parms->num_convs; i++) {
        ipln = parms->num_convs - i;
        radLenEff +=  (conv_rad_len + 2.*parms->si_thickness/9.36 +.003)/fabs(t1z);
        float thetaMS = .014/(m_CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
        float s_total = (i-firstLayer) * deltaS;
        float xSprd = thetaMS * s_total * xFact * 2.021; // = 3.5 sigma/sqrt(3)
        float ySprd = thetaMS * s_total * yFact * 2.021;
        Point trkHit((Point)(disp*t0 + x0));
        if (xSprd > 50) xSprd = 50.;
        if (ySprd > 50) ySprd = 50.;
        float nearHits = _mclsData->numberOfHitsNear(i, xSprd, ySprd, trkHit);
        if(i< 12) shwHits1 +=  nearHits;
        else      shwHits2 +=  nearHits; 
        if( i < lastLayer) {
            sumHits +=  nearHits;
            outHits +=  _mclsData->numberOfHitsNear(i, 50.,   50.,   trkHit) - nearHits;
        }
        disp += deltaS;
    }
    
    norma = (lastLayer - firstLayer);
    //-----------------------------------------------------------
    // end native code
        
    // set the state of this object to reflect the results of the
    // computation
    ReconAnalysis::result (new Result(sumHits, norma, istHitCount, outHits, shwHits1, shwHits2));

    // add the data to the list of archival data elements.
    tr->data.addData (i_1st_count, istHitCount);
    tr->data.addData (i_srp_hit_rat, sumHits/norma);
    tr->data.addData (i_out_hit_rat, outHits/norma);
    tr->data.addData (i_No_SHits, sumHits);
    tr->data.addData (i_ShowerHits1, shwHits1);
    tr->data.addData (i_ShowerHits2, shwHits2);
}

---