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/KtCluster.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 "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"


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


#include "CLHEP/Units/SystemOfUnits.h"

namespace Atlfast {

  //--------------------------------
  // Constructors and destructors
  //-------------------------------
  EventHeaderMaker::EventHeaderMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_tesIO(0),
      m_visibleToAtlas(0)
  {
    
    // 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_outputLocation           = "/Event/AtlfastEventHeader";
    m_missingMomentumLocation  = "/Event/AtlfastMissingMomentum";
    m_beamEnergy               = 7000.0*GeV;
    m_testMode                 = false;
    m_addHalos                 = false;
    m_adjustMissETForIsolation = true;
    
    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( "OutputLocation",          m_outputLocation );
    declareProperty( "BeamEnergy",              m_beamEnergy );
    declareProperty( "AddHalos",                m_addHalos );

    declareProperty( "TestMode",                m_testMode ) ;
    
    
    
  }
  
  
  //-------------------------------------
  // Destructor 
  //-------------------------------------
  EventHeaderMaker::~EventHeaderMaker(){
    if (m_tesIO) 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);


    //get the Global Event Data from singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    m_visibleToAtlas           = ged -> visibleToAtlas();
    m_lumi                     = ged -> lumi();
    m_barrelForwardEta         = ged -> barrelForwardEta();
    m_mcLocation               = ged -> mcLocation();
    m_adjustMissETForIsolation = ged->adjustMissETForIsolation();

    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());

    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_mcLocation);
    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
    
    // Cumulative ET sum
    double sumET = 0.0;

    //missing momentum and sumET
    HepLorentzVector pMiss = missingMomentum(log,sumET);
    
    log << MSG::DEBUG << "Missing momentum " <<pMiss<<endreq; 
    
    HepLorentzVector pEscaped = escapedMomentum(log);
    
    log << MSG::DEBUG << "Escaped momentum " << pEscaped<<endreq;

    //get a container of weight containers for this event
    MCweightContainerCollection weightCollections;
    TesIoStat stat = m_tesIO->getMCweightContainers(weightCollections);
    if(!stat){
      log << MSG::ERROR <<"Problem getting b-quarks" << endreq; 
      return stat;
    }

    EventHeader* header = 
      new EventHeader(
                      nElectrons, 
                      nIsolatedMuons, 
                      nNonIsolatedMuons, 
                      nPhotons, 
                      nJets,
                      nBJets, 
                      nCJets, 
                      nTauJets, 
                      jetCircularity, 
                      eventCircularity,
                      thrust, 
                      oblateness, 
                      pMiss,
                      sumET,
                      pEscaped,
                      weightCollections
                      );
    
    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.status();
    
  }
  
  //-------------------------------------------------------
  // missingMomentum
  //-------------------------------------------------------
  HepLorentzVector EventHeaderMaker::missingMomentum(MsgStream& log, double &sumET) {

    HepLorentzVector unused(0., 0., 0., 0.);
    const MissingMomentum* mm(0);
    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);

    
    double unused_sumET = mm->sumET();
    sumET = unused_sumET +
      +  transverseEnergySum<ReconstructedParticleCollection> (m_electronLocation)
      +  transverseEnergySum<ReconstructedParticleCollection> (m_photonLocation)
      +  transverseEnergySum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
      +  transverseEnergySum<JetCollection>                   (m_jetLocation)
      +  transverseEnergySum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
      -  assocTransverseEnergySum<JetCollection, ReconstructedParticle> (m_jetLocation);

    
    // Add halos if requested and if the missing ET has not already been adjusted
    // to treat isolation problem
    if (m_addHalos && !m_adjustMissETForIsolation){
      log << MSG::DEBUG <<"Adding halo values to missing ET, sumET" << endreq; 
      temp += haloVectorSum(m_electronLocation) + haloVectorSum(m_photonLocation);
      sumET += haloSum(m_electronLocation) + haloSum(m_photonLocation);
    }
    
    log << MSG::DEBUG<<"Total Imbalanced momentum = " << temp 
        << ", total sum of scalar ETs = " << sumET << endreq;
    
    if(m_testMode){
      HepLorentzVector totCluster = 
        momentumSum<IClusterCollection> (m_clusterLocation);
      
      HepLorentzVector clusterAssE = 
        assocMomentumSum<ReconstructedParticleCollection, Cluster>
          (m_electronLocation);
        HepLorentzVector clusterAssE2 = 
          assocMomentumSum<ReconstructedParticleCollection, KtCluster>
          (m_electronLocation);

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

        HepLorentzVector clusterAssM = 
          assocMomentumSum<ReconstructedParticleCollection, Cluster>
          (m_isolatedMuonLocation);
        HepLorentzVector clusterAssM2 = 
          assocMomentumSum<ReconstructedParticleCollection, KtCluster>
          (m_isolatedMuonLocation);

        HepLorentzVector clusterAssJ = 
          assocMomentumSum<JetCollection, Cluster>(m_jetLocation);
        HepLorentzVector clusterAssJ2 = 
          assocMomentumSum<JetCollection, KtCluster>(m_jetLocation);

      HepLorentzVector clusterDiff = totCluster-
        clusterAssJ-
        clusterAssE-
        clusterAssM-
        clusterAssP-
        clusterAssJ2-
        clusterAssE2-
        clusterAssM2-
        clusterAssP2;

      HepLorentzVector totCell = 
        momentumSum<ITwoCptCellCollection>(m_cellLocation);

      HepLorentzVector cellAssCluster = 
        assocMomentumSum<IClusterCollection, TwoCptCell>(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+clusterAssJ2<<endreq;
      log<<MSG::DEBUG<<"Clus Ass El      "<<clusterAssE+clusterAssE2<<endreq;
      log<<MSG::DEBUG<<"Clus Ass Mu      "<<clusterAssM+clusterAssM2<<endreq;
      log<<MSG::DEBUG<<"Clus Ass Ph      "<<clusterAssP+clusterAssP2<<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<IClusterCollection> (m_clusterLocation)
        +  momentumSum<ITwoCptCellCollection> (m_cellLocation)
        +  momentumSum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_electronLocation)   
        - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
        (m_electronLocation)   
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_photonLocation)
        - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
        (m_photonLocation)
        
        - assocMomentumSum<ReconstructedParticleCollection, Cluster>
        (m_nonIsolatedMuonLocation)
        - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
        (m_nonIsolatedMuonLocation)
        
        -  assocMomentumSum<JetCollection, Cluster>
        (m_jetLocation)
        -  assocMomentumSum<JetCollection, KtCluster>
        (m_jetLocation)
        
        -  assocMomentumSum<JetCollection, ReconstructedParticle>
        (m_jetLocation)

        - assocMomentumSum<IClusterCollection, 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

---