#include <FarmMonitorAlg.h>

Inheritance diagram for FarmMonitorAlg:

Public Member Functions
	FarmMonitorAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()

Detailed Description

Definition at line 15 of file FarmMonitorAlg.h.

Constructor & Destructor Documentation

◆ FarmMonitorAlg()

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

Definition at line 35 of file FarmMonitorAlg.cxx.

    : Algorithm( name, pSvcLocator ) {
 
  // Declare the properties
  declareProperty( "PrintRunEventFreq", m_RunEventFreq = 10 );
 
  declareProperty( "Ecm", m_ecm = 3.686 );
  declareProperty( "Vr0cut", m_vr0cut = 1.0 );
  declareProperty( "Vz0cut", m_vz0cut = 5.0 );
}

Referenced by FarmMonitorAlg().

Member Function Documentation

◆ execute()

StatusCode FarmMonitorAlg::execute ( )

Variables for charged track FOR loop

Charged track FOR loop

EMC energy associated with track

MUC information

Tracks after IP cuts

dE/dx information

TOF information

For the 2 highest momentum charged tracks

If the event has 1 IP charged track

If the event has 2 IP charged tracks

Shower (aka Neutral track) FOR loop

Energy of most energetic showers

Histograms filled per event

Total Energy Histograms

Total Charged Momentum Histograms

Total Neutral Momentum Histograms

Total Momentum Histograms

VeeVertex information

Loop over VeeVertex candidates

