EmcRecSplitWeighted.cxx

Go to the documentation of this file.

//
//  Spliter
//  Simple weighted method
//
//  Zhe Wang 2004, 8, 31
//
#include <iostream>
 
#include "EmcRecGeoSvc/IEmcRecGeoSvc.h"
#include "EmcRec/EmcRecParameter.h"
#include "EmcRec/EmcRecShowerEnergy.h"
#include "EmcRec/EmcRecShowerPosLin.h"
#include "EmcRec/EmcRecShowerPosLinShMax.h"
#include "EmcRec/EmcRecShowerPosLog.h"
#include "EmcRec/EmcRecShowerPosLogShMax.h"
#include "EmcRec/EmcRecShowerPosLoglin.h"
#include "EmcRec/EmcRecShowerShape.h"
#include "EmcRec/EmcRecSplitWeighted.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
 
// Constructors and destructors
EmcRecSplitWeighted::EmcRecSplitWeighted() {
  fShowerE     = new EmcRecShowerEnergy;
  fShowerShape = new EmcRecShowerShape;
 
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
  //  if(Para.PositionMode1()=="lin") {
  //    fShowerPos=new EmcRecShowerPosLin;
  //  } else {
  //    fShowerPos=new EmcRecShowerPosLog;
  //  }
 
  if ( Para.PositionMode1() == "lin" ) { fShowerPos = new EmcRecShowerPosLin; }
  else if ( Para.PositionMode1() == "log" ) { fShowerPos = new EmcRecShowerPosLog; }
  else if ( Para.PositionMode1() == "loglin" ) { fShowerPos = new EmcRecShowerPosLoglin; }
  else if ( Para.PositionMode1() == "linsh" ) { fShowerPos = new EmcRecShowerPosLinShMax; }
  else if ( Para.PositionMode1() == "logsh" ) { fShowerPos = new EmcRecShowerPosLogShMax; }
}
 
EmcRecSplitWeighted::~EmcRecSplitWeighted() {
  delete fShowerE;
  delete fShowerPos;
  delete fShowerShape;
}
 
