EmcRecTofMatch.cxx

Go to the documentation of this file.

//
//  Match Emc Shower and Tof Hit
//
//  He Miao 2006, 9, 15
//
#include "EmcRec/EmcRecTofMatch.h"
#include "EmcRec/EmcRecParameter.h"
 
EmcRecTofMatch::EmcRecTofMatch()
{}
 
EmcRecTofMatch::~EmcRecTofMatch()
{}
 
void EmcRecTofMatch::Match(RecEmcShowerMap& aShowerMap,
    RecEmcTofHitMap& tofHitMap)
{
  RecEmcShowerMap::iterator iShowerMap;
  for(iShowerMap=aShowerMap.begin();
      iShowerMap!=aShowerMap.end();
      iShowerMap++) {
 
    // calculate matching window
    double eShower = iShowerMap->second.e5x5();
    double matWindow = 0;
    if(eShower>0) {
      matWindow = 0.01277+0.004268/sqrt(eShower);
    }
    
    double phi = iShowerMap->second.phi();
    phi = phi < 0 ? phi+CLHEP::twopi : phi;
 
    int nphi1 = (int)(phi*88./CLHEP::twopi);      //tof layer1 number
    int nphi2 = (int)(phi*88./CLHEP::twopi+0.5);  //tof layer2 number
    nphi2 += 88;
 
    //find max energy in both layer of tof
    int nphi1max=nphi1;
    double emax=tofHitMap[nphi1].Energy();
    if(tofHitMap[nphi1-1].Energy() > emax) {
      emax = tofHitMap[nphi1-1].Energy();
      nphi1max = nphi1-1;
    }
    if(tofHitMap[nphi1+1].Energy() > emax) {
      emax = tofHitMap[nphi1+1].Energy();
      nphi1max = nphi1+1;
    }
    nphi1 = nphi1max;
 
    int nphi2max=nphi2;
    emax=tofHitMap[nphi2].Energy();
    if(tofHitMap[nphi2-1].Energy() > emax) {
      emax = tofHitMap[nphi2-1].Energy();
      nphi2max = nphi2-1;
    }
    if(tofHitMap[nphi2+1].Energy() > emax) {
      emax = tofHitMap[nphi2+1].Energy();
      nphi2max = nphi2+1;
    }
    nphi2 = nphi2max;
     
    double etof2x1=0;
    double etof2x3=0;
    int nphi[6] = { nphi1-1, nphi1, nphi1+1, nphi2-1, nphi2, nphi2+1 };
    
    for(int i=0;i<6;i++) {
      
      if(tofHitMap[nphi[i]].Energy()<=0) continue;
      
      HepPoint3D pos=tofHitMap[nphi[i]].Position();
      /*if(i<1) {
        //tof layer1 center
        pos = HepPoint3D(0,838.5,tofHitMap[nphi[i]].ZPosition());
      } else {
        //tof layer2 center
        pos = HepPoint3D(0,895.5,tofHitMap[nphi[i]].ZPosition());
      }*/
      
      //matching window
      if(fabs(pos.theta()-iShowerMap->second.theta())<3*matWindow) {
        etof2x3 += tofHitMap[nphi[i]].Energy();
        if( i==1 || i==4 ) {
          etof2x1 += tofHitMap[nphi[i]].Energy();
        }
      }
 
      tofHitMap[nphi[i]].Energy(0);
    }
 
    double ecorr = iShowerMap->second.e5x5() + etof2x3/GeV;
    
    iShowerMap->second.ETof2x3(etof2x3/GeV);
    iShowerMap->second.ETof2x1(etof2x1/GeV);
    iShowerMap->second.setEnergy(ECorrection(ecorr));
  }
}
 
RecEmcEnergy EmcRecTofMatch::ECorrection(const RecEmcEnergy eIn)
{
  EmcRecParameter& Para=EmcRecParameter::GetInstance();
 
  RecEmcEnergy eOut=0;
  double par[4];
  for(int i=0;i<4;i++) {
    par[i]=Para.ECorr(i);
  }
  
  eOut = eIn/(par[0]+par[1]*eIn+par[2]*eIn*eIn+par[3]*eIn*eIn*eIn);
  return eOut;
}