---

TrackerRecon.cxx

Go to the documentation of this file.

//-------------------------------------------------------------------
//  $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/src/TrackerRecon.cxx,v 1.13 2001/09/05 15:52:43 lsrea Exp $   */
//
//     TrackerRecon:
//
//          Steers the Silicon-Tracker Reconstruction    
//
//                    Bill Atwood
//                    B. Atwood, JA Hernando, Santa Cruz, 02/05/99
//
//-------------------------------------------------------------------


#include "reconstruction/TrackerRecon.h"
#include "reconstruction/GlastRecon.h"
/*
#include "instrument/Glast.h"
#include "instrument/Tower.h"
#include "instrument/SiTracker.h"
*/
#include "DetectorData.h"


#include "reconstruction/ReconVisitor.h"
#include "reconstruction/SiClusters.h"
#include "reconstruction/GlastFit.h"
#include "reconstruction/CalRecon.h"
#include "reconstruction/MCRecon.h"

// analysis function includes
#include "reconstruction/ActiveDistance.h"
#include "reconstruction/EnergyCorrection.h"
#include "reconstruction/ExtraHits.h"
#include "reconstruction/TowerBoundaries.h"
#include "reconstruction/TrackerVeto.h"

#include "gui/DisplayRep.h"

#include <algorithm>

#define MAX_NO_FITS 15
#define Z_DISPLAY_LEN 50
#define TEST_MIN 5.0
#define IST_CHI2SELECT 10.0
#define CHI2SELECT 250.
#define RES1_CUT 10.
#define RES2_CUT 5

//#define MAX_NO_FITS 10
//#define Z_DISPLAY_LEN 50
//#define TEST_MIN 8.0
//#define IST_CHI2SELECT 6.0
//#define CHI2SELECT 50.
//#define RES1_CUT 10.
//#define RES2_CUT 5.

//unused: static float testMax;
//unused: static float firstPlaneWt;
//unused: static int numFits;
/*
static bool kalPrivateTup = false;
static int kalglast_event=0;
static int kaltrack_event=0;
static int kalplane_event=0; */

// Update TrackerRecon
static bool CONTROL_writeOut = false;
static int CONTROL_debug = 2;
static double CONTROL_sigmaCut = 8.; // 10.
// following duplicates what is in GlastFit.cxx
//unused: static bool                CONTROL_newReconstruction  = true;
static double              CONTROL_minEnergy          = 0.03;
// static GFdata::particle CONTROL_particle = GFdata::ELECTRON;
// unused: static double CUTOFF_ENERGY=0.02f;
// unused: static double FENE = 1.;   // 0.75-1;
static double XENE = 0.6; // 0.66666-0.75;
static unsigned int MAX_CANDIDATES = 4;
// unused: static bool VETO = true;
// static double MIN_QUALITY = 3.;
// unused:static int MAX_ATTEMPTS = 16;

static int CONTROL_maxGammaTries = 1;
static int CONTROL_maxParticleTries = 15;
static int CONTROL_maxConsecutiveGaps = 2;

static double CONTROL_minQ = -100.;

static double CONTROL_FEne = 2./3.;  // defoult 2/3.
static double CONTROL_FEneGamma = 1.;
static double CONTROL_FEneParticle = 1./3.;
//unused: static double CUT_minQualityParticle = 0.;
//unused: static double CONTROL_maxMoreTracks = 3;
//unused: static double CONTROL_eneMoreTracks = 0.5;  // fraction of the energy given to the more tracks

inline static double sqr(double x){return x*x;}

using namespace std;

// declare static class data members here
TrackerRecon::Parameters TrackerRecon::s_parms; 

// define a pointer to the s_parms data member in the local namespace
namespace {
    const TrackerRecon::Parameters*   parms = &TrackerRecon::params();
    
};

// implement the TrackerRecon::Parameters class

// initialize: sets up parameters from the various GLAST instrument
//  and medium classes.
void    TrackerRecon::Parameters::initialize ()
{
    // initialize kluge parameter structure
    num_convs    = DetectorData::SiTracker::numTrays()-1;
    conv_gap     = DetectorData::SiTracker::traySpacing();
    si_thickness = DetectorData::SiTracker::siThickness();
    conv_rad_len = DetectorData::SiTracker::convRadLen(1);
    tower_width  = DetectorData::Tower::mod_width + 2.*DetectorData::Tower::wall_thickness + DetectorData::Glast::wall_gap;

}