void EmcRecSplitWeighted::Split( RecEmcCluster& aCluster, const RecEmcIDVector& aMaxVec,
                                 RecEmcShowerMap& aShowerMap ) {
  // A exponential function will be used to describe the shape of a shower.
  // Energy fall off with distance from centre of shower is exp(-a*dist/r)
  // dist: distance from center
  // r: moliere radius
  // a: a parameter, needed to be optimised
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
 
  // Calculate the cut of second moment
  double eCluster = aCluster.getEnergy();
  double cut      = 0;
  if ( eCluster > Para.SmCut( 3 ) )
  { cut = Para.SmCut( 0 ) * exp( -( eCluster ) / Para.SmCut( 1 ) ) + Para.SmCut( 2 ); }
  else
  {
    cut = Para.SmCut( 0 ) * exp( -( Para.SmCut( 3 ) ) / Para.SmCut( 1 ) ) +
          Para.SmCut( 2 ); // constant
  }
  cut /= 10;
 
  RecEmcShower                 aShower;
  RecEmcHit                    aHit;
  RecEmcFraction               aFraction;
  RecEmcHitMap::const_iterator ciHitMap;
  if ( aMaxVec.size() == 0 )
  {
    // the cluster is abandoned
  }
 
  // only one seed or second moment belows cut
  else if ( aMaxVec.size() == 1 || aMaxVec.size() > 20 || aCluster.getSecondMoment() < cut )
  {
    aShower.Clear();
    aShower.setStatus( 1 );
    // ID
    aShower.ShowerId( aMaxVec[0] );
    aShower.ClusterId( aCluster.getClusterId() );
    aCluster.InsertShowerId( aMaxVec[0] );
    // each fraction
    double time = aCluster.Find( aMaxVec[0] )->second.getTime();
    if ( time >= Para.TimeMin() && time <= Para.TimeMax() )
    {
      for ( ciHitMap = aCluster.Begin(); ciHitMap != aCluster.End(); ++ciHitMap )
      {
        aHit = ciHitMap->second;
        // if((aHit.time()>=Para.TimeMin()) &&
        //(aHit.time()<=Para.TimeMax())){
        aFraction = aHit;
        aFraction.Fraction( 1 );
        aShower.Insert( aFraction );
        //}
      }
 
      aShower.setTime( time );
      fShowerE->Energy( aShower );
      fShowerPos->Position( aShower );
      fShowerShape->CalculateMoment( aShower );
      aShower.ThetaGap( 9 );
      aShower.PhiGap( 9 );
      // push it into map
      aShowerMap[aMaxVec[0]] = aShower;
    }
  }
 
  // two or more seeds and second moment beyonds cut
  else if ( aMaxVec.size() > 1 && aMaxVec.size() <= 20 && aCluster.getSecondMoment() > cut )
  {
    // cout<<"In EmcRecSplitWeighted: ";
    // cout<<aMaxVec.size();
    // cout<<" seeds in a cluster"<<endl;
 
    RecEmcCluster                tmpCluster = aCluster;
    RecEmcHitMap::const_iterator ci_HitMap;
    ci_RecEmcIDVector            ciMax, ciMax1;
    ;
    RecEmcShower              aShower;
    RecEmcShowerMap           tmpShowerMap;
    RecEmcShowerMap::iterator i_ShowerMap, i2_ShowerMap;
 
    RecEmcFraction aFrac;
    unsigned int   iterations = 0; // iteration time
    double         centroidShift;  // center shift between 2 iterations
    map<RecEmcID, HepPoint3D, less<RecEmcID>>                 showerCentroid; // shower center
    map<RecEmcID, HepPoint3D, less<RecEmcID>>::const_iterator ci_showerCentroid;
 
    IEmcRecGeoSvc* iGeoSvc;
    ISvcLocator*   svcLocator = Gaudi::svcLocator();
    StatusCode     sc         = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
    if ( sc != StatusCode::SUCCESS )
    {
      // cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
    }
 
    // Preparation
    for ( ciMax = aMaxVec.begin(); ciMax != aMaxVec.end(); ++ciMax )
    {
      // initialize: seed's front center
      showerCentroid[*ciMax] = iGeoSvc->GetCFrontCenter( *ciMax );
    }
    do {
      centroidShift = 0.;
      tmpShowerMap.clear();
      for ( ciMax = aMaxVec.begin(); ciMax != aMaxVec.end(); ++ciMax )
      {
        double weightTotal = 0., weight = 0.;
 
        aShower.Clear();
        aShower.ShowerId( *ciMax );
 
        // figure out the weight for each cell for each seed
        for ( ci_HitMap = tmpCluster.Begin(); ci_HitMap != tmpCluster.End(); ++ci_HitMap )
        {
 
          aFrac                  = ci_HitMap->second;
          double       aDistance = 0;
          RecEmcEnergy aESeed    = 0;
 
          bool isSeed = false;
          for ( ciMax1 = aMaxVec.begin(); ciMax1 != aMaxVec.end(); ++ciMax1 )
          {
            HepPoint3D seedPoint( showerCentroid[*ciMax1] );
            HepPoint3D hitPoint( iGeoSvc->GetCFrontCenter( aFrac.getCellId() ) );
 
            RecEmcEnergy theESeed = tmpCluster.Find( *ciMax1 )->second.getEnergy();
            double       theDistance;
            if ( *ciMax1 == aFrac.getCellId() )
            {
              isSeed      = true;
              theDistance = 0.;
            }
            else { theDistance = ( hitPoint - seedPoint ).mag(); }
 
            if ( *ciMax1 == *ciMax )
            {
              aDistance = theDistance;
              aESeed    = theESeed;
            }
 
            weightTotal +=
                theESeed * exp( -Para.LateralProfile() * theDistance / Para.MoliereRadius() );
          }
 
          weight = aESeed * exp( -Para.LateralProfile() * aDistance / Para.MoliereRadius() ) /
                   weightTotal;
          aFrac.Fraction( weight );
          // if((aFrac.time()>=Para.TimeMin() &&
          // aFrac.time()<=Para.TimeMax()) ||
          // isSeed) {
          if ( aFrac.getEnergy() * aFrac.getFraction() > Para.ElectronicsNoiseLevel() )
          { aShower.Insert( aFrac ); }
          //}
          weightTotal = 0;
        }
 
        fShowerE->Energy( aShower );
        fShowerPos->Position( aShower );
        HepPoint3D newCentroid( aShower.position() );
        HepPoint3D oldCentroid( showerCentroid[*ciMax] );
        centroidShift += ( newCentroid - oldCentroid ).mag();
        tmpShowerMap[aShower.getShowerId()] = aShower;
        showerCentroid[*ciMax]              = newCentroid;
      }
 
      centroidShift /= (double)aMaxVec.size();
      for ( ci_showerCentroid = showerCentroid.begin();
            ci_showerCentroid != showerCentroid.end(); ci_showerCentroid++ )
      {
        showerCentroid[ci_showerCentroid->first] =
            tmpShowerMap[ci_showerCentroid->first].position();
      }
      iterations++;
      // cout<<"iterations: "<<iterations<<"\tcentroidShift: "<<centroidShift<<endl;
    } while ( ( iterations < 8 ) && ( centroidShift > 0.01 ) );
 
    unsigned int tht, phi;
    unsigned int tht2, phi2;
    unsigned int dtht, dphi;
    unsigned int thetagap = 0, phigap = 0;
    double       distmin, dist;
    RecEmcID     id, id2, nearest;
    // shower attributes
    for ( i_ShowerMap = tmpShowerMap.begin(); i_ShowerMap != tmpShowerMap.end();
          ++i_ShowerMap )
    {
      aCluster.InsertShowerId( i_ShowerMap->second.getShowerId() );
      i_ShowerMap->second.ClusterId( aCluster.getClusterId() );
      i_ShowerMap->second.setStatus( 2 );
      fShowerE->Energy( i_ShowerMap->second );
      fShowerPos->Position( i_ShowerMap->second );
      fShowerShape->CalculateMoment( i_ShowerMap->second );
      //
      // min theta, phi gap
      id      = ( i_ShowerMap->second ).getShowerId();
      tht     = EmcID::theta_module( id );
      phi     = EmcID::phi_module( id );
      distmin = 999999;
      for ( i2_ShowerMap = tmpShowerMap.begin(); i2_ShowerMap != tmpShowerMap.end();
            ++i2_ShowerMap )
      {
        id2  = ( i2_ShowerMap->second ).getShowerId();
        tht2 = EmcID::theta_module( id2 );
        phi2 = EmcID::phi_module( id2 );
        if ( id != id2 )
        {
          dtht = tht > tht2 ? tht - tht2 : tht2 - tht;
          dphi = phi > phi2 ? phi - phi2 : phi2 - phi;
          if ( dphi > ( EmcID::getPHI_BARREL_MAX() + 1 ) / 2 )
            dphi = EmcID::getPHI_BARREL_MAX() + 1 - dphi;
          dist = sqrt( double( dtht * dtht + dphi * dphi ) );
          if ( dist < distmin )
          {
            distmin = dist;
            nearest = id2;
            if ( dtht > 6 ) dtht = 6;
            if ( dphi > 6 ) dphi = 6;
            thetagap = dtht;
            phigap   = dphi;
          }
        }
      }
      // cout<<"+++++"<<id<<" "<<thetagap<<" "<<phigap<<endl;
      i_ShowerMap->second.NearestSeed( nearest );
      i_ShowerMap->second.ThetaGap( thetagap );
      i_ShowerMap->second.PhiGap( phigap );
 
      // save it
      double time =
          i_ShowerMap->second.Find( i_ShowerMap->second.getShowerId() )->second.getTime();
      if ( time >= Para.TimeMin() && time <= Para.TimeMax() &&
           ( i_ShowerMap->second.getEAll() > Para.EThresholdCluster() ) )
      {
        i_ShowerMap->second.setTime( time );
        aShowerMap[i_ShowerMap->first] = i_ShowerMap->second;
      }
    }
    tmpShowerMap.clear();
  }
}