#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 16 of file FarmMonitorAlg.h.

Constructor & Destructor Documentation

◆ FarmMonitorAlg()

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

Definition at line 37 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);
 
}

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 173 of file FarmMonitorAlg.cxx.

                                   {
  
  double eBeam = m_ecm/2;
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
 
  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() <<endreq;
 
  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" << endreq;
    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 526 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