void TrackerRecon::initialize()
{
    
    // Silicon Hits
    i_Conv_Layer_Hits     = data.addElement("Cnv_Lyr_Hits");
    i_max_controller_hits = data.addElement("max_controller_hits");
    i_First_Conv_Layer    = data.addElement("fst_Cnv_Lyr");
    i_nConv_Layers_Hit    = data.addElement("nCnv_Lyrs_Hit");
    i_Conv_Layer_Hits_Gamma = data.addElement("Cnv_Hits_Gamma");
    i_nConv_Layers_Gamma    = data.addElement("nCnv_Lyrs_Gamma");
    i_First_Conv_Layer_Gamma = data.addElement("fst_Cnv_Lyr_Gamma");
    
    // Tracker Reconstruction
    i_Kal_Energy         = data.addElement("KalEne_Gamma");
    i_Kal_Energy_X       = data.addElement("KalEne_X");
    i_Kal_Energy_Y       = data.addElement("KalEne_Y");
    
    // Event Reconstruction
    i_No_XTrks            = data.addElement("No_X_Trks");
    i_No_YTrks            = data.addElement("No_Y_Trks");
    i_No_Trks             = data.addElement("No_Tracks");
    
    // Gamma quality
    i_fit_type            = data.addElement("Fit_Type");
    i_fit_topo            = data.addElement("Fit_Topo");
    i_Chisq_1st           = data.addElement("Chisq_1st");
    i_Chisq               = data.addElement("Chisq");
    i_No_Hits             = data.addElement("No_Trk_Hits");
    i_No_SHits            = data.addElement("Sum_Hits");
    i_No_Gaps             = data.addElement("Fit_Gaps");
    i_No_Gaps_1st         = data.addElement("First_Fit_Gaps");
    i_No_Noise            = data.addElement("Noise_Hits");
    i_No_Noise_1st        = data.addElement("First_Noise_Hits");
    i_qual_X              = data.addElement("X_Quality_Parm");
    i_qual_Y              = data.addElement("Y_Quality_Parm");
    i_qual                = data.addElement("Quality_Parm");
    
    // Pair reconstruction
    i_Fit_XNhits          = data.addElement("Fit_XNhits");
    i_Fit_YNhits          = data.addElement("Fit_YNhits");
    i_Fit_XChisq          = data.addElement("Fit_XChisq");
    i_Fit_YChisq          = data.addElement("Fit_YChisq");
    i_Fit_XChisq_1st      = data.addElement("Fit_XChisq_1st");
    i_Fit_YChisq_1st      = data.addElement("Fit_YChisq_1st");
    i_Fit_XKalThetaMS     = data.addElement("Fit_XKalTheta");
    i_Fit_YKalThetaMS     = data.addElement("Fit_YKalTheta");
    i_Fit_XKalEne         = data.addElement("Fit_XKalEne");
    i_Fit_YKalEne         = data.addElement("Fit_YKalEne");
    i_Fit_x0              = data.addElement("Fit_x0");
    i_Fit_y0              = data.addElement("Fit_y0");
    i_Fit_z0              = data.addElement("Fit_z0");
    i_Fit_xdir            = data.addElement("Fit_x_dir");
    i_Fit_ydir            = data.addElement("Fit_y_dir");
    i_Fit_zdir            = data.addElement("Fit_z_dir");
    
    i_Pair_XNhits         = data.addElement("Pair_XNhits");
    i_Pair_YNhits         = data.addElement("Pair_YNhits");
    i_Pair_XChisq         = data.addElement("Pair_XChisq");
    i_Pair_YChisq         = data.addElement("Pair_YChisq");
    i_Pair_XChisq_1st     = data.addElement("Pair_XChisq_1st");
    i_Pair_YChisq_1st     = data.addElement("Pair_YChisq_1st");
    i_Pair_XKalThetaMS    = data.addElement("Pair_XKalTheta");
    i_Pair_YKalThetaMS    = data.addElement("Pair_YKalTheta");
    i_Pair_XKalEne        = data.addElement("Pair_XKalEne");
    i_Pair_YKalEne        = data.addElement("Pair_YKalEne");
    i_Pair_x0             = data.addElement("Pair_x0");
    i_Pair_y0             = data.addElement("Pair_y0");
    i_Pair_z0             = data.addElement("Pair_z0");
    i_Pair_xdir           = data.addElement("Pair_x_dir");
    i_Pair_ydir           = data.addElement("Pair_y_dir");
    i_Pair_zdir           = data.addElement("Pair_z_dir");
    
    i_e_frac              = data.addElement("Energy_Split");
    i_errSlopeX           = data.addElement("KalErr_Xslope");
    i_errSlopeY           = data.addElement("KalErr_Yslope");
    i_weightXSlope        = data.addElement("weight_Xslope");
    i_weightYSlope        = data.addElement("weight_Yslope");
    i_xeneXSlope          = data.addElement("xene_Xslope");
    i_xeneYSlope          = data.addElement("yene_Yslope");      
    
    // gamma reconstruction
    i_First_XHit          = data.addElement("fst_X_Lyr");
    i_First_YHit          = data.addElement("fst_Y_Lyr");
    i_Diff_1st_Hit        = data.addElement("Diff_1st_XY_Lyr");
    i_Gamma_x0            = data.addElement("Gamma_x0");
    i_Gamma_y0            = data.addElement("Gamma_y0");
    i_Gamma_z0            = data.addElement("Gamma_z0");
    i_Gamma_xdir          = data.addElement("Gamma_x_dir");
    i_Gamma_ydir          = data.addElement("Gamma_y_dir");
    i_Gamma_zdir          = data.addElement("Gamma_z_dir");
    
    // analysis
    i_Gamma_DLT           = data.addElement("t_Angle");
    i_Zdiff_XY            = data.addElement("Zdiff_XY");
    i_Zdiff_plane         = data.addElement("Zdiff_plane");
    
    // HMA (GSFC) silicon hit clusters (within a tray) ntuple elements 6/3/1999
    i_FitSide_nXused      = data.addElement("FitSide_nXused");
    i_FitSide_nYused      = data.addElement("FitSide_nYused");
    
    // extra analysis JAH 7/9/99
    i_fit_kink = data.addElement("Fit_Kink");
    i_pair_kink = data.addElement("Pair_Kink");
    i_fit_kinkN = data.addElement("Fit_KinkN");
    i_pair_kinkN = data.addElement("Pair_KinkN");
    
    i_hits_previous = data.addElement("Hits_Previous");
    i_hits_first = data.addElement("Hits_First");
    i_hits_next = data.addElement("Hits_Next");
    
    _mclsData = new SiClusters(parms->num_convs);
    
    // now that the 'data' member is initialized,
    // set the data for all of the analysis variable objects.
    // note that the tracker recon object still 
    // creates these values as well - however only
    // to keep the tuple items in the same position in
    // the tuple 
    if(m_activeDist!=0) { // skip if not doing this junk
    m_activeDist->setData (&data);
    m_Ecorr->setData (&data);
    m_extraHits->setData (&data);
    m_twrBounds->setData (&data);
    m_tkrVeto->setData (&data);
    }
    _mGFgamma = 0;

}

TrackerRecon::TrackerRecon()
: Recon(), cal(0), 
m_activeDist(0), 
m_Ecorr(0), 
m_extraHits(0), 
m_twrBounds(0),
m_tkrVeto(0) 
{
    setTitle("TrackerRecon");
    
    // initialize kluge parameter structure
    s_parms.initialize();
    initialize();
}
 
//######################################################
TrackerRecon::TrackerRecon (GlastRecon* r)
//###################################################
: Recon(), cal(r->getCal())
,  m_activeDist(new ActiveDistance(r)), 
m_Ecorr(new EnergyCorrection(r)), 
m_extraHits(new ExtraHits(r)), 
m_twrBounds(new TowerBoundaries(r)),
m_tkrVeto(new TrackerVeto(r)) 
{
    setTitle("TrackerRecon");
    
    // initialize kluge parameter structure
    s_parms.initialize();
    initialize();

}

//##############################################
TrackerRecon::~TrackerRecon ()
//##############################################
{
    delete _mclsData;
    
    delete m_Ecorr;
    delete m_activeDist;
    delete m_extraHits;
    delete m_twrBounds;
    delete m_tkrVeto;
}

//##########################################################
void TrackerRecon::accept (ReconVisitor& d) {d.visit(this);}
//###########################################################

//###################################
void TrackerRecon::clear ()
//###################################
{
    Recon::clear();
    data.clear();
    
    if (_mGFgamma != 0) delete _mGFgamma;
    _mGFgamma = 0;
    
    xtrackList.clear();
    ytrackList.clear();
    
    _mclsData->clear();
}

