---

CalRecon.cxx

Go to the documentation of this file.

// $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/src/CalRecon.cxx,v 1.11 2001/01/22 15:53:55 burnett Exp $


#include "reconstruction/CalRecon.h"
/*
#include "instrument/Calorimeter.h"
#include "instrument/Glast.h"
#include "instrument/Tower.h"
*/

#include "reconstruction/ReconVisitor.h"
#include "reconstruction/CalModules.h"

#include "data/GlastData.h"
#include "data/CsIData.h"
#include "gui/DisplayRep.h"

#include <algorithm> // for sort
#include <cmath>

using std::vector;


CalRecon::CalRecon ()
: Recon()
{
    getParameters();
    setTitle("Glast CalRecon: ver 2.0");

    layerLabel.reserve(num_layers);

    // Initialize the list; using position from Left/Right recombination
    for(unsigned int i = 0; i < num_layers; i++) {
        LRpos.push_back(new vector<Point>);
        layerLabel[i] = new char[13];
        sprintf(layerLabel[i], "%s%d%c", "CsI_eLayer", (i+1), '\0');
    }

    i_CsI_eLayer.reserve(num_layers);

    // Initialize Labelled Data Array...
    i_CsI_Energy  = data.addElement("CsI_Energy_Sum");
    for(unsigned int j = 0; j < num_layers; j++){ 
      i_CsI_eLayer[j] = data.addElement(layerLabel[j]); 
    }

    i_CsI_Etop      = data.addElement("CsI_Etop_Sum");
    i_CsI_Ebottom   = data.addElement("CsI_Ebot_Sum");
    i_CsI_Esides    = data.addElement("CsI_Esid_Sum");
    i_CsI_S_or_B    = data.addElement("CsI_S_or_B");
    i_CsI_Xtals   = data.addElement("CsI_No_Xtals");
    i_CsI_Xtals_trunc = data.addElement("CsI_No_Xtals_Trunc");
    i_CsI_Xtals_keep = data.addElement("CsI_No_Xtals_keep");
    i_CsI_Rsq_mean = data.addElement("CsI_Rsq_mean");
    i_CsI_D3Mean = data.addElement("CsI_D3Mean");
    i_CsI_E_moment = data.addElement("CsI_E_moment");
    i_CsI_E_disp = data.addElement("CsI_E_disp");
    i_CsI_mean_X  = data.addElement("CsI_X");
    i_CsI_mean_Y  = data.addElement("CsI_Y");
    i_CsI_mean_Z  = data.addElement("CsI_Z");
    i_CsI_moment1 = data.addElement("CsI_moment1");
    i_CsI_moment2 = data.addElement("CsI_moment2");
    i_CsI_moment3 = data.addElement("CsI_moment3");
    i_CsI_dir_X   = data.addElement("CsI_Xdir");
    i_CsI_dir_Y   = data.addElement("CsI_Ydir");
    i_CsI_dir_Z   = data.addElement("CsI_Zdir");
    i_CsI_chisq   = data.addElement("CsI_Chisq");
    i_CsI_mips    = data.addElement("CsI_No_MIPs");
    // module with the max # of xtals with > 100 MeV apiece
    i_CsI_MaxMod = data.addElement("CsI_MaxModGT100MeV");
    // number of xtals, with > 100 MeV apiece, in any module 
    i_CsI_ModXtals = data.addElement("CsI_ModXtalsGT100MeV");
        
        i_Diode_Energy =  data.addElement("Diode_Energy");
    // Total energy in all the diodes
    
        // Make a Modules Object to analyze cal data, tower by tower
    calModules = new CalModules();

    clear();
}


CalRecon::~CalRecon ()
{
    unsigned int i;
    for(i = 0; i < LRpos.size(); i++) delete LRpos[i];
    for(i = 0; i < num_layers; i++) delete layerLabel[i];

    delete calModules;
}


void CalRecon::getParameters ()
{
#if 0 // remove external dependence
    num_layers = Calorimeter::numberOfLayers(); 
    calmat = CsIXtal::getMaterial(); 
    s_xtal_gap = Diode::getmaxHeight();
    s_rad_len = Calorimeter::radLen(); 
    mod_width = Tower::mod_width; 
    s_wall_thickness = Tower::wall_thickness;
    s_wall_gap = Tower::wall_gap;
    s_xNum = Glast::xNum; 
    light_att = Calorimeter::light_att();
#else
    num_layers =    8;
    calmat =        "CsI";
    s_xtal_gap =    0.018;
    s_rad_len =     9.0300732267775;  
    mod_width =     38.97;
    s_wall_thickness = 0.15;
    s_wall_gap =    0.1;
    s_xNum =        4;
    light_att =     0.5;

#endif
}

