initial commit

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py

@@ -102,6 +102,9 @@ def getFastCaloSimSvcV2(name="ISF_FastCaloSimSvcV2", **kwargs):
     kwargs.setdefault("CaloCellMakerTools_release"       , [ 'ISF_CaloCellContainerFCSFinalizerTool',
                                                            'ISF_FastHitConvertTool' ])
     kwargs.setdefault("FastCaloSimCaloExtrapolation"     , 'FastCaloSimCaloExtrapolation')
+    kwargs.setdefault("PunchThroughTool"                 , 'ISF_PunchThroughTool'             )
+    kwargs.setdefault("DoPunchThroughSimulation"         , True                                             )
+    kwargs.setdefault("ParticleBroker"                   , 'ISF_ParticleBrokerSvc'               )
 
     kwargs.setdefault("ParamSvc", "ISF_FastCaloSimV2ParamSvc")
     # register the FastCaloSim random number streams

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx

@@ -27,6 +27,8 @@
 // StoreGate
 #include "StoreGate/StoreGateSvc.h"
 #include "StoreGate/StoreGate.h"
+#include "ISF_Interfaces/IParticleBroker.h"
+
 
 #include "CaloEvent/CaloCellContainer.h"
 #include "CaloDetDescr/CaloDetDescrElement.h"
@@ -35,6 +37,7 @@
 
 #include "PathResolver/PathResolver.h"
 
+
 #include "TFile.h"
 #include <fstream>
 
