---

CalModules.cxx

Go to the documentation of this file.

// $Id: CalModules.cxx,v 1.4 2000/12/11 16:38:44 burnett Exp $
//
// Implementation for  access to CsI  data organized by Towers
//

#include "reconstruction/CalModules.h"

#include "data/CsIData.h"


CalModules::CalModules ()
{
}


CalModules::~CalModules ()
{
   clear();
}



float CalModules::energy (unsigned int n) const
{
  return calorList[n].energy;
}

int CalModules::nMods () const
{
  return (int)calorList.size();
}

Point CalModules::calPos (unsigned int n) const
{
  return  calorList[n].pos;
}

enum CalModules::ShowerType CalModules::responseType (unsigned int n) const
{
  return (ShowerType)calorList[n].type;
}

idents::ModuleId CalModules::moduleId (unsigned int n) const
{
  return calorList[n].module;
}

void CalModules::load (const CsIData& data)
{
    using namespace idents;

    static double energyThreshold = .001; // 1 Mev...

    int numXtals = data.nHits(0);  // Hardwired for 1 layer...
    ModuleId lastMod(0,999999);
    Point pos;
    int sum =0;
    float xsum = 0., xxsum = 0., ysum = 0., yysum = 0., xysum = 0., zsum= 0.; 
    //, xCenter, yCenter;
    float esum = 0, xtalEnergy, corlCoef;
    ShowerType eventType;
    for(int i=0; i<numXtals; i++) {
       if(lastMod != data.moduleId(0, i)) {
          if(sum > 0) {
             float detx   =  esum*xxsum - xsum*xsum;
             float dety   =  esum*yysum - ysum*ysum;
             float denom = detx * dety;
             float esumsq= esum*esum;
             float numer = (esum*xysum - xsum*ysum);
             if(denom/(esumsq*esumsq) > .1)
                corlCoef = numer/sqrt(denom);
             else
                corlCoef = 1;

  // Following formula drived by plotting corlCoef vs numTowers for mu's
            double minCorr = .95 - .04*(sum -5)/(1.+exp(sum-10.));
            if(minCorr > .92) minCorr = .92;
            if(fabs(corlCoef) > minCorr && sum > 4)  eventType = MIP;
            else if(sum > 2)                         eventType = QED_SHOWER;
            else                                     eventType = SPATTER;

  // Added data to list
            calorList.push_back(CalModules_Data(Point(xsum/esum, ysum/esum, zsum/esum),
                                 esum, sum, lastMod, corlCoef, eventType));
         }
         eventType = NOT_HIT;
         xsum  = 0.; xxsum = 0.;
         ysum  = 0.; yysum = 0.;
         zsum  = 0.;
         xysum = 0.;
         esum  = 0;
         sum   = 0;
         lastMod = data.moduleId(0, i);
      }
    xtalEnergy = data.energy(0, i);
    if(xtalEnergy >= energyThreshold) {
       sum += 1;
       esum += xtalEnergy;
       pos = data.xtalPos(0, i);
       xsum  += pos.x() * xtalEnergy;
       ysum  += pos.y() * xtalEnergy;
       zsum  += pos.z() * xtalEnergy;
       xxsum += pos.x() * pos.x() * xtalEnergy;
       xysum += pos.x() * pos.y() * xtalEnergy;
       yysum += pos.y() * pos.y() * xtalEnergy;
    }
  }
 // Pick up the last module..
  if(sum > 0) {
     float detx   =  esum*xxsum - xsum*xsum;
     float dety   =  esum*yysum - ysum*ysum;
     float denom = detx * dety;
     float esumsq= esum*esum;
     float numer = (esum*xysum - xsum*ysum);
     if(denom/(esumsq*esumsq) > .1)
        corlCoef = numer/sqrt(denom);
     else
        corlCoef = 1;

     double minCorr = .95 - .04*(sum -5)/(1.+exp(sum-10.));
     if(minCorr > .92) minCorr = .92;
     if(fabs(corlCoef) > minCorr && sum > 4)  eventType = MIP;
     else if(sum > 2)                         eventType = QED_SHOWER;
     else                                     eventType = SPATTER;

     calorList.push_back(CalModules_Data(Point(xsum/esum, ysum/esum, zsum/esum),
                                 esum, sum, lastMod, corlCoef, eventType));
  }
}

void CalModules::clear ()
{
   calorList.clear();
}

void CalModules::printOn (std::ostream& cout) const
{
  cout << "\nCalModules: "<< calorList.size() << " Towers \n";
  for(unsigned i=0; i < calorList.size(); i++)
  {
      cout << '\t' << energy(i) << '\t' << calPos(i) << '\n';
  }
}

// Additional Declarations

---