EmcSelBhaEvent Class Reference

#include <EmcSelBhaEvent.h>

Inheritance diagram for EmcSelBhaEvent:

Public Types
enum	{ m_oneProng = 1 , m_twoProng = 2 }

Public Member Functions
	EmcSelBhaEvent (const std::string &name, ISvcLocator *pSvcLocator)
	~EmcSelBhaEvent ()
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
bool	passed ()
void	setPassed (bool passed)
int	selectedType () const
int	selectedTrkID1 () const
int	selectedTrkID2 () const
int	index (int theta, int phi) const
void	initGeom ()
void	SelectBhabha ()
StatusCode	SelectFillBhabha ()
void	FillBhabha ()
void	CollectBhabha ()
void	OutputMV ()
double	findPhiDiff (double phi1, double phi2)
void	readCorFun ()
void	readEsigma ()
void	readDepEneFun ()
void	readSigmaExp ()
void	readRawPeakIxtal ()
void	readDataAndMCIxtal ()
double	Angle2ClosestShower (int ShowerID)

Detailed Description

Definition at line 54 of file EmcSelBhaEvent.h.

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
m_oneProng
m_twoProng

Definition at line 58 of file EmcSelBhaEvent.h.

{ m_oneProng = 1, m_twoProng = 2 };

Constructor & Destructor Documentation

EmcSelBhaEvent()

EmcSelBhaEvent::EmcSelBhaEvent	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 75 of file EmcSelBhaEvent.cxx.

    : Algorithm( name, pSvcLocator )
    , m_vr0cut( 1.0 )
    , m_vz0cut( 5.0 )
    , m_lowEnergyShowerCut( 0.1 )
    , m_highEnergyShowerCut( 0.5 )
    , m_matchThetaCut( 0.2 )
    , m_matchPhiCut( 0.2 )
    , m_highMomentumCut( 0.5 )
    , m_EoPMaxCut( 1.3 )
    , m_EoPMinCut( 0.6 )
    , m_minAngShEnergyCut( 0.2 )
    , m_minAngCut( 0.3 )
    , m_acolliCut( 0.03 )
    , m_eNormCut( 0.5 )
    , m_pNormCut( 0.5 )
    , m_oneProngMomentumCut( 1.2 )
    ,
 
    m_digiRangeCut( 6 )
    , m_ShEneThreshCut( 0.02 )
    , m_ShEneLeptonCut( 1. )
    , m_minNrXtalsShowerCut( 10 )
    , m_maxNrXtalsShowerCut( 70 )
    , m_phiDiffMinCut( 0.05 )
    , m_phiDiffMaxCut( 0.2 )
    , m_nrShThreshCut( 20 )
    , m_thetaDiffCut( 0.04 )
    , m_LATCut( 0.8 )
    ,
 
    m_showersAccepted( 0 )
    ,
    //--------------------
    m_writeMVToFile( true )
    , m_fileExt( "" )
    , m_inputFileDir( "../InputData/" )
    , m_fileDir( "/home/besdata/public/liucx/Calib/" )
    , m_selMethod( "Ixtal" )
    , m_nXtals( 6240 )
    , m_sigmaCut( 1. )
    , m_beamEnergy( 1.843 )
    , m_MsgFlag( 0 )
    , m_Num( 0 )
    , m_fileNum( 0 )
 
{
 
  // Declare the properties
 
  declareProperty( "Vr0cut", m_vr0cut = 1.0 ); // suggest cut: |z0|<5cm && r0<1cm
  declareProperty( "Vz0cut", m_vz0cut = 5.0);
 
  declareProperty( "lowEnergyShowerCut", m_lowEnergyShowerCut );
  declareProperty( "highEnergyShowerCut", m_highEnergyShowerCut );
  declareProperty( "matchThetaCut", m_matchThetaCut );
  declareProperty( "matchPhiCut", m_matchPhiCut );
 
  declareProperty( "highMomentumCut", m_highMomentumCut );
  declareProperty( "EoPMaxCut", m_EoPMaxCut );
  declareProperty( "EoPMinCut", m_EoPMinCut );
  declareProperty( "minAngShEnergyCut", m_minAngShEnergyCut );
  declareProperty( "minAngCut", m_minAngCut );
  declareProperty( "acolliCut", m_acolliCut );
  declareProperty( "eNormCut", m_eNormCut );
  declareProperty( "pNormCut", m_pNormCut );
  declareProperty( "oneProngMomentumCut", m_oneProngMomentumCut );
 
  // *
 
  declareProperty( "digiRangeCut", m_digiRangeCut );
 
  declareProperty( "ShEneThreshCut", m_ShEneThreshCut );
  declareProperty( "ShEneLeptonCut", m_ShEneLeptonCut );
 
  declareProperty( "minNrXtalsShowerCut", m_minNrXtalsShowerCut );
  declareProperty( "maxNrXtalsShowerCut", m_maxNrXtalsShowerCut );
  declareProperty( "phiDiffMinCut", m_phiDiffMinCut );
  declareProperty( "phiDiffMaxCut", m_phiDiffMaxCut );
  declareProperty( "nrShThreshCut", m_nrShThreshCut );
  declareProperty( "thetaDiffCut", m_thetaDiffCut );
  declareProperty( "LATCut", m_LATCut );
 
  //--------------------
  declareProperty( "writeMVToFile", m_writeMVToFile );
  declareProperty( "fileExt", m_fileExt );
  declareProperty( "fileDir", m_fileDir );
  declareProperty( "inputFileDir", m_inputFileDir );
  declareProperty( "selMethod", m_selMethod );
  declareProperty( "sigmaCut", m_sigmaCut );
  declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
 
  declareProperty( "beamEnergy", m_beamEnergy );
  declareProperty( "Dbeam", m_dbeam = 0.00 );
 
  declareProperty( "UseBeamDetectorEnergy", m_useBeamDetectorEnergy = false );
 
  declareProperty( "elecSaturation", m_elecSaturation = true );
 
  declareProperty( "MsgFlag", m_MsgFlag );
  declareProperty( "Num", m_Num );
  declareProperty( "fileNum", m_fileNum );
 
  int j   = 0;
  m_index = new int*[56];
  for ( j = 0; j < 56; j++ )
  {
    m_index[j] = new int[120];
    for ( int k = 0; k < 120; k++ ) { m_index[j][k] = -1; }
  }
 
  m_iGeoSvc = 0;
 
  for ( int i = 0; i < 6240; i++ ) { m_inputConst[i] = 1.0; }
 
  m_irun = 0;
}

Referenced by EmcSelBhaEvent().

~EmcSelBhaEvent()

EmcSelBhaEvent::~EmcSelBhaEvent

