EmcRecShowerPosLin.cxx

Go to the documentation of this file.

//
//    Linear weighted position attribute reconstruction
//
//    Created by Zhe Wang 2004, 5, 16
//
#include "EmcRec/EmcRecShowerPosLin.h"
#include "EmcRec/EmcRecParameter.h"
 
#include <fstream>
#include <iostream>
 
#include "CLHEP/Matrix/Matrix.h"
#include "EmcRecGeoSvc/IEmcRecGeoSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
 
 
 
void EmcRecShowerPosLin::Position( RecEmcShower& aShower ) {
  RecEmcFractionMap::const_iterator cit;
  HepPoint3D                        pos( 0, 0, 0 );
  HepPoint3D                        possum( 0, 0, 0 );
  double                            w, wsum;
  double                            etot;
  // cout<<"EmcRecShowerPosLin::Position Begin"<<endl;
 
  IEmcRecGeoSvc* iGeoSvc;
  ISvcLocator*   svcLocator = Gaudi::svcLocator();
  StatusCode     sc         = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
  if ( sc != StatusCode::SUCCESS )
  {
    //    cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
  }
 
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
 
  double e5x5 = aShower.e5x5();
 
  etot = 0;
 
  for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
  { etot += ( cit->second.getEnergy() * cit->second.getFraction() ); }
  wsum = 0;
  for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
  {
    w = ( cit->second.getEnergy() * cit->second.getFraction() );
    wsum += w;
    pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
    possum += pos * w;
  }
 
  if ( wsum <= 0 ) { cout << "[[[[[[[[[[[[[[[" << wsum; }
 
  pos = possum / wsum;
  //  aShower.setPosition(pos);
 
  // PosCorr=0 without position correction
  // PosCorr=1 with position correction
  if ( Para.PosCorr() == 0 ) { aShower.setPosition( pos ); }
 
  //----------------------------------------------------------------------
  // position correction
  RecEmcID           id     = aShower.getShowerId();
  const unsigned int module = EmcID::barrel_ec( id );
  ;
  const unsigned int thetaModule = EmcID::theta_module( id );
  const unsigned int phiModule   = EmcID::phi_module( id );
  HepPoint3D         posCorr( pos );
 
  if ( module == 1 )
  { // barrel
    unsigned int theModule;
    if ( thetaModule > 21 )
    {
      theModule = 43 - thetaModule;
      id        = EmcID::crystal_id( module, theModule, phiModule );
      posCorr.setZ( -posCorr.z() );
    }
    else { theModule = thetaModule; }
    // cout<<"EmcID: "<<id<<" theta: "<<theModule<<endl;
 
    EmcRecParameter& Para = EmcRecParameter::GetInstance();
    double           IRShift, parTheta[5], parPhi[5];
 
    if ( theModule == 21 ) { IRShift = 0; }
    else if ( theModule == 20 ) { IRShift = 2.5; }
    else { IRShift = 5.; }
 
    if ( Para.MethodMode() == 0 )
    {
 
      for ( int i = 0; i < 5; i++ )
      {
        if ( theModule == 21 )
        {
          double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
          parTheta[i]   = par[i];
        }
        else if ( theModule == 20 )
        {
          double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
          parTheta[i]   = par[i];
        }
        else
        {
          double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
          parTheta[i]   = par[i];
        }
      }
    }
    else if ( Para.MethodMode() == 1 )
    {
 
      for ( int i = 0; i < 5; i++ )
      {
 
        if ( e5x5 > 1.0 ) parTheta[i] = Para.BarrLinThetaPara( theModule, i );
        else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
          parTheta[i] = Para.BarrLinThetaPara( theModule + 22, i );
        else if ( e5x5 <= 0.5 ) parTheta[i] = Para.BarrLinThetaPara( theModule + 44, i );
 
        if ( e5x5 > 1.0 ) parPhi[i] = Para.BarrLinPhiPara( 0, i );
        else if ( e5x5 <= 1.0 && e5x5 > 0.5 ) parPhi[i] = Para.BarrLinPhiPara( 1, i );
        else if ( e5x5 <= 0.5 ) parPhi[i] = Para.BarrLinPhiPara( 2, i );
      }
    }
 
    HepPoint3D IR( 0, 0, IRShift );
    HepPoint3D center01, center23; // center of point0,1 and point2,3 in crystal
    center01 = ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
    center23 = ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
    // cout<<"before correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
 
    double theta01, theta23, length, d;
    theta01 = ( center01 - IR ).theta();
    theta23 = ( center23 - IR ).theta();
    length  = ( center01 - IR ).mag();
    d       = length * tan( theta01 - theta23 ); // length in x direction
    // cout<<"theta01: "<<theta01<<"\t"<<" theta23: "<<theta23<<"\t"
    //<<" length: "<<length<<" d: "<<d<<endl;
 
    double xIn, xOut, deltaTheta;
    xIn  = length * tan( theta01 - ( posCorr - IR ).theta() ) - d / 2;
    xOut = xIn - ( parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
                   parTheta[3] );
    deltaTheta = atan( ( xOut + d / 2 ) / length ) - atan( ( xIn + d / 2 ) / length );
    // cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
 
    double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 }; // e-
    double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 }; // e+
 
    //    double parPhi[5];// = { 0.7723, -3.1, 0.6992, -0.1441, -0.01012 }; //by photon, not
    //    used
 
    //    for(int i=0;i<5;i++) {
    //      parPhi[i] = (parPhi1[i]+parPhi2[i])/2;  //average of e- and e+
    //    }
 
    double yIn, yOut, deltaPhi;
    //    yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
    //      : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
    if ( phiModule == 0 && posCorr.phi() < 0 )
    { yIn = posCorr.phi() * 180. / CLHEP::pi + 360. - 119 * 3 - 1.843; }
    else if ( phiModule == 119 && posCorr.phi() > 0 )
    { yIn = posCorr.phi() * 180. / CLHEP::pi - 1.843; }
    else
    {
 
      yIn = posCorr.phi() < 0
                ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
                : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
    }
 
    yOut     = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
    deltaPhi = yOut * CLHEP::pi / 180.;
    // cout<<"yIn="<<yIn<<"\tyOut="<<yOut<<"\tdeltaPhi="<<deltaPhi<<endl;
 
    HepPoint3D rotateVector( 0, 1, 0 ); // y axis
    rotateVector.rotateZ( posCorr.phi() );
    posCorr.setZ( posCorr.z() - IRShift );
    posCorr.rotate( -deltaTheta, rotateVector );
    posCorr.setZ( posCorr.z() + IRShift );
    if ( Para.MethodMode() == 0 ) { posCorr.rotateZ( -0.002994 ); }
    else if ( Para.MethodMode() == 1 ) { posCorr.rotateZ( -deltaPhi ); }
 
    if ( thetaModule > 21 ) { posCorr.setZ( -posCorr.z() ); }
  }
  else
  { // endcap
 
    EmcRecParameter& Para = EmcRecParameter::GetInstance();
 
    if ( Para.MethodMode() == 1 )
    {
 
      double     IRShift = 0;
      HepPoint3D IR( 0, 0, IRShift );
      double theta12 = -9999, theta03 = -9999, theta23 = -9999, theta14 = -9999, delta = -9999;
      double phi = -9999, phi01 = -9999, phi23 = -9999, phi12 = -9999, phi34 = -9999,
             deltaphi = -9999;
      double xIn, xOut, deltaTheta, yIn, yOut, deltaPhi;
      double posphi = posCorr.phi();
 
      double parTheta[5], parPhi[5];
 
      HepPoint3D point0, point1, point2, point3, point4;
      point0 = iGeoSvc->GetCrystalPoint( id, 0 );
      point1 = iGeoSvc->GetCrystalPoint( id, 1 );
      point2 = iGeoSvc->GetCrystalPoint( id, 2 );
      point3 = iGeoSvc->GetCrystalPoint( id, 3 );
      point4 = iGeoSvc->GetCrystalPoint( id, 4 );
 
      HepPoint3D center, center01, center23, center12, center03, center34, center14;
      center01 = ( point0 + point1 ) / 2;
      center23 = ( point2 + point3 ) / 2;
      center12 = ( point1 + point2 ) / 2;
      center34 = ( point3 + point4 ) / 2;
      center03 = ( point0 + point3 ) / 2;
      center14 = ( point1 + point4 ) / 2;
 
      phi01 = center01.phi();
      phi23 = center23.phi();
      phi12 = center12.phi();
      phi34 = center34.phi();
 
      if ( ( thetaModule == 1 &&
             ( ( phiModule + 3 ) % 4 == 0 || ( phiModule + 2 ) % 4 == 0 ) ) ||
           ( thetaModule == 3 && ( ( phiModule + 4 ) % 5 == 0 || ( phiModule + 3 ) % 5 == 0 ||
                                   ( phiModule + 2 ) % 5 == 0 ) ) )
      {
        center12 = center23;
        center03 = center14;
 
        center01 = center12;
        center23 = center34;
 
        phi01 = phi12;
        phi23 = phi34;
      }
 
      if ( phiModule == 0 ) { posphi = posphi; }
      else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
                ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
                ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
      { posphi = posphi + 2 * CLHEP::pi; }
      else { posphi = posphi < 0 ? posphi + 2 * CLHEP::pi : posphi; }
 
      phi01 = phi01 < 0 ? phi01 + 2 * CLHEP::pi : phi01;
      phi23 = phi23 < 0 ? phi23 + 2 * CLHEP::pi : phi23;
 
      if ( module == 0 )
      { // east endcap
 
        IRShift = 10;
        center  = center23;
        phi     = phi23;
 
        for ( int i = 0; i < 5; i++ )
        {
          parTheta[i] = Para.EastLinThetaPara( thetaModule, i );
          parPhi[i]   = Para.EastLinPhiPara( 0, i );
        }
      }
      else if ( module == 2 )
      { // west endcap
 
        IRShift = -10;
        center  = center01;
        phi     = phi01;
 
        for ( int i = 0; i < 5; i++ )
        {
          parTheta[i] = Para.WestLinThetaPara( thetaModule, i );
          parPhi[i]   = Para.WestLinPhiPara( 0, i );
        }
      }
 
      theta12 = ( center12 - IR ).theta();
      theta03 = ( center03 - IR ).theta();
      delta   = theta03 - theta12;
      xIn     = ( ( ( posCorr - IR ).theta() - theta12 ) - delta * 0.5 ) / delta;
      xOut    = parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
             parTheta[3];
      deltaTheta = xOut * delta;
 
      deltaphi = fabs( phi23 - phi01 );
      yIn      = ( ( posphi - phi ) - deltaphi * 0.5 ) / deltaphi;
      yOut     = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
      deltaPhi = yOut * deltaphi;
 
      HepPoint3D rotateVector( 0, 1, 0 ); // y axis
      rotateVector.rotateZ( center.phi() );
      posCorr.setZ( posCorr.z() - IRShift );
      posCorr.rotate( -deltaTheta, rotateVector );
      posCorr.setZ( posCorr.z() + IRShift );
      posCorr.rotateZ( -deltaPhi );
    }
  }
 
  //  aShower.setPosition(posCorr);
  if ( Para.PosCorr() == 1 ) { aShower.setPosition( posCorr ); }
 
  //----------------------------------------------------------------------
  if ( aShower.energy() < 1e-5 ) return;
 
  double r, theta, phi;
  double dtheta, dphi, dx, dy, dz;
 
  r     = posCorr.mag();
  theta = posCorr.theta();
  phi   = posCorr.phi();
  // cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
 
  //  EmcRecParameter& Para=EmcRecParameter::GetInstance();
  if ( aShower.energy() > 0 )
  {
    dtheta = Para.SigTheta( 0 ) / sqrt( aShower.energy() ) + Para.SigTheta( 1 );
    dphi   = Para.SigPhi( 0 ) / sqrt( aShower.energy() ) + Para.SigPhi( 1 );
  }
  else
  {
    dtheta = 0;
    dphi   = 0;
  }
  // cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
 
  dx = fabs( iGeoSvc->GetBarrelR() * sin( phi ) * dphi );
  dy = fabs( iGeoSvc->GetBarrelR() * cos( phi ) * dphi );
  if ( theta > 0 )
    dz = fabs( iGeoSvc->GetBarrelR() * dtheta / ( sin( theta ) * sin( theta ) ) );
  else dz = 0;
  // cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
 
  aShower.setDtheta( dtheta );
  aShower.setDphi( dphi );
 
  //----------------------------------------------------------------------
  // Error matrix
  HepSymMatrix matrix( 3 ); // error matrix of r, theta, phi
  matrix[0][0] = 0;
  matrix[1][1] = dtheta * dtheta;
  matrix[2][2] = dphi * dphi;
  matrix[0][1] = 0;
  matrix[0][2] = 0;
  matrix[1][2] = 0;
  // cout<<"matrix: \n"<<matrix;
 
  HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
  matrix1[0][0] = sin( theta ) * cos( phi );
  matrix1[0][1] = r * cos( theta ) * cos( phi );
  matrix1[0][2] = -r * sin( theta ) * sin( phi );
  matrix1[1][0] = sin( theta ) * sin( phi );
  matrix1[1][1] = r * cos( theta ) * sin( phi );
  matrix1[1][2] = r * sin( theta ) * cos( phi );
  matrix1[2][0] = cos( theta );
  matrix1[2][1] = -r * sin( theta );
  matrix1[2][2] = 0;
  // cout<<"matrix1: \n"<<matrix1;
 
  for ( int i = 0; i < 3; i++ )
  {
    for ( int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
  }
  // cout<<"matrix2: \n"<<matrix2;
 
  HepMatrix matrix4( 3, 3 );
  matrix4 = matrix1 * matrix * matrix2;
  // cout<<"matrix4: \n"<<matrix4;
 
  HepSymMatrix matrix5( 3 ); // error matrix of x, y, z
  for ( int i = 0; i < 3; i++ )
  {
    for ( int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
  }
  // cout<<"matrix5: \n"<<matrix5;
 
  aShower.setErrorMatrix( matrix5 );
}