//#########################################
void TrackerRecon::draw (gui::DisplayRep& v)
//#########################################
{
    int nxtracks = xtrackList.size();
    int nytracks = ytrackList.size();
    
    // X-Tracks
    KalHit::TYPE type = KalHit::SMOOTH;
    int itrack=0;
    for (; itrack<nxtracks; itrack++){
        KalTrack xtrack = xtrackList[itrack];
        if( itrack ==0) v.setColor("blue");
        else v.setColor("navy");
        xtrack.drawTrack(v, SiData::X, type);
        xtrack.drawChiSq(v, SiData::X, type);
        
        if (CONTROL_writeOut) xtrack.writeOut(CONTROL_writeOut);
    }
    
    // Y-Tracks
    for (itrack=0; itrack<nytracks; itrack++){
        KalTrack ytrack = ytrackList[itrack];
        if( itrack ==0) v.setColor("blue");
        else v.setColor("navy");
        ytrack.drawTrack(v, SiData::Y, type);
        ytrack.drawChiSq(v, SiData::Y, type);
        
        if (CONTROL_writeOut) ytrack.writeOut(CONTROL_writeOut);
    }

    //Gamma
    if(_mGFgamma) {
        v.setColor("yellow");
        v.moveTo(_mGFgamma->vertex());
        v.lineTo(_mGFgamma->vertex() - 300.*_mGFgamma->direction());
    }

    v.setColor("black");
}
//---------------------------------------
//   Voila, monsieur le reconstructeur
//---------------------------------------
//###########################################################
void TrackerRecon::reconstruct(const SiData* stripData)
//###########################################################
{
    
    if(state() == DONE) return;
    
    RCreconstruct(stripData);
    
    
}
//--------------------------------------------------------
//     Segmented Original TrackerRecon - Flow private fuctions
//--------------------------------------------------------
//########################################################
void TrackerRecon::iniTrackerRecon()
//########################################################
{
    // Data Members for Associations
    x1 = Point(0.,0.,0.);
    t1 = Vector(0.,0.,0.);
    t0 = Vector(0.,0.,0.);
    x0 = Point(0.,0.,0.);
    
    m_totalHits = 0;
    m_firstConvLayer = -1;
    m_numLayersHit = 0;
    
    m_CsIEnergy = 0.;
    m_CsICorrEnergy = 0.;
    m_aveCalHit = Point(0.,0.,0.);
    
    
    m_firstLayer = -1;
    m_fitType = -1;
    m_xFactor = 0.;
    m_yFactor = 0.; 
    
    // clear all of the analysis objects here (if they were made)
    if( m_activeDist !=0 ) {
        m_activeDist->clear();
        m_Ecorr->clear();
        m_extraHits->clear();
        m_twrBounds->clear();
        m_tkrVeto->clear();
    }
}
//########################################################
void TrackerRecon::iniClusters(const SiData* stripData)
//########################################################
{
    
    m_totalHits = 0;
    m_numLayersHit = 0;
    m_firstConvLayer = -1;
    
    int i;
    for(i=0; i<parms->num_convs; i++) {
        int hitCount = stripData->nHits(SiData::X, i) + stripData->nHits(SiData::Y, i);
        if(hitCount > 0) {
            if(m_firstConvLayer < 0) {
                m_firstConvLayer = i;
            }
            m_totalHits+= hitCount;
            m_numLayersHit++;
        }
    }
    // Now make clusters for tracking
    _mclsData->loadFrom(*stripData);

    data.addData(i_Conv_Layer_Hits, m_totalHits);
//THB    data.addData(i_max_controller_hits, stripData->maxControllerHits());
    data.addData(i_First_Conv_Layer, m_firstConvLayer);
    data.addData(i_nConv_Layers_Hit, m_numLayersHit);
    
    int iGammaLayer = 100;
    int nGammaHits = 0;
    int nGammaLayers = 0;
    int indexhit = -1;
    int jndexhit = -1;
    Point Dummy;
    for (i = m_firstConvLayer; i < parms->num_convs ; i++) {
        bool oklayer = false;
        // initialize X coordinate
        Dummy = _mclsData->nextHit(SiData::X, -1, &indexhit);
        do {
            Dummy = _mclsData->nextHit(SiData::X, i, &indexhit);
            if (indexhit > 0) {
                if (_mclsData->clusterNoise(SiData::X, indexhit) == 0) {
                    nGammaHits++;
                    oklayer = true;
                }
            }
        } while (indexhit >0);
        // initialize Y coordinate
        Dummy = _mclsData->nextHit(SiData::Y,-1, &jndexhit);
        do {
            Dummy = _mclsData->nextHit(SiData::Y, i, &jndexhit);
            if (jndexhit > 0) {
                if (_mclsData->clusterNoise(SiData::Y,jndexhit) == 0) {
                    nGammaHits++;
                    oklayer = true;
                }
            }
        } while (jndexhit >0);
        if (oklayer) {
            nGammaLayers++;
            if (i < iGammaLayer) iGammaLayer = i;
        }
    }
    
    data.addData(i_Conv_Layer_Hits_Gamma, nGammaHits);
    data.addData(i_First_Conv_Layer_Gamma, iGammaLayer);
    data.addData(i_nConv_Layers_Gamma, nGammaLayers);
}

//########################################################
void TrackerRecon::iniCalorimeterInput()
//########################################################
{    
    m_CsIEnergy = 0.;
    m_CsICorrEnergy = 0.;
    m_aveCalHit = Point(0.,0.,0.);
    m_CsIEnergy = (double) cal->data.getData("CsI_Energy_Sum");
    m_CsICorrEnergy = m_CsIEnergy;

    // JAH 7/29/99
    // store this information for the moment
    // if there is no RC event in the tracker, at least the
    // corrected Energy is equal to the measured by the Calorimeter
    m_Ecorr->setCsIEnergy (m_CsICorrEnergy);
    
    // Note: Hit correction coef. depends on the radiator thickness and strip pitch
    //       and noise... present values are for 3.5% Pb, 200 um, and 5x10**-5 occup.
    m_CsICorrEnergy = m_CsIEnergy + 0.0006 * (m_totalHits - 55);
    if(m_CsIEnergy < .005) m_CsICorrEnergy = .025 + 0.0008 * (m_totalHits - 55);
    

    m_aveCalHit =Point(cal->data.getData("CsI_X"),
        cal->data.getData("CsI_Y"),
        cal->data.getData("CsI_Z"));
    
    // Dirty Only for e-quiet
    // m_CsIEnergy = (4./3.)*0.1;
}

//-----------------------------------------------------------
//     Update TrackerRecon - for Pair Fit
//-----------------------------------------------------------

//###########################################################
void TrackerRecon::RCreconstruct(const SiData* stripData)
//###########################################################
{
    
    //unused: int MAX_PAIRTRIES = 1;
    
    if(state() == DONE) return;
    
    iniTrackerRecon();    
    iniClusters(stripData);
    if(!(cal == 0)){
        iniCalorimeterInput();
    }
    
    GFtutor::_DATA = _mclsData;
    
    bool okRC = false;
    
    // Gamma reconstruction
    int ntries = 0;
    bool doit = true;
    while (ntries < CONTROL_maxGammaTries && doit) {
        ntries ++;
        bool seedOK = seedCandidates(GFdata::PAIR);
        bool gammaOK = false;
        if (seedOK) gammaOK = candidates(GFdata::GAMMA, GFdata::PAIR);
        if (gammaOK) selectGamma();
        else doit = false;
        okRC = okRC || gammaOK;
    } 
    
    // Other tracks reconstruction
    ntries = 0;
    doit = true; 
    while (ntries < CONTROL_maxParticleTries && doit) {
        GFtutor::CUT_veto = false;
        ntries ++;
        bool seedOK = seedCandidates(GFdata::TRACK);
        bool particleOK = false;
        if (seedOK) particleOK = candidates(GFdata::PARTICLE, GFdata::TRACK);
        if (particleOK) selectParticle();
        else doit = false;
        okRC = okRC || particleOK;
        GFtutor::CUT_veto = true;
    };
    
    if (okRC) loadTuple(stripData);
    // if (kalPrivateTup) KalPrivateTup();
}

//###########################################################
bool TrackerRecon::seedCandidates(enum GFdata::type type)
//###########################################################
{
    bool ok = false;
    
    m_Xcandidates.clear();
    m_Ycandidates.clear();
    
    bool okX = findSeedCandidates(m_Xcandidates, type, SiData::X);
    bool okY = findSeedCandidates(m_Ycandidates, type, SiData::Y);
    
    if (okX || okY) ok = true; 
    return ok;
}