(

)

Definition at line 196 of file EmcSelBhaEvent.cxx.

                                {
 
  if ( m_index != 0 )
  {
    for ( int i = 0; i < 56; i++ ) delete[] m_index[i];
    delete[] m_index;
    m_index = 0;
  }
 
  ///////////
  for ( int j = 0; j < 6240; j++ ) { m_measure[j] = -1; }
}

Member Function Documentation

Angle2ClosestShower()

double EmcSelBhaEvent::Angle2ClosestShower

(

int

ShowerID

)

Definition at line 1594 of file EmcSelBhaEvent.cxx.

                                                         {
 
  double minDist = 999;
 
  SmartDataPtr<EvtRecEvent>    evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + ShowerID;
 
  RecEmcShower*    theShower = ( *itTrk )->emcShower();
  HepLorentzVector theShowerVec( 1, 1, 1, 1 );
  theShowerVec.setTheta( theShower->theta() );
  theShowerVec.setPhi( theShower->phi() );
  theShowerVec.setRho( theShower->energy() );
  theShowerVec.setE( theShower->energy() );
 
  for ( int j = 0; j < evtRecEvent->totalTracks(); j++ )
  {
    EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
 
    if ( !( *jtTrk )->isEmcShowerValid() ) continue;
    if ( ShowerID == ( *jtTrk )->trackId() ) continue;
 
    RecEmcShower* aShower = ( *jtTrk )->emcShower();
 
    if ( aShower->energy() > m_minAngShEnergyCut )
    {
 
      HepLorentzVector aShowerVec( 1, 1, 1, 1 );
      aShowerVec.setTheta( aShower->theta() );
      aShowerVec.setPhi( aShower->phi() );
      aShowerVec.setRho( aShower->energy() );
      aShowerVec.setE( aShower->energy() );
 
      double dist = theShowerVec.angle( aShowerVec );
 
      if ( dist < minDist ) minDist = dist;
    }
  }
 
  return minDist;
}

CollectBhabha()

void EmcSelBhaEvent::CollectBhabha

(

)

Definition at line 860 of file EmcSelBhaEvent.cxx.

                                   {
 
  MsgStream log( msgSvc(), name() );
 
  // check if the Bhabhas were found
  if ( 0 != myBhaEvt->positron()->found() || 0 != myBhaEvt->electron()->found() )
  {
 
    m_taken++;
    // fill the EmcBhabhas
    double calibEnergy = 0.;
    double energyError = 0.;
 
    //------------- electron found --------------------------------------
    if ( myBhaEvt->electron()->found() )
    {
      /*
      int ixtalIndex = index(myBhaEvt->electron()->thetaIndex(),
                            myBhaEvt->electron()->phiIndex());
      calibEnergy = m_eMcPeak[ixtalIndex];
      */
 
      unsigned int module, ntheta, nphi;
      if ( myBhaEvt->electron()->thetaIndex() <= 5 )
      {
        module = 0;
        ntheta = myBhaEvt->electron()->thetaIndex();
      }
      else if ( myBhaEvt->electron()->thetaIndex() >= 50 )
      {
        module = 2;
        ntheta = 55 - myBhaEvt->electron()->thetaIndex();
      }
      else
      {
        module = 1;
        ntheta = myBhaEvt->electron()->thetaIndex() - 6;
      }
      nphi = myBhaEvt->electron()->phiIndex();
 
      RecEmcID   shId = EmcID::crystal_id( module, ntheta, nphi );
      HepPoint3D SeedPos( 0, 0, 0 );
      SeedPos = m_iGeoSvc->GetCFrontCenter( shId );
      /*
      calibEnergy = myBhaEvt->
        getDepoMCShowerEnergy_lab(myBhaEvt->electron()->theta(),
                                  myBhaEvt->electron()->phi(),
                                  myBhaEvt->electron()->thetaIndex(),
                                  myBhaEvt->electron()->phiIndex(),
                                  m_corFun[myBhaEvt->electron()->thetaIndex()],
                                  m_beamEnergy);
      */
      calibEnergy = myBhaEvt->getDepoMCShowerEnergy_lab(
          SeedPos.theta(), SeedPos.phi(), myBhaEvt->electron()->thetaIndex(),
          myBhaEvt->electron()->phiIndex(), m_corFun[myBhaEvt->electron()->thetaIndex()],
          m_beamEnergy );
 
      if ( m_selMethod == "DataMCIxtal" )
      {
        calibEnergy = myBhaEvt->getDepoMCShowerEnergy_lab(
            SeedPos.theta(), SeedPos.phi(), myBhaEvt->electron()->thetaIndex(),
            myBhaEvt->electron()->phiIndex(),
            m_eMeanMC[index( myBhaEvt->electron()->thetaIndex(),
                             myBhaEvt->electron()->phiIndex() )],
            m_beamEnergy );
      }
      /*
      calibEnergy = myBhaEvt->
         getDepoMCShowerEnergy(myBhaEvt->electron()->thetaIndex(),
                               myBhaEvt->electron()->phiIndex(),
                               m_corFun[myBhaEvt->electron()->thetaIndex()],
                               m_beamEnergy);
       */
 
      if ( calibEnergy > 0 )
      {
 
        energyError = myBhaEvt->getErrorDepoMCShowerEnergy(
            myBhaEvt->electron()->thetaIndex(), myBhaEvt->electron()->phiIndex(),
            m_eSigma[myBhaEvt->electron()->thetaIndex()] * m_beamEnergy / 100. );
        if ( m_selMethod == "DataMCIxtal" )
        {
          energyError = m_eRmsMC[index( myBhaEvt->electron()->thetaIndex(),
                                        myBhaEvt->electron()->phiIndex() )] *
                        m_beamEnergy;
        }
      }
      else
        log << MSG::WARNING << " Did not find MC deposited cluster energy "
            << " for this cluster:  thetaInd: " << myBhaEvt->electron()->thetaIndex()
            << "  phiInd: " << myBhaEvt->electron()->phiIndex() << endl
            << "Will not use this cluster for the Emc "
            << "Bhabha calibration !" << endmsg;
 
      // set all that stuff in an EmcBhabha
      myBhaEvt->setElectron()->setErrorOnCalibEnergy( energyError );
      myBhaEvt->setElectron()->setCalibEnergy( calibEnergy );
 
      // myBhaEvt->electron()->print();
    }
    else log << MSG::INFO << name() << ": Electron was not selected ! " << endmsg;
 
    calibEnergy = 0.;
    energyError = 0.;
 
    //---------------- positron found ----------------------------------
    if ( myBhaEvt->positron()->found() )
    {
      /*
      int ixtalIndex = index(myBhaEvt->positron()->thetaIndex(),
                            myBhaEvt->positron()->phiIndex());
      calibEnergy = m_eMcPeak[ixtalIndex];
      */
 
      unsigned int module, ntheta, nphi;
      if ( myBhaEvt->positron()->thetaIndex() <= 5 )
      {
        module = 0;
        ntheta = myBhaEvt->positron()->thetaIndex();
      }
      else if ( myBhaEvt->positron()->thetaIndex() >= 50 )
      {
        module = 2;
        ntheta = 55 - myBhaEvt->positron()->thetaIndex();
      }
      else
      {
        module = 1;
        ntheta = myBhaEvt->positron()->thetaIndex() - 6;
      }
      nphi = myBhaEvt->positron()->phiIndex();
 
      RecEmcID   shId = EmcID::crystal_id( module, ntheta, nphi );
      HepPoint3D SeedPos( 0, 0, 0 );
      SeedPos = m_iGeoSvc->GetCFrontCenter( shId );
      /*
     calibEnergy = myBhaEvt->
        getDepoMCShowerEnergy_lab(myBhaEvt->positron()->theta(),
                                  myBhaEvt->positron()->phi(),
                                  myBhaEvt->positron()->thetaIndex(),
                                  myBhaEvt->positron()->phiIndex(),
                                  m_corFun[myBhaEvt->positron()->thetaIndex()],
                                  m_beamEnergy);
      */
      calibEnergy = myBhaEvt->getDepoMCShowerEnergy_lab(
          SeedPos.theta(), SeedPos.phi(), myBhaEvt->positron()->thetaIndex(),
          myBhaEvt->positron()->phiIndex(), m_corFun[myBhaEvt->positron()->thetaIndex()],
          m_beamEnergy );
 
      if ( m_selMethod == "DataMCIxtal" )
      {
        calibEnergy = myBhaEvt->getDepoMCShowerEnergy_lab(
            SeedPos.theta(), SeedPos.phi(), myBhaEvt->positron()->thetaIndex(),
            myBhaEvt->positron()->phiIndex(),
            m_eMeanMC[index( myBhaEvt->positron()->thetaIndex(),
                             myBhaEvt->positron()->phiIndex() )],
            m_beamEnergy );
      }
 
      /*
      calibEnergy = myBhaEvt->
        getDepoMCShowerEnergy(myBhaEvt->positron()->thetaIndex(),
                              myBhaEvt->positron()->phiIndex(),
                              m_corFun[myBhaEvt->positron()->thetaIndex()],
                              m_beamEnergy);
      */
 
      if ( calibEnergy > 0. )
      {
 
        energyError = myBhaEvt->getErrorDepoMCShowerEnergy(
            myBhaEvt->positron()->thetaIndex(), myBhaEvt->positron()->phiIndex(),
            m_eSigma[myBhaEvt->positron()->thetaIndex()] * m_beamEnergy / 100. );
 
        if ( m_selMethod == "DataMCIxtal" )
        {
          energyError = m_eRmsMC[index( myBhaEvt->electron()->thetaIndex(),
                                        myBhaEvt->electron()->phiIndex() )] *
                        m_beamEnergy;
        }
      }
      else
        log << MSG::WARNING << " Did not find MC deposited cluster energy "
            << "for this cluster:  thetaInd: " << myBhaEvt->positron()->thetaIndex()
            << "  phiInd: " << myBhaEvt->positron()->phiIndex() << endl
            << "Will not use this cluster for the Emc "
            << "Bhabha calibration !" << endmsg;
 
      // set all that stuff in an EmcBhabha
      myBhaEvt->setPositron()->setErrorOnCalibEnergy( energyError );
      myBhaEvt->setPositron()->setCalibEnergy( calibEnergy );
 
      // myBhaEvt->positron()->print();
    }
    else log << MSG::INFO << name() << ": Positron was not selected ! " << endmsg;
 
    // calculate elements of Matrix M and vector R from Bhabha event,
    // M is in SLAP triad format
    fillMatrix();
  }
  else { log << MSG::WARNING << " No Bhabha data for calibration found in event !" << endmsg; }
}