@@ -46,14 +49,21 @@ ISF::FastCaloSimSvcV2::FastCaloSimSvcV2(const std::string& name, ISvcLocator* sv
   : BaseSimulationSvc(name, svc)
   , m_paramSvc("ISF_FastCaloSimV2ParamSvc", name)
   , m_rndGenSvc("AtRndmGenSvc", name)
+  , m_doPunchThrough(false)
+  , m_punchThroughTool("")
+  , m_particleBroker ("ISF_ParticleBroker",name)
+
 {
   declareProperty("ParamSvc"                       ,       m_paramSvc);
   declareProperty("CaloCellsOutputName"            ,       m_caloCellsOutputName) ;
   declareProperty("CaloCellMakerTools_setup"       ,       m_caloCellMakerToolsSetup) ;
   declareProperty("CaloCellMakerTools_release"     ,       m_caloCellMakerToolsRelease) ;
+  declareProperty("PunchThroughTool"               ,       m_punchThroughTool);
+  declareProperty("DoPunchThroughSimulation"       ,       m_doPunchThrough) ;
   declareProperty("RandomSvc"                      ,       m_rndGenSvc                );
   declareProperty("RandomStream"                   ,       m_randomEngineName         );
   declareProperty("FastCaloSimCaloExtrapolation"   ,       m_FastCaloSimCaloExtrapolation );
+  declareProperty("ParticleBroker"                 ,       m_particleBroker, "ISF ParticleBroker Svc" );
 }
 
 /** framework methods */
@@ -63,11 +73,18 @@ StatusCode ISF::FastCaloSimSvcV2::initialize()
 
   ATH_CHECK(m_rndGenSvc.retrieve());
   m_randomEngine = m_rndGenSvc->GetEngine( m_randomEngineName);
+
   if (!m_randomEngine) {
     ATH_MSG_ERROR("Could not get random number engine from RandomNumberService. Abort.");
     return StatusCode::FAILURE;
   }
 
+  if (m_doPunchThrough && m_punchThroughTool.retrieve().isFailure() ) 
+  {
+    ATH_MSG_ERROR (m_punchThroughTool.propertyName() << ": Failed to retrieve tool " << m_punchThroughTool.type());
+    return StatusCode::FAILURE;
+  } 
+
   ATH_CHECK(m_paramSvc.retrieve());
 
   // Get FastCaloSimCaloExtrapolation
@@ -148,6 +165,18 @@ StatusCode ISF::FastCaloSimSvcV2::simulate(const ISF::ISFParticle& isfp)
   Amg::Vector3D particle_position =  isfp.position();
   Amg::Vector3D particle_direction(isfp.momentum().x(),isfp.momentum().y(),isfp.momentum().z());
 
+  if (m_doPunchThrough) {
+     // call punch-through simulation
+     const ISF::ISFParticleContainer* isfpVec = m_punchThroughTool->computePunchThroughParticles(isfp);
+     if (isfpVec) {
+       ISF::ISFParticleContainer::const_iterator partIt    = isfpVec->begin();
+       ISF::ISFParticleContainer::const_iterator partItEnd = isfpVec->end();
+       for ( ; partIt!=partItEnd; ++partIt) {
+         m_particleBroker->push( *partIt, &isfp);
+       }
+     }
+   }
+
    //int barcode=isfp.barcode(); // isfp barcode, eta and phi: in case we need them
   // float eta_isfp = particle_position.eta();
   // float phi_isfp = particle_position.phi();

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.h

@@ -8,6 +8,12 @@
 // ISF includes
 #include "ISF_Interfaces/BaseSimulationSvc.h"
 
+// Framework includes
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/IChronoStatSvc.h"
+#include "AthenaBaseComps/AthService.h"
+
 // FastCaloSim includes
 #include "IFastCaloSimParamSvc.h"
 #include "ISF_FastCaloSimParametrization/IFastCaloSimCaloExtrapolation.h"
@@ -23,6 +29,9 @@
 #include "CaloIdentifier/LArFCAL_ID.h"
 #include "CaloIdentifier/TileID.h"
 
+#include "ISF_FastCaloSimInterfaces/IPunchThroughTool.h"
+
+
 namespace CLHEP
 {
   class HepRandomEngine;
@@ -34,6 +43,10 @@ class TFCSParametrizationBase;
 
 
 namespace ISF {
+
+  class IParticleBroker;
+  class IPunchThroughTool;
+
   /** @class FastCaloSimSvcV2
 
       @author Elmar.Ritsch -at- cern.ch, Geraldine.Conti -at- cern.ch, Flavia.Dias -at- cern.ch
@@ -75,6 +88,11 @@ namespace ISF {
     std::string                     m_randomEngineName;
 
     std::string  m_caloCellsOutputName;
+
+    bool m_doPunchThrough;
+    ToolHandle< IPunchThroughTool >     m_punchThroughTool;
+    ServiceHandle<ISF::IParticleBroker> m_particleBroker;
+
   };
 
 }

Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx

@@ -27,7 +27,7 @@ double ISF::PDFcreator::getRand(std::vector<double> inputParameters, bool discre
   // we have to choose which bin we use for the current energy & eta regime.
   // It's very unlikely to exactly hit a bin's center, therefore we randomly
   // choose either the lower next or the upper next bin.
-
+ 
   int numInputPar = inputParameters.size();
   TAxis **axis = new TAxis*[numInputPar];
   Int_t *chosenBin = new Int_t[numInputPar];
@@ -35,8 +35,8 @@ double ISF::PDFcreator::getRand(std::vector<double> inputParameters, bool discre
 
   // get the axis from the first (could be any other too, but the first one
   // always exists) histogram (which holds the function's first fit parameter)
-  axis[0] = m_par[0]->GetXaxis();
-  axis[1] = m_par[0]->GetYaxis();
+  axis[0] = m_par[0]->GetXaxis(); //Energy
+  axis[1] = m_par[0]->GetYaxis(); //Eta
 
   // loop over all input inputParameters (e.g. eta & energy)
   for (int inputPar=0; inputPar<numInputPar; inputPar++) {
@@ -47,6 +47,7 @@ double ISF::PDFcreator::getRand(std::vector<double> inputParameters, bool discre
     // get the id of the bin corresponding to the current value
     // (use FindFixBin to avoid the axis from being rebinned)
     Int_t bin = curaxis->FindFixBin(curvalue);
+
     // get the center of the closest bin to the input inputParameter
     double closestCenter = curaxis->GetBinCenter(bin);
     // get the bins edge closest to the current value
@@ -78,10 +79,11 @@ double ISF::PDFcreator::getRand(std::vector<double> inputParameters, bool discre
   // of them)
   Int_t bin = 0;
   if (numInputPar == 1)             bin = m_par[0]->GetBin( chosenBin[0] );
-  else if (numInputPar ==2 )        bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1] );
-  else if (numInputPar == 3)        bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1], chosenBin[2] );
+  else if (numInputPar ==2 )        bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1], 1 );
+  else if (numInputPar == 3)        bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1], chosenBin[2] ); //Select Z bin as 1 to choose pions for now
   // TODO: implement case of >3 input parameters
 
+
   // free some memory
   delete [] axis;
 

Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx

@@ -400,6 +400,8 @@ const ISF::ISFParticleContainer* ISF::PunchThroughTool::computePunchThroughParti
   //  -> the azimutal angle phi
   m_initPhi = m_initPs->position().phi();
 
+  m_initPdg = m_initPs->pdgCode();
+
   // this is the place where the magic is done:
   // test for each registered punch-through pdg if a punch-through
   // occures and create these particles
@@ -487,6 +489,7 @@ int ISF::PunchThroughTool::getAllParticles(int pdg, int numParticles) const
       std::vector<double> parameters;
       parameters.push_back( m_initEnergy );
       parameters.push_back( fabs(m_initEta) );
+
       // the maximum number of particles which should be produced
       // if no maximum number is given, this is -1
       int maxParticles = p->getMaxNumParticles();
@@ -634,6 +637,7 @@ ISF::ISFParticle *ISF::PunchThroughTool::getOneParticle(int pdg, double maxEnerg
   std::vector<double> parInitEnergyEta;
   parInitEnergyEta.push_back( m_initEnergy );
   parInitEnergyEta.push_back( fabs(m_initEta) );
+  //parInitEnergyEta.push_back( fabs(m_initEta) ); ADD initial pdgid here
 
   // (2.1) get the energy
   double energy = p->getExitEnergyPDF()->getRand(
@@ -905,7 +909,7 @@ ISF::PDFcreator *ISF::PunchThroughTool::readLookuptablePDF(int pdg, std::string
       name.str("");
       name << folderName << pdg<< "/parameter" << par;
 
-      // get the histogram from the file
+      // get the histogram from the filenan
       TH1 *hist = (TH1*)m_fileLookupTable->Get(name.str().c_str());
 
       // if histogram does not exist -> error!

Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.h

@@ -10,6 +10,7 @@
 
 // Athena Base
 #include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/IAlgTool.h"
 
 #include "BarcodeEvent/Barcode.h"
 #include "BarcodeEvent/PhysicsProcessCode.h"
@@ -115,6 +116,7 @@ namespace ISF {
     mutable double                       m_initEta{0.};     //!< the incoming particle's eta
     mutable double                       m_initTheta{0.};   //!< the incoming particle's theta
     mutable double                       m_initPhi{0.};     //!< the incoming particle's phi
+    mutable int                       m_initPdg{0};     //!< the incoming particle's pdg
 
     /** calo-MS borders */
     double                               m_R1{0.};