//###########################################################
bool TrackerRecon::candidates(enum GFdata::type type, enum GFdata::type seedType)
//###########################################################
{
    bool OK = false;
    
    m_candidates.clear();
    
    unsigned int ix =0;
    unsigned int iy =0;
    for (; ix < m_Xcandidates.size(); ix++) {
        for (; iy < m_Ycandidates.size(); iy++) {
            GFdata Xcandidate = m_Xcandidates[ix];
            GFdata Ycandidate = m_Ycandidates[iy];
            bool ok = false;
            ok = findCandidates(m_candidates, Xcandidate, Ycandidate, type);
            OK = OK || ok;
        } // Y candidates
    } // X candidates  
    
    if (OK) return OK;
    
    // Force a PairFit with no veto
    bool save_veto = GFtutor::CUT_veto;
    std::vector<GFdata> candidates;
    candidates.clear();
    for (ix = 0 ; ix < m_Xcandidates.size(); ix++) {
        GFdata Xcandidate = m_Xcandidates[ix];
        candidates.clear();
        GFtutor::CUT_veto = false;

        findSeedCandidates(candidates, seedType, SiData::Y,Xcandidate.firstLayer(),Xcandidate.tower());
        for (iy =0 ; iy < candidates.size(); iy++) {
            GFtutor::CUT_veto = save_veto;
            GFdata Ycandidate = candidates[iy];
            bool ok = findCandidates(m_candidates,Xcandidate,Ycandidate, type);
            OK = OK || ok;
        } // Y candidates
    } // X candidates
    
    candidates.clear();
    for (iy = 0 ; iy < m_Ycandidates.size(); iy++) {
        GFdata Ycandidate = m_Ycandidates[iy];
        candidates.clear();
        GFtutor::CUT_veto = false;

        findSeedCandidates(candidates, seedType, SiData::X,Ycandidate.firstLayer(),Ycandidate.tower());
        for (ix =0 ; ix < candidates.size(); ix++) {
            GFtutor::CUT_veto = save_veto;
            GFdata Xcandidate = candidates[ix];
            bool ok = findCandidates(m_candidates,Xcandidate,Ycandidate, type);
            OK = OK || ok;
        } // Y candidates
    } // X candidates
    GFtutor::CUT_veto = save_veto;
    return OK;
}

//###########################################################
bool TrackerRecon::findCandidates(std::vector<GFdata>& candidates,
                                  const GFdata& Xcandidate, const GFdata& Ycandidate,
                                  enum GFdata::type type)
                                  //###########################################################
{
    bool ok = false;
    int naccepted = 0;
    
    int iniLayer = (Xcandidate.firstLayer() < Ycandidate.firstLayer()?
        Xcandidate.firstLayer() : Ycandidate.firstLayer());
    int lastLayer = (Xcandidate.firstLayer()< Ycandidate.firstLayer()?
        Ycandidate.firstLayer() : Xcandidate.firstLayer());
    
    if (lastLayer - iniLayer >= CONTROL_maxConsecutiveGaps) return ok;
    
    Point ver =GFdata::doVertex(Xcandidate.ray(),Ycandidate.ray());
    Vector dir = GFdata::doDirection(Xcandidate.direction(),Ycandidate.direction());
    Ray testRay(ver, dir);
    //unused:   int nhits = Xcandidate.nhits()+Ycandidate.nhits();
    //m_CsICorrEnergy = m_CsIEnergy + corrEnergy(nhits,0.);
    
    for (int ilayer = iniLayer ; ilayer <= lastLayer; ilayer++) {

        GFdata candidateGFdata = GFconstructor(type, ilayer, testRay);
        if (candidateGFdata.Q() > CONTROL_minQ) {
            ok = true;
            naccepted++;
            incorporate(candidates, candidateGFdata);
        }
    }
    
    return ok;
}

//###########################################################
GFdata TrackerRecon::GFconstructor(enum GFdata::type type, int ilayer, 
                                   const Ray testRay, SiData::Axis axis)
                                   //###########################################################
{
    GFdata data;

    if (type == GFdata::PARTICLE) {
        GFparticle* _par = new GFparticle(CONTROL_sigmaCut,
            CONTROL_FEneParticle*m_CsICorrEnergy, ilayer, testRay);
        if (!_par->empty() && _par->accept()) data = _par->getGFdata();
        delete _par;
    } else if (type == GFdata::GAMMA) {
        GFgamma* _gamma = new GFgamma(CONTROL_FEne, CONTROL_sigmaCut,
            m_CsICorrEnergy, ilayer, testRay);
        if (!_gamma->empty() && _gamma->accept()) {
            data = _gamma->getGFdata();
        }
        delete _gamma;
    } else if (type == GFdata::TRACK) {
        GFtrack* _track = new GFtrack(axis, CONTROL_sigmaCut, 
            CONTROL_FEneParticle*m_CsICorrEnergy, ilayer, testRay);
        if (!_track->empty() && _track->accept()) data = _track->getGFdata();
        delete _track;
    } else if (type == GFdata::PAIR) {
        GFpair* _pair = new GFpair(CONTROL_FEne, axis, CONTROL_sigmaCut,
            CONTROL_FEneGamma*m_CsICorrEnergy, ilayer, testRay);
        if (!_pair->empty()) {
            if (_pair->accept()) {
                data = _pair->getBest()->getGFdata();
            }
        }
        delete _pair;
    } 
    
    return data;
}
//###########################################################
void TrackerRecon::selectGamma()
//###########################################################
{
    
    if (_mGFgamma != 0) delete _mGFgamma;
    _mGFgamma = 0;
    
    GFdata gammaCandidate = m_candidates[0];
    
    int iniLayer = gammaCandidate.firstLayer();
    Ray testRay = gammaCandidate.ray();
    
    //unused:   int nhits = gammaCandidate.nhits();
    // m_CsICorrEnergy = m_CsIEnergy + corrEnergy(nhits,0.);
    
    _mGFgamma = new GFgamma(CONTROL_FEne, CONTROL_sigmaCut,
        m_CsICorrEnergy, iniLayer, testRay);
    
    if (CONTROL_writeOut && !_mGFgamma->empty()) _mGFgamma->writeOut(CONTROL_debug);
    
    loadTrackerRecon();
    
    _mGFgamma->flagAllHits();
    
}
//###########################################################
void TrackerRecon::selectParticle()
//###########################################################
{
    
    GFdata candidate = m_candidates[0];
    
    int iniLayer = candidate.firstLayer();
    Ray testRay = candidate.ray();
    
    //unused:   int nhits = candidate.nhits();
    // m_CsICorrEnergy = m_CsIEnergy + corrEnergy(nhits,0.);
    
    GFparticle* _GFparticle = new GFparticle(CONTROL_sigmaCut, 
        CONTROL_FEneParticle*m_CsICorrEnergy, iniLayer, testRay);
    
    if (_GFparticle->Q() > CONTROL_minQ) {
        if (!_GFparticle->getXGFtrack()->empty()) {
            xtrackList.push_back(
                _GFparticle->getXGFtrack()->getKalTrack());
        }
        if (!_GFparticle->getYGFtrack()->empty()) {
            ytrackList.push_back(
                _GFparticle->getYGFtrack()->getKalTrack());
        }
    }
    _GFparticle->flagAllHits();
    
    delete _GFparticle;
}