void CalRecon::accept (ReconVisitor& rv)
{
    rv.visit(this);
}

void CalRecon::reconstruct (const CsIData* csiData)
{
    using namespace idents;

    if(state() == DONE) return;

    // First do module analysis...
    calModules->load( *csiData);


    const double zero = 0.;
    double E_moment = 0., E_disp = 0.;
    E_layer.resize(num_layers, 0.);

    // Determine Total energy deposited in Calorimeter
    unsigned int ilayer; // to be compatible with old/new ANSI for loop scoping
    for(ilayer = 0; ilayer < num_layers; ilayer++) {
        numlayhits.push_back(csiData->nHits(ilayer));
        ntothits += numlayhits[ilayer];
        for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
            double energy = 0.5*(csiData->Lresp(ilayer, ixtl) + csiData->Rresp(ilayer, ixtl));
            eTotal += energy;           // Total energy in CAL
            E_layer[ilayer] += energy;  // Total energy per layer
        }
        E_moment += E_layer[ilayer] * (ilayer - num_layers/2.);
    }
    
    double diode_total=0.;
    // computes total energy in all the diodes
    for(ilayer = 0; ilayer < num_layers; ilayer++) {
          for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
            std::vector<double> energy = csiData->Diodes_Energy(ilayer,ixtl);
            if( energy.size()<4 ) continue; 
            diode_total += energy[0]+energy[1]+energy[2]+energy[3];           // Total energy in diodes
          }
    }



    if(eTotal <= 0 ) { // nothing in Cal: abort reconstruction
        setState(DONE);
        return;
    }

    double onePercent = .01*eTotal;

#if 0 //THB: wire in constants taken from simulation 
    Material *Mat = Material::hatch(calmat);
    const float outzone = Mat->radiationLength() * s_rad_len / num_layers;
    const float G_hafwid = Glast::Width * 0.5 - ACDTopMed::veto_thickness 
        - (ACDSideMed::two_layer_side) * ACDTopMed::veto_thickness - ACDTopMed::veto_standoff;
    const float XYend_breadth = G_hafwid, XYend_longitd = G_hafwid - s_xtal_gap;
    const double E_nXtal = 0.006, E_D3mean = 0.010;
    const double zCalBot = Glast::zModCntr - Glast::zModDepth / 2.; //z-coord of bottom of CAL
    double e_bottom = 0., e_sides = 0.;
    float E_100MeV = 0.100f;  // 100 MeV
    const int maxModId = (Glast::xNum*Glast::yNum);
#else
    const float outzone = 2.10000f;
    const float G_hafwid = 73.9400f;
    const float XYend_breadth = G_hafwid, XYend_longitd = G_hafwid - s_xtal_gap;
    const double E_nXtal = 0.006, E_D3mean = 0.010;
    const double zCalBot = -37.186000000000; //z-coord of bottom of CAL
    double e_bottom = 0., e_sides = 0.;
    float E_100MeV = 0.100f;  // 100 MeV
    const int maxModId = 16;