Referenced by execute().

execute()

StatusCode EmcSelBhaEvent::execute

(

)

Definition at line 289 of file EmcSelBhaEvent.cxx.

                                   {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::DEBUG << "in execute()" << endmsg;
 
  // create the object that store the needed data of the Bhabha event
 
  int event, run;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
 
  run   = eventHeader->runNumber();
  event = eventHeader->eventNumber();
  // cout<<"---------"<<event<<".........."<<run<<endl;
  m_events++;
  m_run = run;
 
  //////////////////
  // get beam energy
  //////////////////
  if ( m_ReadBeamEFromDB && m_irun != run )
  {
    m_irun = run;
    m_BeamEnergySvc->getBeamEnergyInfo();
    if ( m_BeamEnergySvc->isRunValid() ) m_beamEnergy = m_BeamEnergySvc->getbeamE();
    // if(m_readDb.isRunValid(m_irun))
    //  m_beamEnergy=m_readDb.getbeamE(m_irun);
    // if(m_BeamEnergySvc->isRunValid())
    //  m_beamEnergy=m_BeamEnergySvc->getbeamE();
    cout << "beamEnergy=" << m_beamEnergy << endl;
    m_beamEnergy = m_beamEnergy - m_dbeam;
 
    double           the = 0.011;
    double           phi = 0;
    HepLorentzVector ptrk;
    ptrk.setPx( m_beamEnergy * sin( the ) * cos( phi ) );
    ptrk.setPy( m_beamEnergy * sin( the ) * sin( phi ) );
    ptrk.setPz( m_beamEnergy * cos( the ) );
    ptrk.setE( m_beamEnergy );
 
    ptrk = ptrk.boost( -0.011, 0, 0 );
 
    cout << "beamEnergy=" << m_beamEnergy << "  cms " << ptrk.e()
         << "  ratio=" << ( m_beamEnergy - ptrk.e() ) / ptrk.e() << endl;
    m_beamEnergy = ptrk.e();
  }
 
  if ( m_useBeamDetectorEnergy )
  {
 
    IScanEnergySvc* m_ScanEnergySvc = 0;
    StatusCode      scScan;
    scScan = Gaudi::svcLocator()->service( "ScanEnergySvc", m_ScanEnergySvc );
    if ( scScan != StatusCode::SUCCESS ) { cout << "can not use ScanEnergySvc" << endmsg; }
    else
    {
 
      std::cout << "Test ScanEnergySvc  getScanEnergy= " << m_ScanEnergySvc->getScanEnergy()
                << std::endl;
      m_beamEnergy = m_ScanEnergySvc->getScanEnergy() / 2.0 / 1000.;
      cout << "ScanbeamEnergy=" << m_beamEnergy << endl;
    }
  }
 
  ////////////
  // int mmea=0;
 
  for ( int j = 0; j < 6240; j++ ) { m_measure[j] = -1; }
 
  if ( m_elecSaturation == true )
  {
 
    ///////////////////////////find Measure/////////////
    SmartDataPtr<EmcDigiCol> emcDigiCol( eventSvc(), "/Event/Digi/EmcDigiCol" );
    if ( !emcDigiCol )
    {
      log << MSG::FATAL << "Could not find EMC digi" << endmsg;
      return ( StatusCode::FAILURE );
    }
 
    EmcDigiCol::iterator iter3;
    for ( iter3 = emcDigiCol->begin(); iter3 != emcDigiCol->end(); iter3++ )
    {
      RecEmcID id( ( *( iter3 ) )->identify() );
      // cout<<id<<endl;
 
      unsigned int npart, ntheta, nphi;
      npart  = EmcID::barrel_ec( id );
      ntheta = EmcID::theta_module( id );
      nphi   = EmcID::phi_module( id );
 
      unsigned int newthetaInd = -999;
      if ( npart == 0 ) newthetaInd = ntheta;
      if ( npart == 1 ) newthetaInd = ntheta + 6;
      if ( npart == 2 ) newthetaInd = 55 - ntheta;
 
      int ixtal        = index( newthetaInd, nphi );
      m_measure[ixtal] = ( *iter3 )->getMeasure();
      // if ((*iter3)->getMeasure()==3) mmea=9;
    }
  }
 
  ////////////
 
  myBhaEvt = new EmcBhabhaEvent();
 
  // Select Bhabha
  SelectBhabha();
  if ( m_selectedType == BhabhaType::OneProng ) m_OneProng++;
  // number of events with   TwoProngMatched
  if ( m_selectedType == BhabhaType::TwoProngMatched ) m_TwoProngMatched++;
  // number of events with   TwoProngOneMatched
  if ( m_selectedType == BhabhaType::TwoProngOneMatched ) m_TwoProngOneMatched++;
 
  if ( m_selectedType == BhabhaType::Nothing ) { m_Nothing++; }
 
  // retreive shower list of event
 
  if ( m_selectedType == BhabhaType::TwoProngMatched )
  {
    FillBhabha();
 
    // collect bhabha event for digi-calibration of EMC
    // and fill matrix and vector of system of linear equations
    CollectBhabha();
  }
 
  delete myBhaEvt;
  myBhaEvt = 0;
 
  return StatusCode::SUCCESS;
}