//###########################################################
void TrackerRecon::loadTrackerRecon()
//###########################################################
{
    m_firstLayer = _mGFgamma->firstLayer();
    m_fitType = 0; // need some work
    m_xFactor =XENE;
    m_yFactor = XENE;
    
    x1 = GFdata::doVertex(_mGFgamma->getBest(SiData::X)->ray(),
        _mGFgamma->getBest(SiData::Y)->ray());
    t1 = GFdata::doDirection(_mGFgamma->getBest(SiData::X)->direction(),
        _mGFgamma->getBest(SiData::Y)->direction());
    
    if (!_mGFgamma->getPair(SiData::X)->empty() && !_mGFgamma->getPair(SiData::Y)->empty()) {
        x2 = GFdata::doVertex(_mGFgamma->getPair(SiData::X)->ray(),
            _mGFgamma->getPair(SiData::Y)->ray());
        t2 = GFdata::doDirection(_mGFgamma->getPair(SiData::X)->direction(),
            _mGFgamma->getPair(SiData::Y)->direction());
    } else {
        x2 = x1;
        t2 = t1;
    }
    
    x0 = _mGFgamma->vertex();
    t0 = _mGFgamma->direction();
    
    
    if (!_mGFgamma->getBest(SiData::X)->empty()) xtrackList.push_back(
        _mGFgamma->getBest(SiData::X)->getKalTrack());
    if (!_mGFgamma->getPair(SiData::X)->empty()) xtrackList.push_back(
        _mGFgamma->getPair(SiData::X)->getKalTrack());
    
    if (!_mGFgamma->getBest(SiData::Y)->empty()) ytrackList.push_back(
        _mGFgamma->getBest(SiData::Y)->getKalTrack());
    if (!_mGFgamma->getPair(SiData::Y)->empty()) ytrackList.push_back(
        _mGFgamma->getPair(SiData::Y)->getKalTrack());
    
}

//---------------------------------------------------------
//     Update TrackerRecon - for GlastFit - Second level
//---------------------------------------------------------

//###########################################################
bool TrackerRecon::findSeedCandidates(std::vector<GFdata>& candidates, GFdata::type type, SiData::Axis axis)
//###########################################################
{
    bool OK = false;   
    int indexhit = -1;
    Point PIni = _mclsData->nextHit(axis, -1, &indexhit); 
    
    int ilayer=m_firstConvLayer;

    for ( ; ilayer < DetectorData::SiTracker::numTrays()-3; ilayer++){
        bool ok = findSeedCandidates(candidates, type, axis, ilayer);
        OK = OK || ok;
    }
    return OK;
}

//###########################################################
bool TrackerRecon::findSeedCandidates(std::vector<GFdata>& candidates, GFdata::type type, 
                                      SiData::Axis axis, int ilayer, int itower)