#endif
    int *nXtalsGT100;
    nXtalsGT100 = new int[maxModId];
    int i;
    for ( i = 0; i < maxModId; i++) nXtalsGT100[i] = 0;
    int numXtalTrunc = 0, numXtalkeep = 0;

    E_hits.resize(ntothits, zero);

    int ihit = -1;
    for(ilayer = 0; ilayer < num_layers; ilayer++) {    // bottom to top
        for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
            Vector myP = csiData->xtalPos(ilayer, ixtl);

            double LRx, LRy, LRz;
            double LRsum = csiData->Lresp(ilayer, ixtl) + csiData->Rresp(ilayer, ixtl);
            double LRdiff = csiData->Lresp(ilayer, ixtl) - csiData->Rresp(ilayer, ixtl);

//            double f = -3.0*LRdiff/LRsum;  // range: 0.->1. => -16.25+gap -> +16.25-gap cm
            double f = - ((2-light_att)/light_att)*LRdiff/LRsum; 
                        // light attenuation is defined now in xml file
            
                        double longitudinal_pos = f * (0.5*(mod_width-2*s_xtal_gap));

            if( ilayer%2 ) { 
                LRx = myP.x() + longitudinal_pos;
                LRy = myP.y();
            } else {
                LRx = myP.x();
                LRy = myP.y() + longitudinal_pos;
            }
            LRz = myP.z();

            // filling the new list of positions
            LRpos[ilayer]->push_back(Point(LRx, LRy, LRz));

            Vector xPos(LRx, LRy, LRz);

            double energy = LRsum*0.5; 
            E_hits[++ihit] = energy;
            bCenter += xPos * energy;

            // Count number of Xtals with > 1% of total energy, and with > E_nXtal
            if(energy > onePercent) numXtalTrunc++;
            if(energy > E_nXtal) numXtalkeep++;
            // count number of xtals with energy > E_100MeV
            // keeping track for individual towers
            if(energy > E_100MeV) ++nXtalsGT100[csiData->moduleId(ilayer, ixtl)];
            
            // Sum energy in Xtals w/in outer layers, within outzone of sides or bottom.
            if(ilayer == 0) e_bottom += energy;  // sum energy deposited in bottom layer
            if ( (ilayer != num_layers-1) && (ilayer != 0) ) {
                if( ilayer%2 ) {
                    if ( ((XYend_longitd - fabs(xPos.x())) < outzone) ||
                        ((XYend_breadth - fabs(xPos.y())) < outzone) )  {
                        e_sides += energy;
                    }
                } else {
                    if ( ((XYend_breadth - fabs(xPos.x())) < outzone) ||
                        ((XYend_longitd - fabs(xPos.y())) < outzone) ) {
                        e_sides += energy;
                    }
                }
            }
        }
        float Eterm = (E_moment - E_layer[ilayer] * (ilayer - num_layers/2.));
        E_disp += Eterm * Eterm;
    }
    
    E_disp = sqrt(E_disp);
    bCenter *= 1./eTotal;

    // which module has the most xtals with energy > 100 MeV
    int maxMod = 0, nXtals = 0;
    for(int modId = 0; modId < maxModId; ++modId) { // loop through all modules
        if(nXtalsGT100[modId] > nXtals){
            nXtals = nXtalsGT100[modId];
            maxMod = modId;
        }
    }

    // Create Use_R list and initialize to 1
    const int one = 1;
    Use_R.clear();
    Use_R.resize(ntothits, one);

    int counter = 0, maxiter = 2, sid_OR_bot = 0;
    double fracHits = 0.35, D3Mean = 0.;

    if(eTotal < 0.25) fracHits = 0.50;
    while (counter < maxiter) {
        if (ntothits > 3) {  // require at least 3 xtals to perform moments analysis
                        moments_anal(counter);
                } else {
                        chisq = -1;
                        m_dir = Vector(0., 0., 0.);
                }
        if (chisq == -1) break;

        // determine if reconstructed trajectory hits side or bottom. 1 = side, 0 = bottom
        if (counter == (maxiter-1)) {
            float zRat = (zCalBot - m_cen.z())/m_dir.z();
            float xEdge = zRat*m_dir.x() + m_cen.x();
            float yEdge = zRat*m_dir.y() + m_cen.y();
            if ((abs(xEdge) > G_hafwid) || (abs(yEdge) > G_hafwid)) sid_OR_bot = 1;
        }

        // compute perpendicular distances between longest axis and Xtal hits
        // creating a list of the distances.  sort, tagging distances above
        // {fracHits percentile} for nonuse in next iteration of moment_anal.
        // compute vertical energy moment and dispersion
        vector<double>axial_R_D3;
        vector<double>axial_copy;
        vector<double>axial_R;

        int ihit = -1;
        for(ilayer = 0; ilayer < num_layers; ilayer++) {
            if (numlayhits[ilayer] != 0) {
                for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
                    Point xPos = (*LRpos[ilayer])[ixtl];
                    double temp = RLinePnt(xPos);
                    axial_R.push_back(temp);
                    axial_copy.push_back(temp);
                    if(E_hits[++ihit] > E_D3mean) axial_R_D3.push_back(temp);
                }
            }
        }

        if(ntothits != 0) {
            std::sort(axial_R.begin(), axial_R.end());
            double disc_R = axial_R[(int)(ntothits * fracHits)];
            for(unsigned int ihit = 0; ihit < axial_R.size(); ihit++) {
                if( axial_copy[ihit] > disc_R) Use_R[ihit]=0;
            }

            if((axial_R_D3.size() > 2) && (counter == 0)) {
                std::sort(axial_R_D3.begin(), axial_R_D3.end());
                int lastindx = axial_R_D3.size() - 1;
                D3Mean = ( axial_R_D3[lastindx] + axial_R_D3[lastindx-1]
                    + axial_R_D3[lastindx-2] ) / 3.;            
            }
        }
        ++counter;
    }

    float Emax = 0.;
    int numTowers = calModules->nMods(), noMIPS = 0;
    ModuleId maxTower(0,99999);
    for(int iTower = 0; iTower < numTowers; iTower++) {
        float eTower = calModules->energy(iTower);
        if(eTower > Emax ) {
            Emax = eTower;
            maxTower = calModules->moduleId(iTower);
        }
        if(calModules->responseType(iTower) == CalModules::MIP) noMIPS++;
    }

    data[i_CsI_Energy] = eTotal; // total energy deposited in CAL
    for(unsigned int jlayer = 0; jlayer < num_layers; jlayer++) { // From top to bottom
      ilayer = num_layers-1 - jlayer;
        data[i_CsI_eLayer[jlayer]] = E_layer[ilayer];  // Energy per layer
    }
 
        data[i_Diode_Energy] = diode_total;
        

    data[i_CsI_Etop] = E_layer[(num_layers-1)];
    data[i_CsI_Ebottom] = e_bottom;
    data[i_CsI_Esides] = e_sides;
    data[i_CsI_S_or_B] = sid_OR_bot;

    data[i_CsI_Xtals] = ntothits;
    data[i_CsI_Xtals_trunc] = numXtalTrunc;
    data[i_CsI_Xtals_keep] = numXtalkeep;
    data[i_CsI_Rsq_mean] = m_Rsq_mean;
    data[i_CsI_D3Mean] = D3Mean;
    data[i_CsI_E_moment] = E_moment;
    data[i_CsI_E_disp] = E_disp;

    data[i_CsI_mean_X] = m_cen.x();
    data[i_CsI_mean_Y] = m_cen.y();
    data[i_CsI_mean_Z] = m_cen.z();
    data[i_CsI_moment1] = moment1;
    data[i_CsI_moment2] = moment2;
    data[i_CsI_moment3] = moment3;

    data[i_CsI_dir_X] = m_dir.x();
    data[i_CsI_dir_Y] = m_dir.y();
    data[i_CsI_dir_Z] = m_dir.z();

    data[i_CsI_chisq] = chisq;
    data[i_CsI_mips] = noMIPS;

    data[i_CsI_MaxMod] = maxMod;
    data[i_CsI_ModXtals] = nXtalsGT100[maxMod];

    setState(DONE);

    delete [] nXtalsGT100;
    // Destroy the lists
    unsigned int j;
    for(j=0; j < LRpos.size(); j++) {
        LRpos[j]->clear();
    }
    return;

}

