---

TowerArray.cxx

Go to the documentation of this file.

// $Id: TowerArray.cxx,v 1.4 2001/01/23 03:41:32 burnett Exp $
//
// This file is part of Gismo 2
//
//    Implementation of class TowerArray and TowerElement member functions
//


#include "gismo/TowerArray.h"

#include "gismo/MCParticle.h"
#include "geometry/Box.h"
#include "geometry/CoordTransform.h"
#include "geomrep/BoxRep.h"
#include "Random.h"
#include <algorithm>
float
TowerArray::s_display_scale=10.;

void
TowerArray::score(MCParticle* track)
{
  float deltaE = track->energyLoss();
  if( deltaE<=0 || track->mass()==0 && track->status()==MCParticle::ALIVE)
          return;
  Point pos(track->position());
  pos.transform(globalToLocal());
  long  i = calculateIndex(pos);
  if(i < 0) return;

  scoreTower(i, deltaE);
  return;
}

float
TowerArray::scoreTower(long ix, float deltaE)
{
    bool found=false;
    float etot = deltaE;

    for( iterator tower = begin(); tower != end(); ++tower) {
        TowerElement& s = *tower;
        if( ix > s.index()) continue;
        if( ix == s.index() ) {
           etot = (s += deltaE);
        }
        else {
        // add before the next, to keep in order of the index
           insert(tower, TowerElement( ix, deltaE) );
        }
        found = true;
        break;
    }

    //  If the tower isn't in the list yet, and no others have higher ids,
    //  add it to the end of the list.
   if( !found ) {
        push_back(TowerElement( ix, deltaE) );
   }
   return etot;

}

long
TowerArray::calculateIndex(const Point& point) const
{
  int ix = (int)((point.x()/xWidth+0.5)*xSegmentation);
  int iy = (int)((point.y()/yWidth+0.5)*ySegmentation);

  ix = std::min(xSegmentation-1,std::max(ix,0));
  iy = std::min(ySegmentation-1,std::max(iy,0));
  /* Assume now that always is in, or close enough to boundary
  if( ix<0 ||ix>=xSegmentation ) return -1;
  if( iy<0 || iy>=ySegmentation ) return -1;
  */
  return ix*ySegmentation + iy;


}
void
TowerArray::calculateBin(long ind, float& xmin, float& ymin) const
{
  // reverse calculateIndex
  float fy = ind%ySegmentation;
  float fx = ind/ySegmentation;
  xmin = (fx/xSegmentation-0.5)*xWidth + xWidth/xSegmentation/2.;
  ymin = (fy/ySegmentation-0.5)*yWidth + yWidth/ySegmentation/2.;

}
int TowerArray::hits()const
{ return (int)size();}

float TowerArray::energy(unsigned i)const
{   return towers(i)->energy();
}

Point
TowerArray::xtalPos(unsigned i)const
{
   float xmin, ymin;
   calculateBin(towers(i)->index(), xmin, ymin);

   Point center(xmin//+.5*xWidth/xSegmentation
                 ,ymin//+0.5*yWidth/ySegmentation
                 ,0);
   center.transform(localToGlobal());
   return center;
}

long TowerArray::xtalId(unsigned i)const
{  return towers(i)->index(); }

void
TowerArray::clear()
{
    TowerElementList::erase(begin(), end());
}

void
TowerArray::generateResponse()
{
    for(iterator it = begin(); it != end(); ++it){
        TowerElement& theTower = *it;

        theTower.simEnergy = theTower.energy();
        if(energyResolution > 0)
            theTower.simEnergy +=  Random::gauss() * energyResolution *
                                   sqrt(theTower.energy());
    }
}

Point *
TowerArray::meanEnergy(float Ethres, float *Etotal, int *numTowers)
{
    float xCenter, yCenter;
    float eX = 0.;
    float eY = 0.;
    *numTowers = 0;

    *Etotal = 0.;
    for(iterator theTower = begin(); theTower<end(); ++theTower){

        if(theTower->simEnergy >= Ethres) {
           *Etotal += theTower->simEnergy;
           calculateBin(theTower->index(), xCenter, yCenter);
           eX += xCenter * theTower->simEnergy;
           eY += yCenter * theTower->simEnergy;
           *numTowers += 1;
        }
      }
    if(*numTowers > 0) {
       eX = eX/(*Etotal);
       eY = eY/(*Etotal);
    }
    Point *retHit = new Point(eX, eY, zPlane);
    retHit->transform(localToGlobal());
    return retHit;
}

void
TowerArray::draw(gui::DisplayRep& view)const
{
    if(size() == 0) return;

    const CoordTransform T= localToGlobal();

    for(const_iterator t = begin(); t<end(); ++t){
        const TowerElement& theTower = *t;
        float xmin, ymin;
        float height = theTower.energy()* displayScale();
        if( height <= 0 ) continue;
        calculateBin(theTower.index(), xmin, ymin);
        Box box(xWidth/xSegmentation,
                  yWidth/ySegmentation,
                  height );

        box.move(Vector(xmin,ymin,zPlane+height/2));
        box.transform(T);

        view.append(BoxRep(box));
    }
}

void TowerArray::readData(std::istream& is)
{
    clear();
    int count;
    is >> count;
    while (count --)
    {  unsigned id; int ie;
       is >> id >> ie;
       push_back(TowerElement(id, ie/1e6));
    }

}
void TowerArray::writeData(std::ostream& os)
{
    os << size() << '\n';
    for(iterator theTower = begin(); theTower<end(); ++theTower){
       os << theTower->index() << '\t'<< int(1e6*theTower->energy()) << '\n';

    }
}
void TowerArray::printOn(std::ostream& cout)const
{
  float total=0, xmin,ymin;
    for(const_iterator theTower = begin(); theTower<end(); ++theTower){

    calculateBin(theTower->index(),xmin,ymin);
    cout << theTower->index() << ' '
         << xmin << ' ' << ymin << ' '
         << theTower->energy() << '\n';
    total += theTower->energy();
  }
  cout << "total: " << total << '\n';
}


// writeParameters -- just enough to function when reconstructed
void TowerArray::writeParameters(std::ostream& out) const{
   out << xSegmentation <<'\t'
       << ySegmentation <<'\t'
       << xWidth        <<'\t'
       << yWidth        <<'\t'
       << zPlane        <<'\t';
}


TowerArray::~TowerArray()
{
    clear();
}
// =========================================================

TowerElement::~TowerElement()
{
}

---