//###########################################################
{
    //unused:   int nconstructed = 0;
    int naccepted = 0;
    bool ok = false;
    
    if (ilayer >= DetectorData::SiTracker::numTrays()-3) return ok;
    
    // Initialize tray
    int indexhit = -1;
    Point PIni = _mclsData->nextHit(axis, -1, &indexhit); 
    
    Point PRef;
    if (m_CsIEnergy > CONTROL_minEnergy) PRef = m_aveCalHit;
    
    int ihit=0;
    int nhitsInLayer = _mclsData->nHits(axis,ilayer);

    //PRef.printOn(std::cout);
    //std::cout << "\n";

    for ( ; ihit<nhitsInLayer; ihit++) {
        // Get Initial Point
        PIni = _mclsData->nextHit(axis,ilayer,&indexhit);

        if (indexhit<0)  {
            break;
        }        
        if (SiClusters::tower(indexhit) != itower && itower != 0) {
            continue;
        }
        
        // Get Reference Point
        //  if (PRef.mag() == 0. || m_CsIEnergy <= CONTROL_minEnergy)
        PRef = createPRef(axis, ilayer, PIni);

        // Get GFpair within layer, PIni and PRef - Inputs
        if (PIni.mag() == 0. || PRef.mag() == 0.)  {
            continue;
        }

        
        Vector VDir(PRef.x()-PIni.x(),PRef.y()-PIni.y(),PRef.z()-PIni.z());
        Ray testRay = Ray(PIni, VDir.unit());
        
        GFdata candidateGFdata = GFconstructor(type, ilayer, testRay, axis);
        
        if (candidateGFdata.Q() > CONTROL_minQ) {
            ok = true;
            naccepted++;
            incorporate(candidates, candidateGFdata);
        }
    }
    
    if (naccepted > 0) ok = true;
    return ok;
    
}
//###########################################################
void TrackerRecon::incorporate(std::vector<GFdata>& pDatalist, const GFdata pData)
//###########################################################
{
    bool ienter = false;
    for (unsigned int i =0; i < pDatalist.size(); i++) {
        if (pData.Q() >= pDatalist[i].Q()) {
            pDatalist.insert(&pDatalist[i],pData);
            ienter = true;
        }
        if (ienter) break;
    }
    if (!ienter) pDatalist.push_back(pData);
    if (pDatalist.size()>MAX_CANDIDATES) pDatalist.pop_back();
    
}
//########################################################
Point TrackerRecon::createPRef(SiData::Axis axis,int ilayer, const Point& Pini)
//########################################################
{
    double weight = CONTROL_minEnergy/(3.*m_CsICorrEnergy); // weight;
    
    Point PCal = m_aveCalHit;

    if (m_CsIEnergy < CONTROL_minEnergy) weight = 1.;
    if (PCal.mag() == 0) weight = 1.;
    Point PTrk = _mclsData -> meanHitInside(ilayer+1, s_parms.tower_width, Pini);

    if (PTrk.mag() == 0.) {
        PTrk = _mclsData -> meanHitInside(ilayer+2, s_parms.tower_width, Pini);  
    }
    if (PTrk.mag() == 0.) weight = 0.;
    double x = (1.-weight)*PCal.x()+weight*PTrk.x();  
    double y = (1.-weight)*PCal.y()+weight*PTrk.y();  
    double z = (1.-weight)*PCal.z()+weight*PTrk.z();

    Point PRef(x,y,z);

    if (PRef.mag() != 0.) {
        if (axis == SiData::X ) y = Pini.y();
        else x = Pini.x();
    } else {
        x = Pini.x();
        y = Pini.y();
        z = Pini.z()-DetectorData::SiTracker::traySpacing();
    }
    PRef = Point(x,y,z);

    return PRef;
} 
//--------------------------------------------------------
//     Kalman Private Ntuple - Debug Version
//--------------------------------------------------------
//########################################################
void TrackerRecon::loadTuple(const SiData* _stripData)
//########################################################
{           
    // Load results into data array
    //-----------------------------
    
    data.addData(i_No_XTrks, xtrackList.size());
    data.addData(i_No_YTrks, ytrackList.size());
    data.addData(i_No_Trks,  xtrackList.size() +ytrackList.size());
    
    if (_mGFgamma == 0) return;
    if (_mGFgamma->empty()) return;
    
    GFtrack* _Xbest = _mGFgamma->getXpair()->getBest();
    GFtrack* _Ybest = _mGFgamma->getYpair()->getBest();
    
    data.addData(i_Kal_Energy,_mGFgamma->RCenergy());
    data.addData(i_Kal_Energy_X, _mGFgamma->getXpair()->RCenergy());
    data.addData(i_Kal_Energy_Y, _mGFgamma->getYpair()->RCenergy());
    
    // data.add(i_TRC_code, 10*xtrackList.size()+ytrackList.size()); 
    
    // Gamma quality
    double type = 0;
    if (!_mGFgamma->empty()) type = 11;
    if (!_mGFgamma->getPair(SiData::X)->empty()) type += 1;
    if (!_mGFgamma->getPair(SiData::Y)->empty()) type += 10;
    data.addData(i_fit_type, type);
    
    double topo = 0;
    if (_mGFgamma->fix()) topo = 1.;
    if (_mGFgamma->conflictPattern()) topo *= -1.;
    data.addData(i_fit_topo, topo);
    
    data.addData(i_Chisq , _Xbest->chiSquare() + _Ybest->chiSquare());
    data.addData(i_Chisq_1st , _Xbest->chiSquareSegment() + _Ybest->chiSquareSegment());
    data.addData(i_qual_X ,_Xbest->Q());
    data.addData(i_qual_Y ,_Ybest->Q());
    data.addData(i_qual ,  _Xbest->Q() + _Ybest->Q());
    data.addData(i_No_Hits, _Xbest->numDataPoints()+_Ybest->numDataPoints());
    data.addData(i_No_Gaps, _Xbest->numGaps()+_Ybest->numGaps());
    data.addData(i_No_Gaps_1st,  _Xbest->numFirstGaps()+_Ybest->numFirstGaps());
    data.addData(i_No_Noise      ,_Xbest->numNoise()+_Ybest->numNoise());
    data.addData(i_No_Noise_1st  ,_Xbest->numFirstNoise()+_Ybest->numFirstNoise());
    
    int nxhits = _Xbest->numDataPoints();
    int nyhits = _Ybest->numDataPoints();
    
    if (!_mGFgamma->getXpair()->getPair()->empty()) {
        int nhits = _mGFgamma->getXpair()->getPair()->numDataPoints();
        if (nhits < nxhits) nxhits = nhits;
    }
    if (!_mGFgamma->getYpair()->getPair()->empty()) {
        int nhits = _mGFgamma->getYpair()->getPair()->numDataPoints();
        if (nhits< nyhits) nyhits = nhits;
    }

    // Pair reconstruction
    data.addData(i_Fit_XNhits ,_Xbest->nhits());
    data.addData(i_Fit_YNhits ,_Ybest->nhits());
    data.addData(i_Fit_XChisq ,_Xbest->chiSquare());
    data.addData(i_Fit_YChisq ,_Ybest->chiSquare());
    data.addData(i_Fit_XChisq_1st ,_Xbest->chiSquareSegment());
    data.addData(i_Fit_YChisq_1st ,_Ybest->chiSquareSegment());
    //data.addData(i_Fit_XChisq_Com ,_Xbest->computeChiSqSegment(nxhits));
    //data.addData(i_Fit_YChisq_Com ,_Ybest->computeChiSqSegment(nyhits));
    data.addData(i_Fit_XKalEne ,_Xbest->RCenergy());
    data.addData(i_Fit_YKalEne ,_Ybest->RCenergy());
    data.addData(i_Fit_XKalThetaMS ,_Xbest->KalThetaMS());
    data.addData(i_Fit_YKalThetaMS ,_Ybest->KalThetaMS());
    data.addData(i_Fit_xdir,  t1.x());
    data.addData(i_Fit_ydir,  t1.y());
    data.addData(i_Fit_zdir,  t1.z());
    data.addData(i_Fit_x0,    x1.x());
    data.addData(i_Fit_y0,    x1.y());
    data.addData(i_Fit_z0,    x1.z());
    
    GFtrack* _Xpair = _mGFgamma->getXpair()->getPair();
    GFtrack* _Ypair = _mGFgamma->getYpair()->getPair();

    if (!_Xpair->empty()) {
        data.addData(i_Pair_XNhits ,_Xpair->nhits());
        data.addData(i_Pair_XChisq ,_Xpair->chiSquare());
        data.addData(i_Pair_XChisq_1st ,_Xpair->chiSquareSegment());
//        data.addData(i_Pair_XChisq_Com ,_Xpair->computeChiSqSegment(nxhits));
        data.addData(i_Pair_XKalEne ,_Xpair->RCenergy());
        data.addData(i_Pair_XKalThetaMS ,_Xpair->KalThetaMS());
    }
    if (!_Ypair->empty()) {
        data.addData(i_Pair_YNhits ,_Ypair->nhits());
        data.addData(i_Pair_YChisq ,_Ypair->chiSquare());
        data.addData(i_Pair_YChisq_1st ,_Ypair->chiSquareSegment());
//        data.addData(i_Pair_YChisq_Com ,_Ypair->computeChiSqSegment(nyhits));
        data.addData(i_Pair_YKalEne ,_Ypair->RCenergy());        
        data.addData(i_Pair_YKalThetaMS ,_Ypair->KalThetaMS());
    }
    data.addData(i_Pair_x0,  x2.x());
    data.addData(i_Pair_y0,  x2.y());
    data.addData(i_Pair_z0,  x2.z());
    data.addData(i_Pair_xdir,  t2.x());
    data.addData(i_Pair_ydir,  t2.y());
    data.addData(i_Pair_zdir,  t2.z());
    
    data.addData(i_Gamma_xdir,t0.x());
    data.addData(i_Gamma_ydir,t0.y());
    data.addData(i_Gamma_zdir,t0.z());  
    data.addData(i_Gamma_x0,  x0.x());
    data.addData(i_Gamma_y0,  x0.y());
    data.addData(i_Gamma_z0,  x0.z());
    
    data.addData(i_e_frac, (m_xFactor+m_yFactor)/2.);
    data.addData(i_errSlopeX,     _mGFgamma->getXpair()->errorSlope());             
    data.addData(i_errSlopeY,     _mGFgamma->getYpair()->errorSlope());            
    data.addData(i_xeneXSlope,    _mGFgamma->getXpair()->RCenergy());   
    data.addData(i_xeneYSlope,    _mGFgamma->getYpair()->RCenergy());      
    data.addData(i_weightXSlope,  _mGFgamma->getXpair()->weightSlope());        
    data.addData(i_weightYSlope,  _mGFgamma->getYpair()->weightSlope());          
    
    float sinDLT = sqrt(std::max(0.0,(1. - sqr(t0*t1))));
    data.addData(i_Gamma_DLT, sinDLT);
    
    //  First Used Layer and previous Hits...
    int x_layer = _Xbest->firstLayer();
    data.addData(i_First_XHit, x_layer);
    int y_layer = _Ybest->firstLayer();
    int diffXY = x_layer - y_layer;
    data.addData(i_Diff_1st_Hit, diffXY);
    
    // Tracker Recon Analysis Variables
    //----------------------------------
    
    // compute the tracker veto parameters

    // ******* Take Note *******
    if(cal == 0){
        return;
    }
    m_tkrVeto->compute();
    
    // compute the tower/track boundary parameters
    m_twrBounds->compute();
    
    // compute the "Active_Dist" parameter
    m_activeDist->compute();
    
    double Zdiff_XY= _mGFgamma->getXpair()->vertex().z()-_mGFgamma->getYpair()->vertex().z();
    data.addData(i_Zdiff_XY, Zdiff_XY);
    Point firstHit = _mGFgamma->getFirstHit();
    double Zdiff_plane = _mGFgamma->vertex().z()-firstHit.z();
    data.addData(i_Zdiff_plane,Zdiff_plane);
    
    // compute the extra hits parameters
    m_extraHits->compute();
    
    // compute the CsI corrected energy
    m_Ecorr->compute();
    
    // make sure the data member for the CsICorrEnergy is set...
    m_CsICorrEnergy = m_Ecorr->result()->energyCorr();
    
    // kink in the first segment
    double fitkink,pairkink,fitkinkN,pairkinkN;
    fitkink = TAna_kink(pairkink, fitkinkN, pairkinkN);
    data.addData(i_fit_kink,fitkink);
    data.addData(i_pair_kink,pairkink);
    data.addData(i_fit_kinkN,fitkinkN);
    data.addData(i_pair_kinkN,pairkinkN);
    
    // hits in the previos and next plane
    double previousHits, firstHits, nextHits;
    firstHits = TAna_densityHits(previousHits, nextHits);
    data.addData(i_hits_previous,previousHits);
    data.addData(i_hits_first,firstHits);
    data.addData(i_hits_next,nextHits);
}