FillBhabha()

void EmcSelBhaEvent::FillBhabha

(

)

Definition at line 609 of file EmcSelBhaEvent.cxx.

                                {
 
  MsgStream log( msgSvc(), name() );
 
  SmartDataPtr<EvtRecEvent>    evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  EvtRecTrackIterator itTrk1 = evtRecTrkCol->begin() + m_selectedTrkID1;
 
  RecEmcShower* theShower1 = ( *itTrk1 )->emcShower();
  // RecMdcTrack *theMdcTrk1 = (*itTrk1)->mdcTrack();
 
  EvtRecTrackIterator itTrk2 = evtRecTrkCol->begin() + m_selectedTrkID2;
 
  RecEmcShower* theShower2 = ( *itTrk2 )->emcShower();
  // RecMdcTrack *theMdcTrk2 = (*itTrk2)->mdcTrack();
  EvtRecTrackIterator itTrk;
  RecEmcShower*       theShower;
 
  // * * * * * * * * * * * * * * * * * * * * * * * * * ** *  ** *
  // loop all showers of an event set EmcShower and EmcShDigi
 
  EmcShower* aShower = new EmcShower();
 
  double   ene, theta, phi, eseed;
  double   showerX, showerY, showerZ;
  long int thetaIndex, phiIndex, numberOfDigis;
 
  RecEmcID     showerId;
  unsigned int showerModule;
 
  HepPoint3D showerPosition( 0, 0, 0 );
 
  if ( !m_showerList.empty() ) m_showerList.clear();
 
  for ( int ish = 0; ish < 2; ish++ )
  {
 
    if ( ish == 0 )
    {
      itTrk     = itTrk1;
      theShower = theShower1;
    }
    if ( ish == 1 )
    {
      itTrk     = itTrk2;
      theShower = theShower2;
    }
    // ene=theShower->energy(); //corrected shower energy unit GeV
    ene   = theShower->e5x5();  // uncorrected 5x5 energy unit GeV
    eseed = theShower->eSeed(); // unit GeV
 
    /////////////////
    /*
   RecExtTrack *extTrk = (*itTrk)->extTrack();
 
      Hep3Vector extpos = extTrk->emcPosition();
 
      cout<<"Extpos="<<extpos<<endl;
 
 
      RecEmcShower *theShower = (*itTrk)->emcShower();
 
      Hep3Vector emcpos(theShower->x(), theShower->y(), theShower->z());
 
      cout<<"emcpos="<<emcpos<<endl;
 
      RecEmcID shId= theShower->getShowerId();
      unsigned int nprt= EmcID::barrel_ec(shId);
      //RecEmcID cellId= theShower->cellId();
      HepPoint3D SeedPos(0,0,0);
      SeedPos=m_iGeoSvc->GetCFrontCenter(shId);
      cout<<"SeedPos="<<SeedPos<<endl;
    */
    //////////////////
 
    showerPosition = theShower->position();
    theta          = theShower->theta();
    phi            = theShower->phi();
    showerX        = theShower->x();
    showerY        = theShower->y();
    showerZ        = theShower->z();
 
    showerId     = theShower->getShowerId();
    showerModule = EmcID::barrel_ec( showerId );
 
    // The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5)
    // module is defined  by Endcap_east(0),Barrel(1),Endcap_west(2)
 
    thetaIndex = EmcID::theta_module( showerId );
 
    phiIndex = EmcID::phi_module( showerId );
 
    // cout<<thetaIndex<<"   "<<phiIndex<<endl;
    //-------------------
 
    EmcShDigi*      aShDigi = new EmcShDigi();
    EmcShDigi       maxima  = EmcShDigi();
    list<EmcShDigi> shDigiList;
    RecEmcID        cellId;
    unsigned int    module;
    double          digiEne, time, fraction, digiTheta, digiPhi;
    double          digiX, digiY, digiZ;
    long int        digiThetaIndex, digiPhiIndex;
    HepPoint3D      digiPos( 0, 0, 0 );
 
    RecEmcFractionMap::const_iterator ciFractionMap;
 
    if ( !shDigiList.empty() ) shDigiList.clear();
    RecEmcFractionMap fracMap5x5 = theShower->getFractionMap5x5();
    for ( ciFractionMap = fracMap5x5.begin(); ciFractionMap != fracMap5x5.end();
          ciFractionMap++ )
    {
 
      digiEne = ciFractionMap->second.getEnergy(); // digit energy unit GeV
 
      time     = ciFractionMap->second.getTime();
      fraction = ciFractionMap->second.getFraction();
 
      cellId = ciFractionMap->second.getCellId();
 
      digiPos   = m_iGeoSvc->GetCFrontCenter( cellId );
      digiTheta = digiPos.theta();
      digiPhi   = digiPos.phi();
      digiX     = digiPos.x();
      digiY     = digiPos.y();
      digiZ     = digiPos.z();
 
      module = EmcID::barrel_ec( cellId );
      // The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5)
      // module is defined  by Endcap_east(0),Barrel(1),Endcap_west(2)
 
      digiThetaIndex = EmcID::theta_module( cellId );
 
      digiPhiIndex = EmcID::phi_module( cellId );
 
      // set EmcShDigi
      aShDigi->setEnergy( digiEne );
      aShDigi->setTheta( digiTheta );
      aShDigi->setPhi( digiPhi );
      aShDigi->setModule( module );
      aShDigi->setThetaIndex( digiThetaIndex );
      aShDigi->setPhiIndex( digiPhiIndex );
      aShDigi->setTime( time );
      aShDigi->setFraction( fraction );
      aShDigi->setWhere( digiPos );
      aShDigi->setY( digiX );
      aShDigi->setY( digiY );
      aShDigi->setZ( digiZ );
 
      shDigiList.push_back( *aShDigi );
    }
    shDigiList.sort(); // sort energy from small to large
 
    numberOfDigis = shDigiList.size();
 
    maxima = *( --shDigiList.end() );
    // cout<<"maxima="<<maxima.energy()<<endl;
    // cout<<maxima.thetaIndex()<<"  max  "<<maxima.phiIndex()<<endl;
    // set EmcShower
    aShower->setEnergy( ene );
    aShower->setTheta( theta );
    aShower->setPhi( phi );
    aShower->setModule( showerModule );
    aShower->setThetaIndex( thetaIndex );
    aShower->setPhiIndex( phiIndex );
    aShower->setNumberOfDigis( numberOfDigis );
    aShower->setDigiList( shDigiList );
    aShower->setMaxima( maxima );
    aShower->setWhere( showerPosition );
    aShower->setX( showerX );
    aShower->setY( showerY );
    aShower->setZ( showerZ );
    m_showerList.push_back( *aShower );
    delete aShDigi;
  }
  // m_showerList.sort();         //sort energy from small to large
 
  delete aShower;
 
  ///////////////
 
  if ( m_showerList.size() == 2 )
  {
    // iterator for the bump list
    list<EmcShower>::const_iterator iter_Sh;
    int                             itrk = 0;
    for ( iter_Sh = m_showerList.begin(); iter_Sh != m_showerList.end(); iter_Sh++ )
    {
 
      itrk++;
      EmcShower shower;
      shower       = EmcShower();
      double theta = iter_Sh->theta();
      shower       = *iter_Sh;
 
      // fill the Bhabha event
      //  if ( (itrk==1&&theMdcTrk1->charge()>0)
      //    ||(itrk==2&&theMdcTrk2->charge()>0) ){
      if ( itrk == 1 )
      {
        // set properties
        myBhaEvt->setPositron()->setShower( shower );
        myBhaEvt->setPositron()->setTheta( shower.theta() );
        myBhaEvt->setPositron()->setPhi( shower.phi() );
        unsigned int newthetaInd = -999;
        // The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
        // in Emc Bhabha Calibration
        if ( shower.module() == 0 ) newthetaInd = shower.thetaIndex();
        if ( shower.module() == 1 ) newthetaInd = shower.thetaIndex() + 6;
        if ( shower.module() == 2 ) newthetaInd = 55 - shower.thetaIndex();
 
        myBhaEvt->setPositron()->setThetaIndex( newthetaInd );
 
        unsigned int newphiInd = shower.phiIndex();
        myBhaEvt->setPositron()->setPhiIndex( newphiInd );
        myBhaEvt->setPositron()->setFound( true );
 
        log << MSG::INFO << name() << ": Positron: theta,phi energy "
            << myBhaEvt->positron()->theta() << "," << myBhaEvt->positron()->shower().phi()
            << " " << myBhaEvt->positron()->shower().energy() << endmsg;
      }
 
      // if ( (itrk==1&&theMdcTrk1->charge()<0)
      //   ||(itrk==2&&theMdcTrk2->charge()<0) ){
      if ( itrk == 2 )
      {
        // set properties including vertex corrected theta
        myBhaEvt->setElectron()->setShower( shower );
        myBhaEvt->setElectron()->setTheta( shower.theta() );
        myBhaEvt->setElectron()->setPhi( shower.phi() );
        unsigned int newthetaInd = -999;
        // The thetaIndex is defined by Endcap_east(0-5),Barrel(6-49),Endcap_west(50-55)
        // in Emc Bhabha Calibration
        if ( shower.module() == 0 ) newthetaInd = shower.thetaIndex();
        if ( shower.module() == 1 ) newthetaInd = shower.thetaIndex() + 6;
        if ( shower.module() == 2 ) newthetaInd = 55 - shower.thetaIndex();
 
        myBhaEvt->setElectron()->setThetaIndex( newthetaInd );
        unsigned int newphiInd = shower.phiIndex();
        myBhaEvt->setElectron()->setPhiIndex( newphiInd );
        myBhaEvt->setElectron()->setFound( true );
 
        log << MSG::INFO << name() << ": Electron: theta,phi energy "
            << myBhaEvt->electron()->theta() << "," << myBhaEvt->electron()->shower().phi()
            << " " << myBhaEvt->electron()->shower().energy() << endmsg;
      }
    }
  }
}

