---

EventHeaderMaker.cxx

Go to the documentation of this file.

// ================================================
// EventHeaderMaker class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast
//

#include "AtlfastAlgs/EventHeaderMaker.h"
#include "AtlfastAlgs/MissingMomentum.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "AtlfastEvent/EventHeader.h"
#include "AtlfastEvent/IKinematic.h"
#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/Cell.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/MCparticleCollection.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/Jet.h"

#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/TesIO.h"

// standard library includes
#include <cmath> 
#include <string>
#include <vector>
#include <algorithm>
#include <assert.h>

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "StoreGate/DataHandle.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"


//Generator includes
#include "GeneratorObjects/McEventCollection.h"

namespace Atlfast {
  //--------------------------------
  // Constructors and destructors
  //-------------------------------
  EventHeaderMaker::EventHeaderMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){
    
    // Default paths for entities in the TES
    m_electronLocation         = "/Event/AtlfastIsolatedElectrons";
    m_photonLocation           = "/Event/AtlfastIsolatedPhotons";
    m_isolatedMuonLocation     = "/Event/AtlfastIsolatedMuons";
    m_nonIsolatedMuonLocation  = "/Event/AtlfastNonIsolatedMuons";
    m_jetLocation              = "/Event/AtlfastJets";
    m_cellLocation             = "/Event/AtlfastCells";
    m_clusterLocation          = "/Event/AtlfastClusters";
    m_mcTruthLocation          = "/Event/McEventCollection";
    m_outputLocation           = "/Event/AtlfastEventHeader";
    m_missingMomentumLocation  = "/Event/AtlfastMissingMomentum";
    m_beamEnergy               = 7000.0;
    m_testMode                 = false;
    
    
    
    declareProperty( "MissingMomentumLocation", m_missingMomentumLocation ) ;
    declareProperty( "ElectronLocation",        m_electronLocation ) ;
    declareProperty( "PhotonLocation",          m_photonLocation );
    declareProperty( "IsolatedMuonLocation",    m_isolatedMuonLocation );
    declareProperty( "NonIsolatedMuonLocation", m_nonIsolatedMuonLocation );
    declareProperty( "JetLocation",             m_jetLocation );
    declareProperty( "McTruthLocation",         m_mcTruthLocation );
    declareProperty( "OutputLocation",          m_outputLocation );
    declareProperty( "BeamEnergy",              m_beamEnergy );
    
    declareProperty( "TestMode",                m_testMode ) ;
    
    
    
  }
  
  
  //-------------------------------------
  // Destructor 
  //-------------------------------------
  EventHeaderMaker::~EventHeaderMaker(){
    delete m_tesIO;
  }

  //---------------------------------
  // initialise() 
  //---------------------------------
  StatusCode EventHeaderMaker::initialize()
  {
    // .. put any initialisation in here...
    MsgStream log(messageService(), name());
    log << MSG::DEBUG << "in initialize()" << endreq;
    
    //set up the list of particles which are considered as "escaping"
    m_escapingParticles.push_back(ParticleCodes::NU_E);
    m_escapingParticles.push_back(ParticleCodes::NU_MU);
    m_escapingParticles.push_back(ParticleCodes::NU_TAU);

    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO();

    //get the Global Event Data from singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    m_visibleToAtlas   = ged -> visibleToAtlas();
    m_lumi             = ged -> lumi();
    m_barrelForwardEta = ged ->  barrelForwardEta();
    HeaderPrinter hp("Atlfast EventHeader Maker:", log);
    hp.add("TestMode                    ", m_testMode);
    hp.add("EndCap Begins at (eta)      ", m_barrelForwardEta);
    hp.add("Luminosity                  ", m_lumi);
    hp.add("Beam Energy                 ", m_beamEnergy);
    hp.add("MC location                 ", m_mcTruthLocation);
    hp.add("Electron Location           ", m_electronLocation);
    hp.add("Photon   Location           ", m_photonLocation);
    hp.add("Isolated Muon Location      ", m_isolatedMuonLocation);
    hp.add("Non-Isolated Muon Location  ", m_nonIsolatedMuonLocation);
    hp.add("Jet Location                ", m_jetLocation);
    hp.add("Cell Location               ", m_cellLocation);
    hp.add("Cluster Location            ", m_clusterLocation);
    hp.add("Unused cell+cluster sum     ", m_missingMomentumLocation);
    hp.add("Output Location             ", m_outputLocation);    
    hp.print();
    

    
    
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  StatusCode EventHeaderMaker::finalize()
  {
    // .. put any finalisation in here...
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing " << endreq;
    
    return StatusCode::SUCCESS ;
  }
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  StatusCode EventHeaderMaker::execute( )
  {
    
    //................................
    // make a message logging stream
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Executing " << endreq;
    std::string mess;
    
      

    //
    // Check each particle list in turn and get the number of entries
    //
    int nElectrons         = 
      numberInList<ReconstructedParticleCollection>(m_electronLocation);
    int nPhotons           = 
      numberInList<ReconstructedParticleCollection>(m_photonLocation);
    int nIsolatedMuons     = 
      numberInList<ReconstructedParticleCollection>(m_isolatedMuonLocation);
    int nNonIsolatedMuons  = 
      numberInList<ReconstructedParticleCollection>(m_nonIsolatedMuonLocation);
    int nJets              = numberInList<JetCollection>(m_jetLocation);
    int nBJets             = numberJetFlavor(m_jetLocation,
                                             ParticleCodes::BQUARK);
    int nCJets             = numberJetFlavor(m_jetLocation,
                                             ParticleCodes::CQUARK);
    int nTauJets           = numberJetFlavor(m_jetLocation,
                                             ParticleCodes::TAU);
    nTauJets              += numberJetFlavor(m_jetLocation,
                                             -ParticleCodes::TAU);
    
    log << MSG::DEBUG << "Electrons: " << nElectrons << endreq;
    log << MSG::DEBUG << "Photons: " << nPhotons << endreq;
    log << MSG::DEBUG << "IsolatedMuons: " << nIsolatedMuons << endreq;
    log << MSG::DEBUG << "NonIsolatedMuons: " << nNonIsolatedMuons << endreq;
    log << MSG::DEBUG << "Jets: " << nJets << endreq;
    
    
    // Event Shape variable calculations
    
    double jetCircularity = 0.0; //FIXME not written yet
    double eventCircularity = 0.0; //FIXME not written yet
    double thrust = 0.0; //FIXME not written yet
    double oblateness = 0.0; //FIXME not written yet
    
    //missing momentum
    HepLorentzVector pMiss = missingMomentum(log);
    
    log << MSG::DEBUG << "Missing momentum " <<pMiss<<endreq; 
    
    HepLorentzVector pEscaped = escapedMomentum(log);
    
    log << MSG::DEBUG << "Escaped momentum " << pEscaped<<endreq;

    EventHeader* header = 
      new EventHeader(
                      nElectrons, 
                      nIsolatedMuons, 
                      nNonIsolatedMuons, 
                      nPhotons, 
                      nJets,
                      nBJets, 
                      nCJets, 
                      nTauJets, 
                      jetCircularity, 
                      eventCircularity,
                      thrust, 
                      oblateness, 
                      pMiss, 
                      pEscaped);
    
    TesIoStat stat = m_tesIO->store(header,m_outputLocation) ;
    std::string message = stat? "Stored Evt Header":"Failed to Evt Header";
    log<<MSG::DEBUG<<message<<endreq;
     return stat;

  }
  //-------------------------------------------------------
  // missingMomentum
  //-------------------------------------------------------
  HepLorentzVector EventHeaderMaker::missingMomentum(MsgStream& log) {

    HepLorentzVector unused(0., 0., 0., 0.);
    const DataHandle<MissingMomentum> mm;
    TesIoStat stat = m_tesIO->getDH(mm,m_missingMomentumLocation );
    if(!stat){
      log << MSG::WARNING<<"Could Not Find Missing momentum in TES"<<endreq;
      return unused;
    }
    log << MSG::DEBUG<<"Found Missing Momentum in TES"<<endreq;

    unused = mm->momentum();

    HepLorentzVector temp(0.0,0.0,0.0,0.0);
    

    temp = unused
      +  momentumSum<ReconstructedParticleCollection> (m_electronLocation)
      +  momentumSum<ReconstructedParticleCollection> (m_photonLocation)
      +  momentumSum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
      +  momentumSum<JetCollection>                   (m_jetLocation)
      +  momentumSum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
      -  assocMomentumSum<JetCollection, ReconstructedParticle> (m_jetLocation)
      ;

    log<<MSG::DEBUG
       <<"Total Imbalanced momentum     "
       <<temp
       <<endreq;



    if(m_testMode){
        HepLorentzVector totCluster = 
          momentumSum<ClusterCollection> (m_clusterLocation);

        HepLorentzVector clusterAssE = 
          assocMomentumSum<ReconstructedParticleCollection, Cluster>
          (m_electronLocation);

        HepLorentzVector clusterAssP = 
          assocMomentumSum<ReconstructedParticleCollection, Cluster>
          (m_photonLocation);

        HepLorentzVector clusterAssM = 
          assocMomentumSum<ReconstructedParticleCollection, Cluster>
          (m_isolatedMuonLocation);
        
        HepLorentzVector clusterAssJ = 
          assocMomentumSum<JetCollection, Cluster>(m_jetLocation);

      HepLorentzVector clusterDiff = totCluster-
        clusterAssJ-
        clusterAssE-
        clusterAssM-
        clusterAssP;

      HepLorentzVector totCell = momentumSum<CellCollection>(m_cellLocation);

      HepLorentzVector cellAssCluster = 
        assocMomentumSum<ClusterCollection, Cell>(m_clusterLocation);

      HepLorentzVector cellDiff = totCell-cellAssCluster;

      HepLorentzVector diffS  = clusterDiff+cellDiff-unused;
      

      
      log<<MSG::DEBUG<<"Total cluster    "<<totCluster<<endreq;
      log<<MSG::DEBUG<<"Clus Ass Jet     "<<clusterAssJ<<endreq;
      log<<MSG::DEBUG<<"Clus Ass El      "<<clusterAssE<<endreq;
      log<<MSG::DEBUG<<"Clus Ass Mu      "<<clusterAssM<<endreq;
      log<<MSG::DEBUG<<"Clus Ass Ph      "<<clusterAssP<<endreq;
      log<<MSG::DEBUG<<"Clus diff        "<<clusterDiff<<endreq;
      log<<MSG::DEBUG<<"Total Cell       "<<totCell<<endreq;
      log<<MSG::DEBUG<<"Cell Ass Clus    "<<cellAssCluster<<endreq;
      log<<MSG::DEBUG<<"Clus diff        "<<cellDiff<<endreq;
      log<<MSG::DEBUG<<"Cell+Clus diff   "<<cellDiff<<endreq;
      log<<MSG::DEBUG<<"Unused           "<<unused<<endreq;
      log<<MSG::DEBUG<<"diff             "<<diffS<<endreq;

      log<<MSG::DEBUG<<endreq;
      
      assert( diffS.rho() <0.001);





      HepLorentzVector temp2 =  
        momentumSum<ReconstructedParticleCollection>    (m_electronLocation)
        +  momentumSum<ReconstructedParticleCollection> (m_photonLocation)
        +  momentumSum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
        +  momentumSum<JetCollection>  (m_jetLocation)
        +  momentumSum<ClusterCollection> (m_clusterLocation)
        +  momentumSum<CellCollection> (m_cellLocation)
        +  momentumSum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_electronLocation)   
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_photonLocation)
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_nonIsolatedMuonLocation)
        
        -  assocMomentumSum<JetCollection, Cluster>
        (m_jetLocation)
        
        -  assocMomentumSum<JetCollection, ReconstructedParticle>
        (m_jetLocation)

        - assocMomentumSum<ClusterCollection, Cell>
        (m_clusterLocation);

      
      assert( (temp-temp2).rho() <0.001);
    }
    
    // calculate missing momentum from the reconstructed sum of momenta:
    //
    temp.setPx( -temp.px() );
    temp.setPy( -temp.py() );
    temp.setPz( -temp.pz() );
    temp.setE(  2.0*m_beamEnergy - temp.e() );
    
    return temp;
  }
  
  
  
  //-------------------------------------------------------
  // escapedMomentum
  //-------------------------------------------------------
  HepLorentzVector EventHeaderMaker::escapedMomentum(MsgStream& log) {

    HepLorentzVector temp(0.0,0.0,0.0,0.0);

    MCparticleCollection  mcParticles;
    TesIoStat stat = m_tesIO->getMC( mcParticles, m_visibleToAtlas ) ;
    std::string mess = 
      stat? "MC particles retrieved":"MC particles retrieve failed";
    log<<MSG::DEBUG<<mess<<endreq;
    

    MCparticleCollection::const_iterator part = mcParticles.begin();
    for (;part != mcParticles.end(); ++part) {
      temp += (*part)->momentum(); 
    }

    return temp;
  }

} // end of namespace bracket

---