//########################################################
double TrackerRecon::TAna_trackerVeto(double& xsigma, double& ysigma)
//########################################################
{
    // X view
    GFtrack* _Xbest = _mGFgamma->getXpair()->getBest();
    GFtrack* _Ybest = _mGFgamma->getYpair()->getBest();
    
    int indexhit;
    bool Xveto = _Xbest->veto(indexhit,xsigma);
    bool Yveto = _Ybest->veto(indexhit,ysigma);
    
    double icode = 0;
    if (Xveto) icode = 1;
    if (Yveto) icode += 2;
    return icode;
}
//########################################################
double TrackerRecon::TAna_kink(double& pairKink, double& fitKinkN, double& pairKinkN)
//########################################################
{
    double fitKink = 0.;
    pairKink = 0.;
    
    GFtrack* xtrk = _mGFgamma->getBest(SiData::X);
    GFtrack* ytrk = _mGFgamma->getBest(SiData::Y);
    
    if (xtrk->empty() || ytrk->empty()) return fitKink;
    
    double xkink = xtrk->kink(0);
    double ykink = ytrk->kink(0);
    fitKink = (fabs(ykink) > fabs(xkink)? ykink:xkink);
    
    double xkinkN = xtrk->kinkNorma(0);
    double ykinkN = ytrk->kinkNorma(0);
    fitKinkN = (fabs(ykinkN) > fabs(xkinkN)? ykinkN:xkinkN);
    
    // Pair Track
    xtrk = _mGFgamma->getPair(SiData::X);
    ytrk = _mGFgamma->getPair(SiData::Y);
    
    if (xtrk->empty() || ytrk->empty()) return fitKink;
    
    xkink = xtrk->kink(0);
    ykink = ytrk->kink(0);
    pairKink = (fabs(ykink) > fabs(xkink)? ykink:xkink);
    
    xkinkN = xtrk->kinkNorma(0);
    ykinkN = ytrk->kinkNorma(0);
    pairKinkN = (fabs(ykinkN) > fabs(xkinkN)? ykinkN:xkinkN);
    
    return fitKink;
}
//########################################################
double TrackerRecon::TAna_densityHits(double& previousHits, double& nextHits)
//########################################################
{
    double firstHits = 0.;
    previousHits = 0.;
    nextHits = 0.;
    
    //  Find number of hits around first hit: 5 sigma 
    
    // get the geometrical Factors
    int firstLayer = _mGFgamma->firstLayer();
    Point firstHit = _mGFgamma->vertex();
    //unused:   int diffXY = _mGFgamma->getXpair()->firstLayer()-_mGFgamma->getYpair()->firstLayer();
    //unused:   float dz    =  _mGFgamma->getBest(SiData::X)->vertex().z() - 
    //unused:           _mGFgamma->getBest(SiData::Y)->vertex().z();
    
    float hitRadFact = fabs(t0.z()) + sqrt(1. - t0.z()*t0.z())/fabs(t0.z());
    if(hitRadFact > 3.) hitRadFact = 3.;
    int ipln = parms->num_convs - firstLayer; 
    float radLenEff = (DetectorData::SiTracker::convRadLen(ipln) + 2.*parms->si_thickness/9.36 +.003)/fabs(t1.z());
    //unused:    float thetaCone = .014/(m_CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
    //unused:    float minSprd  = 5.* thetaCone * parms->conv_gap; // 5 sigma cut
    //unused:    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;
    }
    //unused:    float dltX = dltSprd*xFact;
    //unused:    float dltY = dltSprd*yFact;
    
    float deltaS = parms->conv_gap/fabs(t0.z());
    
    // previous
    int iplane = -1;
    int i;
    for (i=std::max(firstLayer-1,0);i<=firstLayer+1;i++) {
        ipln = parms->num_convs - i;
        radLenEff +=  (DetectorData::SiTracker::convRadLen(ipln) + 2.*parms->si_thickness/9.36 +.003)/fabs(t1.z());
        float thetaMS = .014/(m_CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
        float s_total = (i-firstLayer) * deltaS;
        float xSprd = thetaMS * deltaS * xFact * 2.021; // = 3.5 sigma/sqrt(3)
        float ySprd = thetaMS * deltaS * yFact * 2.021;
        Point trkHit((Point)(s_total*t0 + x0));
        if (xSprd > 50) xSprd = 50.;
        if (ySprd > 50) ySprd = 50.;
        float nearHits = -1.;
        if (i>=0) nearHits = _mclsData->numberOfHitsNear(i, xSprd, ySprd, trkHit);
        iplane++;
        if (iplane == 0) previousHits = nearHits;
        else if (iplane == 1) firstHits = nearHits;
        else if (iplane == 2) nextHits = nearHits;
    }
    
    return firstHits;
}
//########################################################
double TrackerRecon::TAna_towerBoundary(double& xborder, double& yborder)
//########################################################
{
    
/*
double x_ist=0;
double y_ist=0;
double z_ist=0;
z_ist = (_Xbest->getHit(0).z() > _Ybest->getHit(0).z()) ?
_Xbest->getHit(0).z() : _Ybest->getHit(0).z();

  float dz    =  _Xbest->getHit(0).z() - _Ybest->getHit(0).z();
  if(dz > 0.) {
  x_ist = _Xbest->getHit(0).x();
  y_ist = _Ybest->getHit(0).x() + t0.y()*dz/t0.z();
  }
  else {
  x_ist = _Xbest->getHit(0).x() - t0.x()*dz/t0.z();
  y_ist = _Ybest->getHit(0).x();
  }
  Point firstHit = Point(x_ist, y_ist, z_ist);
  Point vertex = x0 + t0*(z_ist/t0.z());
    m_firstLayer = (x_layer < y_layer) ? x_layer : y_layer; */
    // unused:  double border = 0.;
    
    Point vertex = _mGFgamma->getFirstHit();
    float x_ist = vertex.x();
    float y_ist = vertex.y();
    
    float num_2 = DetectorData::Glast::xNum/2.; 
    float x_twr_bnd = fmod((x_ist + s_parms.tower_width*num_2), s_parms.tower_width);
    x_twr_bnd = x_twr_bnd - s_parms.tower_width*.5;
    float y_twr_bnd = fmod((y_ist + s_parms.tower_width*num_2), s_parms.tower_width);
    y_twr_bnd = y_twr_bnd - s_parms.tower_width*.5;
    
    x_twr_bnd = fabs(x_twr_bnd);
    y_twr_bnd = fabs(y_twr_bnd);
    float twr_bnd   = (x_twr_bnd > y_twr_bnd) ? x_twr_bnd : y_twr_bnd;
    
    xborder = x_twr_bnd;
    yborder = y_twr_bnd;
    
    return twr_bnd;
}