Referenced by execute().

finalize()

StatusCode EmcSelBhaEvent::finalize

(

)

Definition at line 423 of file EmcSelBhaEvent.cxx.

                                    {
 
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in finalize()" << endmsg;
 
  // output the data of Matrix and vector to files
  OutputMV();
  /*
  for (int i=1000; i < 1500; i++) {
        int xtalIndex=myCalibData->xtalInd(i);
 
        int nhitxtal=myCalibData->xtalHitsDir(xtalIndex);
        cout<<"ixtal, Nhitdir="<< xtalIndex<<"   "<<nhitxtal<<endl;
 
  }
  */
  if ( m_writeMVToFile )
  {
    delete myCalibData;
    myCalibData = 0;
  }
 
  cout << "Event=" << m_events << endl;
 
  cout << "m_Nothing =" << m_Nothing << endl;
  cout << "m_OneProng=" << m_OneProng << endl;
 
  cout << "m_TwoProngMatched=" << m_TwoProngMatched << endl;
 
  cout << "m_TwoProngOneMatched=" << m_TwoProngOneMatched << endl;
 
  cout << "m_showersAccepted=" << m_showersAccepted << endl;
 
  return StatusCode::SUCCESS;
}

findPhiDiff()

double EmcSelBhaEvent::findPhiDiff	(	double	phi1,
		double	phi2 )