Definition at line 197 of file FarmMonitorAlg.cxx.

                                   {
 
  double eBeam = m_ecm / 2;
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  int                              run   = eventHeader->runNumber();
  int                              event = eventHeader->eventNumber();
  if ( event % m_RunEventFreq == 0 )
    std::cout << "Run " << run << ", event " << event << std::endl;
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  ////////////////////////////////////////
  /// Variables for charged track FOR loop
  int    nChargedTracks = 0, nChargedTracksIP = 0;
  int    nCharge = 0, nChargeIP = 0;
  double totalVisibleEnergy = 0, totalChargedEnergy = 0, totalEMCEnergy = 0;
  double totalChargedPX = 0, totalChargedPY = 0, totalChargedPZ = 0;
 
  double highestIPTrackP = -1, secondHighestIPTrackP = -2;
  int    highestPIPTrackId = -1, secondHighestPIPTrackId = -1;
 
  /// Charged track FOR loop
  for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isMdcTrackValid() ) continue;
    RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
 
    int    trackId = mdcTrk->trackId();
    int    charge  = mdcTrk->charge();
    double r0      = mdcTrk->r();
    double z0      = mdcTrk->z();
    h_trackR0->fill( r0 );
    h_trackZ0->fill( z0 );
 
    nChargedTracks++;
    nCharge += charge;
 
    double pX       = mdcTrk->px();
    double pY       = mdcTrk->py();
    double pZ       = mdcTrk->pz();
    double pMag     = mdcTrk->p();
    double cosTheta = cos( mdcTrk->theta() );
    double phi      = mdcTrk->phi();
 
    double chargedEnergy = sqrt( pMag * pMag + mPi * mPi );
    totalVisibleEnergy += chargedEnergy;
    totalChargedEnergy += chargedEnergy;
 
    totalChargedPX += pX;
    totalChargedPY += pY;
    totalChargedPZ += pZ;
 
    /// EMC energy associated with track
    if ( ( *itTrk )->isEmcShowerValid() )
    {
      RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
      totalEMCEnergy += emcChargedTrk->energy();
    }
 
    /// MUC information
    if ( ( *itTrk )->isMucTrackValid() )
    {
      RecMucTrack* mucTrk   = ( *itTrk )->mucTrack();
      double       mucDepth = mucTrk->depth();
      if ( mucDepth > 0 )
      {
        h_mucDepth->fill( mucDepth );
        h_mucDepthVsCosTheta->fill( cosTheta, mucDepth );
        h_mucDepthVsPhi->fill( phi, mucDepth );
      }
    } // End of "isMucShowerValid()" IF
 
    ////////////////////////
    /// Tracks after IP cuts
    if ( fabs( z0 ) >= m_vz0cut ) continue;
    if ( r0 >= m_vr0cut ) continue;
 
    nChargedTracksIP++;
    nChargeIP += charge;
 
    /// dE/dx information
    double dedxProbPH       = -10;
    int    dedxNumTotalHits = -10;
    int    dedxNumGoodHits  = -10;
 
    if ( ( *itTrk )->isMdcDedxValid() )
    {
      RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
 
      dedxNumTotalHits = dedxTrk->numTotalHits();
      dedxNumGoodHits  = dedxTrk->numGoodHits();
      h_dedxTotalHitsIP->fill( dedxNumTotalHits );
      h_dedxGoodHitsIP->fill( dedxNumGoodHits );
 
      h_dedxElecChiIP->fill( dedxTrk->chiE() );
      h_dedxMuonChiIP->fill( dedxTrk->chiMu() );
      h_dedxPionChiIP->fill( dedxTrk->chiPi() );
      h_dedxKaonChiIP->fill( dedxTrk->chiK() );
      h_dedxProtonChiIP->fill( dedxTrk->chiP() );
 
      dedxProbPH = dedxTrk->probPH();
      h_dedxProbPHIP->fill( dedxProbPH );
      h_dedxProbPHVsMomentumIP->fill( pMag, dedxProbPH );
 
    } // End of "isMdcDedxValid()" IF
 
    /// TOF information
    if ( ( *itTrk )->isTofTrackValid() )
    {
      SmartRefVector<RecTofTrack>           tofTrkCol = ( *itTrk )->tofTrack();
      SmartRefVector<RecTofTrack>::iterator iter_tof  = tofTrkCol.begin();
 
      for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
      {
        TofHitStatus* status = new TofHitStatus;
        status->setStatus( ( *iter_tof )->status() );
 
        if ( status->is_barrel() )
        {
          if ( !( status->is_counter() ) ) continue; // ?
 
          double tofPH = ( *iter_tof )->ph();
          double tof   = ( *iter_tof )->tof();
 
          h_tofElecIP_Barrel->fill( tof - ( *iter_tof )->texpElectron() );
          h_tofMuonIP_Barrel->fill( tof - ( *iter_tof )->texpMuon() );
          h_tofPionIP_Barrel->fill( tof - ( *iter_tof )->texpPion() );
          h_tofKaonIP_Barrel->fill( tof - ( *iter_tof )->texpKaon() );
          h_tofProtonIP_Barrel->fill( tof - ( *iter_tof )->texpProton() );
          h_tofVsMomentumIP->fill( pMag, tof );
 
          if ( status->layer() == 1 )
          {
            h_tofPHIP_BarrelLayer1->fill( tofPH );
            h_tofIP_BarrelLayer1->fill( tof );
          }
          if ( status->layer() == 2 )
          {
            h_tofPHIP_BarrelLayer2->fill( tofPH );
            h_tofIP_BarrelLayer2->fill( tof );
          }
        } // End of TOF barrel IF
 
        else
        {
          if ( !( status->is_counter() ) ) continue; // ?
 
          double tofPH = ( *iter_tof )->ph();
          double tof   = ( *iter_tof )->tof();
 
          h_tofElecIP_Endcap->fill( tof - ( *iter_tof )->texpElectron() );
          h_tofMuonIP_Endcap->fill( tof - ( *iter_tof )->texpMuon() );
          h_tofPionIP_Endcap->fill( tof - ( *iter_tof )->texpPion() );
          h_tofKaonIP_Endcap->fill( tof - ( *iter_tof )->texpKaon() );
          h_tofProtonIP_Endcap->fill( tof - ( *iter_tof )->texpProton() );
          h_tofVsMomentumIP->fill( pMag, tof );
 
          if ( status->is_east() )
          {
            h_tofPHIP_EastEndcap->fill( tofPH );
            h_tofIP_EastEndcap->fill( tof );
          }
          else
          {
            h_tofPHIP_WestEndcap->fill( tofPH );
            h_tofIP_WestEndcap->fill( tof );
          }
        } // End of TOF endcap IF
 
      } // End of "iter_tof" FOR
    }   // End of "isTofTrackValid()" IF
 
    /// For the 2 highest momentum charged tracks
    if ( pMag > highestIPTrackP )
    {
      secondHighestPIPTrackId = highestPIPTrackId;
      secondHighestIPTrackP   = highestIPTrackP;
      highestPIPTrackId       = trackId;
      highestIPTrackP         = pMag;
    }
    if ( ( pMag > secondHighestIPTrackP ) && ( pMag < highestIPTrackP ) )
    {
      secondHighestPIPTrackId = trackId;
      secondHighestIPTrackP   = pMag;
    }
 
  } // End of charged track FOR
 
  ///////////////////////////////////////
  /// If the event has 1 IP charged track
  if ( nChargedTracksIP > 0 )
  {
    EvtRecTrackIterator itTrk        = evtRecTrkCol->begin() + highestPIPTrackId;
    RecMdcTrack*        mdcTrk       = ( *itTrk )->mdcTrack();
    double              highestPPhi0 = mdcTrk->phi();
    h_pIPTrack1DivEb->fill( mdcTrk->p() / eBeam );
 
    if ( ( *itTrk )->isEmcShowerValid() )
    {
      RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
      h_eEMCIPTrack1DivEb->fill( emcChargedTrk->energy() / eBeam );
    }
 
    /// If the event has 2 IP charged tracks
    if ( nChargedTracksIP > 1 )
    {
      itTrk                           = evtRecTrkCol->begin() + secondHighestPIPTrackId;
      RecMdcTrack* mdcTrk             = ( *itTrk )->mdcTrack();
      double       secondHighestPPhi0 = mdcTrk->phi();
      h_pIPTrack2DivEb->fill( mdcTrk->p() / eBeam );
 
      if ( ( *itTrk )->isEmcShowerValid() )
      {
        RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
        h_eEMCIPTrack2DivEb->fill( emcChargedTrk->energy() / eBeam );
      }
 
      h_acoplanarity_2HighestPIPTracks->fill(
          fabs( CLHEP::pi - fabs( highestPPhi0 - secondHighestPPhi0 ) ) );
    } // End of "nChargedTracksIP > 1" IF
  }   // End of "nChargedTracksIP > 0" IF
 
  ///////////////////////////////////////
  /// Shower (aka Neutral track) FOR loop
  int    nNeutralTracks = 0, nNeutralTracksGT30MeV = 0;
  double totalNeutralEnergy = 0;
  double totalNeutralPX = 0, totalNeutralPY = 0, totalNeutralPZ = 0;
 
  for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
 
    nNeutralTracks++;
    double eraw = emcTrk->energy();
    if ( eraw > 0.030 ) nNeutralTracksGT30MeV++;
 
    totalVisibleEnergy += eraw;
    totalEMCEnergy += eraw;
    totalNeutralEnergy += eraw;
 
    double theta = emcTrk->theta();
    double phi   = emcTrk->phi();
 
    double pX = eraw * cos( phi ) * sin( theta );
    double pY = eraw * sin( phi ) * sin( theta );
    double pZ = eraw * cos( theta );
 
    totalNeutralPX += pX;
    totalNeutralPY += pY;
    totalNeutralPZ += pZ;
 
    /// Energy of most energetic showers
    if ( nNeutralTracks == 1 ) h_eEMCShower1DivEb->fill( eraw / eBeam );
    if ( nNeutralTracks == 2 ) h_eEMCShower2DivEb->fill( eraw / eBeam );
    if ( nNeutralTracks == 3 ) h_eEMCShower3DivEb->fill( eraw / eBeam );
 
  } // End of neutral track FOR
 
  ///////////////////////////////
  /// Histograms filled per event
 
  h_nChargedTracks->fill( nChargedTracks );
  h_nChargedTracksIP->fill( nChargedTracksIP );
 
  h_netCharge->fill( nCharge );
  h_netChargeIP->fill( nChargeIP );
 
  h_nNeutralTracks->fill( nNeutralTracks );
  h_nNeutralTracksGT30MeV->fill( nNeutralTracksGT30MeV );
 
  /// Total Energy Histograms
  h_eVisibleDivEcm->fill( totalVisibleEnergy / m_ecm );
  h_eEMCDivEcm->fill( totalEMCEnergy / m_ecm );
  h_eNeutralDivEcm->fill( totalNeutralEnergy / m_ecm );
  h_eChargedDivEcm->fill( totalChargedEnergy / m_ecm );
 
  /// Total Charged Momentum Histograms
  double totalChargedPXY =
      sqrt( totalChargedPX * totalChargedPX + totalChargedPY * totalChargedPY );
  double totalChargedP =
      sqrt( totalChargedPXY * totalChargedPXY + totalChargedPZ * totalChargedPZ );
 
  h_netMomentumDivEcm_AllChargedTracks->fill( totalChargedP / m_ecm );
  h_netTransMomentumDivEcm_AllChargedTracks->fill( totalChargedPXY / m_ecm );
  if ( totalChargedP > 0 )
  { h_cosTheta_AllChargedTracks->fill( totalChargedPZ / totalChargedP ); }
  else
  {
    if ( nChargedTracks > 0 )
    {
      log << MSG::INFO << "Run " << run << ", event " << event << ": totalChargedP <= 0! "
          << endmsg;
    }
  }
 
  /// Total Neutral Momentum Histograms
  double totalNeutralPXY =
      sqrt( totalNeutralPX * totalNeutralPX + totalNeutralPY * totalNeutralPY );
  double totalNeutralP =
      sqrt( totalNeutralPXY * totalNeutralPXY + totalNeutralPZ * totalNeutralPZ );
 
  h_netMomentumDivEcm_AllNeutralTracks->fill( totalNeutralP / m_ecm );
  h_netTransMomentumDivEcm_AllNeutralTracks->fill( totalNeutralPXY / m_ecm );
  if ( totalNeutralP > 0 )
  { h_cosTheta_AllNeutralTracks->fill( totalNeutralPZ / totalNeutralP ); }
  else
  {
    if ( nNeutralTracks > 0 )
    {
      log << MSG::INFO << "Run " << run << ", event " << event << ": totalNeutralP <= 0! "
          << endmsg;
    }
  }
 
  /// Total Momentum Histograms
  double totalEventPX = totalChargedPX + totalNeutralPX;
  double totalEventPY = totalChargedPY + totalNeutralPY;
  double totalEventPZ = totalChargedPZ + totalNeutralPZ;
 
  double totalEventPXY = sqrt( totalEventPX * totalEventPX + totalEventPY * totalEventPY );
  double totalEventP   = sqrt( totalEventPXY * totalEventPXY + totalEventPZ * totalEventPZ );
 
  h_netMomentumDivEcm_AllTracks->fill( totalEventP / m_ecm );
  h_netTransMomentumDivEcm_AllTracks->fill( totalEventPXY / m_ecm );
  if ( totalEventP > 0 ) { h_cosTheta_AllTracks->fill( totalEventPZ / totalEventP ); }
  else
  {
    log << MSG::INFO << "Run " << run << ", event " << event << ": totalEventP <= 0! "
        << endmsg;
  }
 
  /// VeeVertex information
  SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol( eventSvc(),
                                                       "/Event/EvtRec/EvtRecVeeVertexCol" );
  if ( !evtRecVeeVertexCol )
  {
    log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endmsg;
    return StatusCode::FAILURE;
  }
 
  /// Loop over VeeVertex candidates
  int numKs = 0, numLambda = 0;
  for ( EvtRecVeeVertexCol::iterator veeItr = evtRecVeeVertexCol->begin();
        veeItr < evtRecVeeVertexCol->end(); veeItr++ )
  {
 
    h_ksMass->fill( ( *veeItr )->mass() );
    if ( fabs( ( *veeItr )->mass() - mKs ) < 0.010 ) ++numKs;
 
    h_lambdaMass->fill( ( *veeItr )->mass() );
    if ( fabs( ( *veeItr )->mass() - mLambda ) < 0.010 ) ++numLambda;
 
  } // End of "evtRecVeeVertexCol" FOR
 
  h_nKs->fill( numKs );
  h_nLambda->fill( numLambda );
 
  return StatusCode::SUCCESS;
}

◆ finalize()

StatusCode FarmMonitorAlg::finalize ( )

Definition at line 569 of file FarmMonitorAlg.cxx.

                                    {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
 
  return StatusCode::SUCCESS;
}

◆ initialize()

StatusCode FarmMonitorAlg::initialize ( )

Total energies

Net momenta

Charged Tracks

2 highest momentum charged tracks

Neutral Tracks

MUC information

PID (dE/dx) information

PID (TOF) information