//########################################################
double TrackerRecon::TAna_activeDistance()
//########################################################
{  
    
    float activeDist = -20.;
#if 0 //THB: This depends on simulation environment, functionality needs to be restored
    Point firstHit = _mGFgamma->getFirstHit();
    Ray prjRay(firstHit, -t0);
    Point x_prj;
    
    GlastMed *glast = GlastMed::instance();
    
    int istLayer = _mGFgamma->firstLayer();
    
    float firstProp = 0.;
    if (istLayer > 0) { //go to next gap up + espsilon
        firstProp = SiTracker::traySpacing()/fabs(t0.z());
        x_prj = prjRay.position(firstProp);
    }
    else {
        x_prj = firstHit;
        firstProp = 0.;
    }
    
    const Medium *med0 = glast->inside(x_prj);
    prjRay = Ray(x_prj, -t0);
    const Medium *detMed = 0;
    if(med0){
        firstProp = med0->distanceToLeave(prjRay, detMed, 100.);
        activeDist = -2. - firstProp;
    }
    
    Detector *det =  0;
    if(detMed) {
        med0 = detMed;
        x_prj = prjRay.position(firstProp);
        prjRay = Ray(x_prj, -t0);
        firstProp = med0->distanceToLeave(prjRay, detMed, 100.);
        activeDist = -4. - firstProp;
        if(detMed) det = detMed->detector();
        if(det) {
            if(!strcmp(det->nameOf(), "SiDetector")) {
                x_prj = prjRay.position(firstProp);
                activeDist = static_cast<const SiDetectorMed*>(det)->data().insideActiveArea(x_prj);
            }
        }
    }
#endif    
    return activeDist;
}

//########################################################
double TrackerRecon::TAna_extraHits(double& norma, double& istHitCount, double& OutHits,
                                    double& ShwHits1, double& ShwHits2)
                                    //########################################################
{  
    
    double SumHits = 0.;
    norma = 1.;
    OutHits = 0.;
    ShwHits1 = 0.;
    ShwHits2 = 0.;
    
    //  Find number of hits around first hit: 5 sigma out
    int firstLayer = _mGFgamma->firstLayer();
    Point firstHit = _mGFgamma->vertex();
    int diffXY = _mGFgamma->getXpair()->firstLayer()-_mGFgamma->getYpair()->firstLayer();
    float dz    =  _mGFgamma->getBest(SiData::X)->vertex().z() - 
        _mGFgamma->getBest(SiData::Y)->vertex().z();
    
    float hitRadFact = fabs(t0.z()) + sqrt(1. - t0.z()*t0.z())/fabs(t0.z());
    if(hitRadFact > 3.) hitRadFact = 3.;
    int ipln = parms->num_convs - firstLayer; 
    float radLenEff = (DetectorData::SiTracker::convRadLen(ipln) + 2.*parms->si_thickness/9.36 +.003)/fabs(t1.z());
    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;
    }
    if(istHitCount == 1 && 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 == 1) 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(m_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 +=  (DetectorData::SiTracker::convRadLen(ipln) + 2.*parms->si_thickness/9.36 +.003)/fabs(t1.z());
        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;
    }
    
    SumHits = sumHits;
    norma = lastLayer - firstLayer;
    OutHits = outHits;
    ShwHits1 = shwHits1;
    ShwHits2 = shwHits2;
    return sumHits;
}

//########################################################
double TrackerRecon::TAna_corrEnergy(double ShwHits1, double ShwHits2, double twr_bnd)
//########################################################
{ 
    double eneCorr = 0;
    GFtrack* _Xbest = _mGFgamma->getXpair()->getBest();
    GFtrack* _Ybest = _mGFgamma->getYpair()->getBest();
    float hit_energy = 0.0003*ShwHits1  + .0013*ShwHits2;  //Energy in Hits
    float num_Layers_Crossed = _Xbest->numDataPoints()+ _Xbest->numGaps();
    if(_Xbest->numDataPoints() < _Ybest->numDataPoints()) 
        num_Layers_Crossed = _Ybest->numDataPoints()+ _Ybest->numGaps();
    float x_ing_energy = .0016*num_Layers_Crossed; // Take out observed layer correlation
    float edg_corr   =  (twr_bnd < 10.) ? 1 : .68 + (.069  - .0036* twr_bnd )*twr_bnd; // Edge of Tower correction 
    float trk_energy = hit_energy + x_ing_energy;  // Total Tracker energy  
    trk_energy /= (fabs(t0.z()) > .2) ? fabs(t0.z()) : .2; // Note: t0.z < 0
    eneCorr = (m_CsIEnergy + trk_energy)/edg_corr;
    return eneCorr;    
}
/*
//########################################################
void TrackerRecon::KalPrivateTup()
//########################################################
{

  int iwrite=0; 
  // Planes Tuple - Information of the XBest planes
  if (xtrackList.size() >0) {
                GlastFit* _xtrack=_mbestXFit;
                std::string title("kalplanes.tup");
                for (int iplane=0; iplane<_xtrack->kplanelist.size(); iplane++){
                kalplane_event==0 && iplane ==0? iwrite=1: iwrite=0; 
                KalTup planetup(kalplane_event, iplane, *_xtrack);
                planetup.output(iwrite,title);
                }
                kalplane_event++;
                }
                
                  // Tracks tuple
                  if (xtrackList.size() >0 && ytrackList.size() >0 ) {
                  std::string title_e("kalevent.tup");
                  kaltrack_event==0? iwrite=1: iwrite=0;
                  KalTup eventtup(kaltrack_event, xtrackList[0], ytrackList[0]);
                  eventtup.output(iwrite, title_e);
                  kaltrack_event++;
                  } 
} */
//--------------------------------------------------------
//     Correction of the Energy
//--------------------------------------------------------
//########################################################
double TrackerRecon::corrEnergy(int level, enum configuration, int nShowerHits, double slope)
//########################################################
{   
    // B. Atwood, A. Djannati, JA Hernando.
    
    // level for the Ecorrection
    //          1 - total number of hits in the tracker
    //          2 - number of hits in the gamma shower
    //          3 - and correction for the geometry leakeages
    
    double deltaE = 0.;
    
    if (level == 1) {
        // this correction is Hardwire for the baseline 4x4 and 5x10^-5 occupancy
        //  JA hernando, D. Melton 6/30/99
        // deltaE = 0.0004*(nShowerHits - 55);
        // previous code by Bill Atwood (5x5 but not relevant changes are expected
        deltaE = 0.0006 * (m_totalHits - 55);
        deltaE = 0.025 + 0.0008 * (m_totalHits - 55);
    } else if (level == 2) {
    } else {
    }
    
    return deltaE;
}


GFgamma*    TrackerRecon::gammaFit ()
{
    return _mGFgamma;
}

---