Definition at line 1424 of file EmcSelBhaEvent.cxx.

                                                             {
  double diff;
  diff = phi1 - phi2; // rad
 
  while ( diff > pi ) diff -= twopi;
  while ( diff < -pi ) diff += twopi;
 
  return diff;
}

index()

int EmcSelBhaEvent::index ( int theta, int phi ) const

inline

Definition at line 83 of file EmcSelBhaEvent.h.

                                        {
    int val = ( ( m_index )[theta][phi] );
    return ( val );
  }

Referenced by CollectBhabha(), and execute().

initGeom()

void EmcSelBhaEvent::initGeom

(

)

Definition at line 1307 of file EmcSelBhaEvent.cxx.

                              {
 
  MsgStream log( msgSvc(), name() );
 
  int numberOfXtalsRing;
 
  EmcStructure* theEmcStruc = new EmcStructure();
  theEmcStruc->setEmcStruc();
 
  // number Of Theta Rings
  int nrOfTheRings = theEmcStruc->getNumberOfTheRings();
 
  m_nXtals = theEmcStruc->getNumberOfXtals();
 
  // init geometry
  for ( int the = theEmcStruc->startingTheta(); the < nrOfTheRings; the++ )
  {
 
    numberOfXtalsRing = theEmcStruc->crystalsInRing( (unsigned int)the );
 
    for ( int phi = 0; phi < numberOfXtalsRing; phi++ )
    { m_index[the][phi] = theEmcStruc->getIndex( (unsigned int)the, (unsigned int)phi ); }
  }
 
  log << MSG::INFO << " Emc geometry for Bhabha calibration initialized !! " << endl
      << "Number of Xtals: " << m_nXtals << endmsg;
  delete theEmcStruc;
}

Referenced by initialize().

initialize()

StatusCode EmcSelBhaEvent::initialize

(

)

Definition at line 209 of file EmcSelBhaEvent.cxx.

                                      {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  //---------------------------------------<<<<<<<<<<
  m_showersAccepted = 0;
  m_events          = 0;
  m_Nothing         = 0;
  m_OneProng        = 0;
  // number of events with   TwoProngMatched
  m_TwoProngMatched = 0;
  // number of events with   TwoProngOneMatched
  m_TwoProngOneMatched = 0;
 
  //--------------------------------------->>>>>>>>>>>
  // initialize Emc geometry to convert between index <=> theta,phi
  initGeom();
 
  // create the object that holds the calibration data
  if ( m_writeMVToFile ) myCalibData = new EmcBhaCalibData( m_nXtals, m_MsgFlag );
  else myCalibData = 0;
 
  // get EmcRecGeoSvc
 
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode   sc         = svcLocator->service( "EmcRecGeoSvc", m_iGeoSvc );
  if ( sc != StatusCode::SUCCESS ) { cout << "Error: Can't get EmcRecGeoSvc" << endl; }
 
  /*
  // use EmcCalibConstSvc
  StatusCode scCalib;
  scCalib = Gaudi::svcLocator() -> service("EmcCalibConstSvc", m_emcCalibConstSvc);
  if( scCalib != StatusCode::SUCCESS){
    log << MSG::ERROR << "can not use EmcCalibConstSvc" << endmsg;
  }
  else {
    std::cout << "Test EmcCalibConstSvc   DigiCalibConst(0)=  "
              << m_emcCalibConstSvc -> getDigiCalibConst(0) << std::endl;
  }
  */
 
  StatusCode scBeamEnergy;
  scBeamEnergy = Gaudi::svcLocator()->service( "BeamEnergySvc", m_BeamEnergySvc );
 
  if ( scBeamEnergy != StatusCode::SUCCESS )
  { log << MSG::ERROR << "can not use BeamEnergySvc" << endmsg; }
 
  // read correction function from the file of 'cor.conf'
  // from MC Bhabha data,
  // expected depostion energy for bhabha calibration at cms. system
  // it is as a function of thetaID(0-55)
  readCorFun();
  // read Esigma(Itheta)
  // from MC Bhabha data,
  // it is as a function of thetaID(0-55)  from the file of 'sigma.conf'
  readEsigma();
 
  // read peak of bhabha rawdata before bhabha calibration,
  // it is as a function of thetaID(0-55) from the file of "peakI.conf"
  readDepEneFun();
 
  // read sigma of bhabha rawdata before bhabha calibration,
  // it is as a function of thetaID(0-55) from the file of "sigmaI.conf"
  readSigmaExp();
  readRawPeakIxtal();
  // read data/mc(expected) with ixtal
  readDataAndMCIxtal();
  /*
  ofstream out;
  out.open("expectedE.conf");
  for(int i=0; i<6240;i++){
    out<<i<<"\t"<<expectedEnergy(i)<<endl;
  }
  out.close();
  */
  return StatusCode::SUCCESS;
}

OutputMV()

void EmcSelBhaEvent::OutputMV

(

)

