TruthClusterizationFactory.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
 
 
 
#include "SiClusterizationTool/TruthClusterizationFactory.h"
 
//for position estimate and clustering
#include "InDetIdentifier/PixelID.h"
#include "InDetPrepRawData/PixelCluster.h"
 
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/Incident.h"
#include "InDetSimData/InDetSimData.h"
#include "AtlasHepMC/GenParticle.h"
#include "AthenaKernel/RNGWrapper.h"
#include "CxxUtils/checker_macros.h"
 
#include "TrkEventPrimitives/ParamDefs.h"
 
#include "CLHEP/Random/RandFlat.h"
 
#include <TMath.h>
 
namespace InDet {
 
  TruthClusterizationFactory::TruthClusterizationFactory(const std::string& name,
                                                         const std::string& n, const IInterface* p):
    AthAlgTool(name, n, p)
  {
    declareInterface<TruthClusterizationFactory>(this);
  } 
  
  StatusCode TruthClusterizationFactory::initialize() {
 
    ATH_MSG_DEBUG(m_truthClusterSplittingEff);
 
    // random svc
    CHECK(m_rndmSvc.retrieve());
  
    // get the random stream
    ATH_MSG_DEBUG ( "Getting random number engine : <" << m_rndmEngineName << ">" );
    m_rndmEngine = m_rndmSvc->getEngine(this, m_rndmEngineName);
    if (!m_rndmEngine) {
      ATH_MSG_ERROR("Could not find RndmEngine : " << m_rndmEngineName);
      return StatusCode::FAILURE;
    }
    else {
      ATH_MSG_DEBUG("Found RndmEngine : " << m_rndmEngineName);
    }
 
    ATH_CHECK( m_simDataCollectionName.initialize() );
 
    return StatusCode::SUCCESS;
  }
  
 
 
  std::vector<double> TruthClusterizationFactory::estimateNumberOfParticles(const InDet::PixelCluster& pCluster) const
  {
    const EventContext& ctx = Gaudi::Hive::currentContext();
    if (ctx.evt() != m_rndmEngine->evtSeeded()) {
       ATHRNG::RNGWrapper* wrapper ATLAS_THREAD_SAFE = m_rndmEngine;
       wrapper->setSeed (this->name(), ctx);
    }
    CLHEP::HepRandomEngine* engine = m_rndmEngine->getEngine (ctx);
    
    std::vector<double> probabilities(3,0.);
    const auto &rdos = pCluster.rdoList();
    unsigned int nPartContributing = 0;
    if (m_discardPUHits) {
      // Initialize set for a list of distinct uniqueIDs for the cluster
      std::set<int> uniqueIDs;
      SG::ReadHandle<InDetSimDataCollection> pixSdoColl(m_simDataCollectionName);
      //Loop over all elements (pixels/strips) in the cluster
      if(pixSdoColl.isValid()){
        for (auto rdoIter :  rdos){
          auto simDataIter = pixSdoColl->find(rdoIter);
          if (simDataIter != pixSdoColl->end()){
            // get the SimData and count the individual contributions
            auto simData = (simDataIter->second);
            for( const auto& deposit : simData.getdeposits() ){
              //If deposit exists
              if (!deposit.first){ATH_MSG_DEBUG("No deposits found"); continue;}
              // This should only be used for samples without pile-up
              if (deposit.first.eventIndex() != 0) continue;
              uniqueIDs.insert(deposit.first.id());
            }
          }
        }
      }
      //uniqueIDs lists the unique truth particles contributing to the
      //cluster
      nPartContributing = uniqueIDs.size();
    }
    else {
      //Initialize set for a list of distinct uniqueIDs for the cluster
      //- if we are taking into account pile-up, then we need
      //GenEvent::event_number() + GenParticle::id() to uniquely
      //idenftify a GenParticle from a group of GenEvents.
      std::set< std::pair<HepMcParticleLink::index_type, int> > uniqueIDs;
      SG::ReadHandle<InDetSimDataCollection> pixSdoColl(m_simDataCollectionName);
      //Loop over all elements (pixels/strips) in the cluster
      if(pixSdoColl.isValid()){
        for (auto rdoIter :  rdos){
          auto simDataIter = pixSdoColl->find(rdoIter);
          if (simDataIter != pixSdoColl->end()){
            // get the SimData and count the individual contributions
            auto simData = (simDataIter->second);
            for( const auto& deposit : simData.getdeposits() ){
              //If deposit exists
              if (!deposit.first){ATH_MSG_DEBUG("No deposits found"); continue;}
              // This should only be used for samples with pile-up
              uniqueIDs.insert(std::make_pair(deposit.first.eventIndex(), deposit.first.id()));
            }
          }
        }
      }
      //uniqueIDs lists the unique truth particles contributing to the
      //cluster
      nPartContributing = uniqueIDs.size();
    }
    ATH_MSG_VERBOSE("n Part Contributing: " << nPartContributing);
    ATH_MSG_VERBOSE("Smearing TruthClusterizationFactory probability output for TIDE studies");
    //If only 1 truth particles found
    //For pure PU case nPartContributing=0, assume that there is a single particle contributing as well
    if (nPartContributing<=1) {
      //NN will always return 100% chance of there being only 1 particle
      probabilities[0] = 1.0;
    }
    //If two unique truth particles found in cluster
    else if (nPartContributing==2) {
      //90% chance NN returns high probability of there being 2 particles
      if (CLHEP::RandFlat::shoot( engine, 0, 1 ) < m_truthClusterSplittingEff) probabilities[1] = 1.0;
      //Other 10% NN returns high probability of there being 1 particle
      else probabilities[0] = 1.0;
    }
    //If greater than 2 unique truth particles in cluster
    else if (nPartContributing>2) {
      //90% chance NN returns high probability of there being >2 particles
      if (CLHEP::RandFlat::shoot( engine, 0, 1 ) < m_truthClusterSplittingEff) probabilities[2] = 1.0;
      //Other 10% NN returns high probability of there being 1 particle
      else probabilities[0] = 1.0;
    }
    
    return probabilities;
 
  }
 
  std::vector<Amg::Vector2D> TruthClusterizationFactory::estimatePositions(const InDet::PixelCluster& ) const
  { 
    ATH_MSG_ERROR("TruthClusterizationFactory::estimatePositions called for ITk ambiguity setup, should never happen! Digital clustering should be run for positions & errors.");
    return std::vector<Amg::Vector2D>(2,Amg::Vector2D (2,0.));
  }
  
}//end InDet namespace