void CalRecon::clear ()
{
    Recon::clear();
    data.clear();
    numlayhits.clear();
    E_hits.clear();
    E_layer.clear();
    calModules->clear();
    bCenter = Point(0,0,0);
    chisq = 0;
    eTotal = 0.0;
    ntothits = 0;
    m_dir = Vector(0., 0., 0.);
}

void CalRecon::draw (gui::DisplayRep& v)
{
    static double scale=10;
    static const char* number[]={"0","1","2"};

    v.setColor("blue");
    v.markerAt(bCenter);

    if( chisq > 0 ) {// no solution, no display

        for( int i = 0; i<3; i++) {
            double d = scale*moment[1]/moment[i];
            v.moveTo(bCenter + principal_axis[i]*d);
            v.lineTo(bCenter - principal_axis[i]*d);
            v.drawText(number[i]);
        }
    }

    v.setColor("black");
}

void CalRecon::moments_anal (int counter)
{
    float Ethres = .001f;
    double xprm = 0.,  yprm = 0.,  zprm = 0., Rsq = 0.,
        Ixx = 0., Iyy = 0., Izz = 0.,
        Ixy = 0., Ixz = 0., Iyz = 0.;

    numXtals = 0;
    if (counter == 0)   m_Rsq_mean = 0.;

    // Construct elements of (symmetric) "Inertia" Tensor:
    // See Goldstein, 1965, Chapter 5 (especially, eqs. 5-6,7,22,26).
    // Analysis easy when translated to energy centroid.

    m_cen = Point(bCenter.x(), bCenter.y(), bCenter.z());

    int ihit = -1;
    unsigned int ilayer; // for compatible old/new ANSI for-loop scoping
    for(ilayer = 0; ilayer < num_layers; ilayer++) {
        for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
            ++ihit;
            if((E_hits[ihit] > Ethres) && (Use_R[ihit] > 0)) {
                numXtals++;

                Point myPos = (*LRpos[ilayer])[ixtl];
                xprm = myPos.x() - m_cen.x();
                yprm = myPos.y() - m_cen.y();
                zprm = myPos.z() - m_cen.z();

                Rsq = xprm*xprm + yprm*yprm + zprm*zprm;
                if (counter == 0) m_Rsq_mean += Rsq;
                Ixx += (Rsq - xprm*xprm) * E_hits[ihit];
                Iyy += (Rsq - yprm*yprm) * E_hits[ihit];
                Izz += (Rsq - zprm*zprm) * E_hits[ihit];
                Ixy -= xprm*yprm * E_hits[ihit];
                Ixz -= xprm*zprm * E_hits[ihit];
                Iyz -= yprm*zprm * E_hits[ihit];
            }
        }
    }

    if ((counter == 0) && (numXtals > 0)) m_Rsq_mean = m_Rsq_mean / numXtals;

    // Fit Calorimeter Track via principal moments approach.

    double intcp, slope, resid = 0.;
    double p, q, r, a, b, rad_test, m, psi;

    // Render determinant of Inertia Tensor into cubic form.

    p = - (Ixx + Iyy + Izz);
    q = Iyy*Izz + Iyy*Ixx + Izz*Ixx - (Ixy*Ixy + Iyz*Iyz + Ixz*Ixz);
    r = - Ixx*Iyy*Izz + Ixx*Iyz*Iyz + Iyy*Ixz*Ixz +
        Izz*Ixy*Ixy - 2.*Ixy*Iyz*Ixz;

    // See CRC's Standard Mathematical Tables (19th edition), pp 103-105.
    // The substitution, y = x - p/3 converts  y^3 + p*y^2 + q*y + r = 0
    // to the form  x^3 + a*x + b = 0 .  Then, if b^2/4 + a^3/27 < 0 ,
    // there will be three real roots -- guaranteed since the Inertia Tensor
    // is symmetric.  A second substitution, x = m*cos(psi) , yields the roots.

    a = (3.*q - p*p)/3.;
    b = (2.*p*p*p - 9.*p*q + 27.*r)/27.;

    rad_test = b*b/4. + a*a*a/27.;

    float dircos[3][3];

    if ((rad_test < 0.) && (Ixy != 0.) && (Ixz != 0.) && (Iyz != 0.)  ){

        // Construct the roots, which are the principal moments.

        m = 2. * sqrt(-a/3.);
        psi = acos( 3.*b/(a*m) ) / 3.;
        moment[0] = m * cos(psi) - p/3.;
        moment[1] = m * cos(psi + 2.*M_PI/3.) - p/3.;
        moment[2] = m * cos(psi + 4.*M_PI/3.) - p/3.;

        if (counter == 0) {
            moment1 = moment[0];
            moment2 = moment[1];
            moment3 = moment[2];
        }

        // Construct direction cosines; dircos for middle root is parallel to
        // longest principal axis.

        double A, B, C, D;  // needs to be double, else FPE under OSF

        for(int iroot=0; iroot<3; iroot++) {
            A = Iyz * (Ixx - moment[iroot]) - Ixy*Ixz;
            B = Ixz * (Iyy - moment[iroot]) - Ixy*Iyz;
            C = Ixy * (Izz - moment[iroot]) - Ixz*Iyz;

            D = sqrt( 1. / ( 1./(A*A) + 1./(B*B) + 1./(C*C) ) ) / C;
            dircos[0][iroot] = D * C / A;
            dircos[1][iroot] = D * C / B;
            dircos[2][iroot] = D;

            // Create vector corresponding to roots

            principal_axis[iroot] = Vector(D*C/A, D*C/B, D);
        }

        m_dir = Vector(dircos[0][1], dircos[1][1], dircos[2][1]);

        // Chisquared = sum of residuals about principal axis, through centroid

        slope = m_dir.y() / m_dir.x();
        intcp = m_cen.y()  - slope * m_cen.x();
        chisq = 0.;

        int ihit = -1;
        for(ilayer = 0; ilayer < num_layers; ilayer++) {
            for(int ixtl = 0; ixtl < numlayhits[ilayer]; ixtl++) {
                //Point xPos = csiData->xtalPos(ilayer, ixtl);
                Point xPos = (*LRpos[ilayer])[ixtl];
                resid = (xPos.y() - (intcp + slope*xPos.x()));
                chisq += resid*resid*E_hits[++ihit];
            }
        }
        chisq = ntothits>2 ? chisq/(ntothits - 2): 0;
    }

    // ill-conditioned cubic solution:

    else {
        chisq = -1.;
        m_dir = Vector(0., 0., 0.);
    }
}

double CalRecon::RLinePnt (Point xPos)
{    
    // compute distance between a line and a point
    Vector r(xPos - m_cen);
    Vector LinePt = r.cross(m_dir);
    return( (LinePt.mag() / m_dir.mag()) );
}

---