Definition at line 1336 of file EmcSelBhaEvent.cxx.

                              {
 
  MsgStream log( msgSvc(), name() );
 
  // check this first because I sometimes got two endJob transitions
  if ( myCalibData != 0 )
 
    // if set write the matrix and vector to a file
    if ( m_writeMVToFile )
    {
 
      int  count;
      char cnum[20];
      if ( m_run < 0 ) { count = sprintf( cnum, "MC%d_%i", -m_run, m_Num ); }
      if ( m_run >= 0 ) { count = sprintf( cnum, "%d_%i", m_run, m_fileNum ); }
 
      std::string runnum = "";
      runnum.append( cnum );
 
      ofstream    runnumberOut;
      std::string runnumberFile = m_fileDir;
      runnumberFile += m_fileExt;
      runnumberFile += "runnumbers.dat";
      runnumberOut.open( runnumberFile.c_str(), ios::out | ios::app );
 
      ifstream runnumberIn;
      runnumberIn.open( runnumberFile.c_str() );
      bool status = false;
      while ( !runnumberIn.eof() )
      {
 
        std::string num;
        runnumberIn >> num;
        if ( runnum == num )
        {
          status = true;
          log << MSG::INFO << " the runnumber: " << runnum
              << " exists in the file runnumbers.dat" << endmsg;
          break;
        }
        else
        {
          status = false;
          log << MSG::INFO << " the runnumber: " << runnum
              << " does not exist in the file runnumbers.dat" << endmsg;
        }
      }
      runnumberIn.close();
 
      ofstream    vectorOut;
      std::string vectorFile = m_fileDir;
      vectorFile += m_fileExt;
      vectorFile += runnum;
      vectorFile += "CalibVector.dat";
      vectorOut.open( vectorFile.c_str() );
 
      ofstream    matrixOut;
      std::string matrixFile = m_fileDir;
      matrixFile += m_fileExt;
      matrixFile += runnum;
      matrixFile += "CalibMatrix.dat";
      matrixOut.open( matrixFile.c_str() );
 
      if ( vectorOut.good() && matrixOut.good() && runnumberOut.good() )
      {
 
        myCalibData->writeOut( matrixOut, vectorOut );
 
        log << MSG::INFO << " Wrote files " << ( vectorFile ) << " and " << ( matrixFile )
            << endmsg;
 
        if ( !status )
        {
          runnumberOut << runnum << "\n";
          log << MSG::INFO << "Wrote files " << runnumberFile << endmsg;
        }
      }
      else
        log << MSG::WARNING << " Could not open vector and/or matrix file !" << endl
            << "matrix file : " << matrixFile << endl
            << "vector file : " << vectorFile << endmsg;
 
      vectorOut.close();
      matrixOut.close();
      runnumberOut.close();
    }
}

Referenced by finalize().

passed()

bool EmcSelBhaEvent::passed ( )

inline

Definition at line 72 of file EmcSelBhaEvent.h.

{ return m_passed; }

Referenced by setPassed().

readCorFun()

void EmcSelBhaEvent::readCorFun

(

)

Definition at line 1434 of file EmcSelBhaEvent.cxx.

                                {
 
  ifstream    corFunIn;
  std::string corFunFile = m_inputFileDir;
  corFunFile += m_fileExt;
  corFunFile += "cor.conf";
  corFunIn.open( corFunFile.c_str() );
  for ( int i = 0; i < 56; i++ )
  {
    corFunIn >> m_corFun[i];
 
    cout << "energy corFun=" << m_corFun[i] << endl;
  }
  corFunIn.close();
}

Referenced by initialize().

readDataAndMCIxtal()

void EmcSelBhaEvent::readDataAndMCIxtal

(

)

Definition at line 1531 of file EmcSelBhaEvent.cxx.

                                        {
  // read mean of e5x5cms
  ifstream    EIn;
  std::string EFile = m_inputFileDir;
  EFile += m_fileExt;
  EFile += "X4680-meaNo3-e5x5cmsCut0.7-2mc.conf";
  EIn.open( EFile.c_str() );
 
  double MeanData[6240], MeanMC[6240];
  double RmsData[6240], RmsMC[6240];
  int    ixtal;
  for ( int i = 0; i < 6240; i++ )
  {
    EIn >> ixtal >> MeanData[i] >> MeanMC[i] >> RmsData[i] >> RmsMC[i];
    m_eMeanData[ixtal] = MeanData[i];
    m_eMeanMC[ixtal]   = MeanMC[i];
    // m_eRmsData[ixtal]=RmsData[i];
    // m_eRmsMC[ixtal]=RmsMC[i];
  }
  EIn.close();
 
  /////////////////// read peak of e5x5cms
 
  ifstream EIn1, EIn2, EIn3, EIn4;
 
  std::string EFile1 = m_inputFileDir;
  std::string EFile2 = m_inputFileDir;
  std::string EFile3 = m_inputFileDir;
  std::string EFile4 = m_inputFileDir;
 
  EFile1 += "findpeak-all-IV.conf";
  EFile2 += "findrms-all-IV.conf";
  EFile3 += "peakIxtalMC.conf";
  EFile4 += "findrms-all-MC.conf";
 
  EIn1.open( EFile1.c_str() );
  EIn2.open( EFile2.c_str() );
  EIn3.open( EFile3.c_str() );
  EIn4.open( EFile4.c_str() );
 
  // double MeanData[6240],MeanMC[6240];
  // double RmsData[6240],RmsMC[6240];
  int ixtal1, ixtal2, ixtal3, ixtal4;
  for ( int i = 0; i < 6240; i++ )
  {
    EIn1 >> ixtal1 >> MeanData[i];
    EIn2 >> ixtal2 >> RmsData[i];
    EIn3 >> ixtal3 >> MeanMC[i];
    EIn4 >> ixtal4 >> RmsMC[i];
    // m_eMeanData[ixtal1]=MeanData[i];
    m_eRmsData[ixtal2] = RmsData[i];
    // m_eMeanMC[ixtal3]=MeanMC[i];
    m_eRmsMC[ixtal4] = RmsMC[i];
    // cout<<"..........."<<
    // m_eMeanData[ixtal1]<<"\t"<<m_eRmsData[ixtal2]<<"\t"<<m_eMeanMC[ixtal3]<<"\t"<<m_eRmsMC[ixtal4]<<endl;
  }
 
  EIn1.close();
  EIn2.close();
  EIn3.close();
  EIn4.close();
}

Referenced by initialize().

readDepEneFun()

void EmcSelBhaEvent::readDepEneFun

(

)

Definition at line 1466 of file EmcSelBhaEvent.cxx.

                                   {
 
  ifstream    EDepEneIn;
  std::string EDepEneFile = m_inputFileDir;
  EDepEneFile += m_fileExt;
  EDepEneFile += "peakI.conf";
  EDepEneIn.open( EDepEneFile.c_str() );
  for ( int i = 0; i < 56; i++ )
  {
    EDepEneIn >> m_eDepEne[i];
    // m_eDepEne[i]=m_eDepEne[i]-0.020;
    // m_eDepEne[i]=m_eDepEne[i]/1.843*m_beamEnergy;
    cout << "DepEne=" << m_eDepEne[i] << endl;
  }
  EDepEneIn.close();
}

Referenced by initialize().

readEsigma()

void EmcSelBhaEvent::readEsigma

(

)

Definition at line 1450 of file EmcSelBhaEvent.cxx.

                                {
 
  ifstream    EsigmaIn;
  std::string EsigmaFile = m_inputFileDir;
  EsigmaFile += m_fileExt;
  EsigmaFile += "sigma.conf";
  EsigmaIn.open( EsigmaFile.c_str() );
  for ( int i = 0; i < 56; i++ )
  {
    EsigmaIn >> m_eSigma[i];
 
    cout << "Sigma=" << m_eSigma[i] << endl;
  }
  EsigmaIn.close();
}

Referenced by initialize().

readRawPeakIxtal()

void EmcSelBhaEvent::readRawPeakIxtal

(

)

Definition at line 1497 of file EmcSelBhaEvent.cxx.

                                      {
 
  ifstream    EIn;
  std::string EFile = m_inputFileDir;
  EFile += m_fileExt;
  EFile += "findpeak.conf";
  EIn.open( EFile.c_str() );
 
  double Peak[6240];
  int    ixtal;
  for ( int i = 0; i < 6240; i++ )
  {
    EIn >> ixtal >> Peak[i];
    m_eRawPeak[ixtal] = Peak[i];
  }
  EIn.close();
 
  /*
  ifstream EIn1;
  std::string EFile1 = m_inputFileDir;
  EFile1 += m_fileExt;
  EFile1 += "findpeak-mc.conf";
  EIn1.open(EFile1.c_str());
 
  double Peak1[6240];
  int ixtal1;
  for(int i=0;i<6240;i++) {
    EIn1>>ixtal1>>Peak1[i];
  m_eMcPeak[ixtal1]=Peak1[i];
  }
  EIn1.close();
  */
}

Referenced by initialize().

readSigmaExp()

void EmcSelBhaEvent::readSigmaExp

(

)

Definition at line 1483 of file EmcSelBhaEvent.cxx.

                                  {
  ifstream    ESigmaExpIn;
  std::string ESigmaExpFile = m_inputFileDir;
  ESigmaExpFile += m_fileExt;
  ESigmaExpFile += "sigmaI.conf";
  ESigmaExpIn.open( ESigmaExpFile.c_str() );
  for ( int i = 0; i < 56; i++ )
  {
    ESigmaExpIn >> m_eSigmaExp[i];
    cout << "SigmaExp=" << m_eSigmaExp[i] << endl;
  }
  ESigmaExpIn.close();
}

Referenced by initialize().

SelectBhabha()

void EmcSelBhaEvent::SelectBhabha

(

)

Definition at line 500 of file EmcSelBhaEvent.cxx.

                                  {
  MsgStream log( msgSvc(), name() );
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  m_selectedType = BhabhaType::Nothing;
  m_events++;
 
  /*
  Vint iGood;
  iGood.clear();
  int nCharge = 0;
  int numberOfTrack = 0;
 
  //the accumulated event information
  double totalEnergy = 0.;
  int numberOfShowers = 0;
  int numberOfShowersWithTrack = 0;
  bool secondUnMTrackFound=false;
  float eNorm = 0.;
  float pNorm = 0.;
  float acolli_min = 999.;
  double minAngFirstSh = 999., minAngSecondSh = 999.;
  float theFirstTracksTheta = 0;
  float theFirstTracksMomentum = 0;
 
  //the selection candidates
  RecMdcTrack *theFirstTrack = 0;
  RecMdcTrack *theSecondTrack = 0;
  RecMdcTrack *theKeptTrack = 0;
  RecEmcShower *theFirstShower = 0;
 
  RecEmcShower *theSecondShower = 0;
  RecEmcShower *theKeptShower = 0;
  double keptTheta = 0.;
  int theFirstShID = 9999;
  int theSecondShID = 9999;
  int theKeptShID = 9999;
  int theTrkID = 9999;
  */
 
  Vint iGam;
  iGam.clear();
 
  // double ene5x5,theta,phi,eseed;
  // double showerX,showerY,showerZ;
  // long int thetaIndex,phiIndex;
  // HepPoint3D showerPosition(0,0,0);
  //  RecEmcID showerId;
  // unsigned int npart;
 
  for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
  {
 
    if ( i >= evtRecTrkCol->size() ) break;
 
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
 
    if ( ( *itTrk )->isEmcShowerValid() )
    {
 
      RecEmcShower* theShower = ( *itTrk )->emcShower();
      int           TrkID     = ( *itTrk )->trackId();
      double        Shower5x5 = theShower->e5x5();
      // cout<<"Shower5x5="<<Shower5x5<<endl;
      /*
      ene5x5=theShower->e5x5();   //uncorrected 5x5 energy unit GeV
      eseed=theShower->eSeed(); //unit GeV
 
      showerPosition = theShower->position();
      theta = theShower->theta();
      phi= theShower->phi();
      showerX=theShower->x();
      showerY=theShower->y();
      showerZ=theShower->z();
 
      showerId = theShower->getShowerId();
 
      npart = EmcID::barrel_ec(showerId);
 
      //The theta number is defined by Endcap_east(0-5),Barrel(0-43),Endcap_west(0-5)
      //module is defined  by Endcap_east(0),Barrel(1),Endcap_west(2)
 
      thetaIndex=EmcID::theta_module(showerId);
 
      phiIndex=EmcID::phi_module(showerId);
      */
 
      if ( Shower5x5 > ( 0.6 * m_beamEnergy ) ) { iGam.push_back( TrkID ); }
    }
  }
  int nGam = iGam.size();
  // log << MSG::DEBUG << "num Good Photon " << nGam  << " , " <<evtRecEvent->totalNeutral()
  // <<endmsg;
 
  if ( nGam == 2 )
  {
 
    m_selectedType   = BhabhaType::TwoProngMatched;
    m_selectedTrkID1 = iGam[0];
    m_selectedTrkID2 = iGam[1];
  }
 
  //return StatusCode::SUCCESS;
}

Referenced by execute().

selectedTrkID1()

int EmcSelBhaEvent::selectedTrkID1 ( ) const

inline

Definition at line 77 of file EmcSelBhaEvent.h.

{ return m_selectedTrkID1; }

selectedTrkID2()

int EmcSelBhaEvent::selectedTrkID2 ( ) const

inline

Definition at line 79 of file EmcSelBhaEvent.h.

{ return m_selectedTrkID2; }

selectedType()

int EmcSelBhaEvent::selectedType ( ) const

inline

Definition at line 75 of file EmcSelBhaEvent.h.

{ return m_selectedType; }

SelectFillBhabha()

StatusCode EmcSelBhaEvent::SelectFillBhabha

(

)

setPassed()

void EmcSelBhaEvent::setPassed ( bool passed )

inline

Definition at line 73 of file EmcSelBhaEvent.h.

{ m_passed = passed; }

---

The documentation for this class was generated from the following files:

• EmcSelBhaEvent.h
• EmcSelBhaEvent.cxx