MaterialEffects.cc

Go to the documentation of this file.

//Framework Headers
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

//TrackingTools Headers

// Famos Headers
#include "FastSimulation/Event/interface/FSimEvent.h"
#include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
#include "FastSimulation/TrackerSetup/interface/TrackerLayer.h"
#include "FastSimulation/MaterialEffects/interface/MaterialEffects.h"

#include "FastSimulation/MaterialEffects/interface/PairProductionSimulator.h"
#include "FastSimulation/MaterialEffects/interface/MultipleScatteringSimulator.h"
#include "FastSimulation/MaterialEffects/interface/BremsstrahlungSimulator.h"
#include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
#include "FastSimulation/MaterialEffects/interface/NuclearInteractionSimulator.h"
#include "FastSimulation/MaterialEffects/interface/NuclearInteractionFTFSimulator.h"
#include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"

#include <list>
#include <map>
#include <string>

MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
  : PairProduction(nullptr), Bremsstrahlung(nullptr),MuonBremsstrahlung(nullptr),
    MultipleScattering(nullptr), EnergyLoss(nullptr), 
    NuclearInteraction(nullptr),
    pTmin(999.), use_hardcoded(true)
{
  // Set the minimal photon energy for a Brem from e+/-

  use_hardcoded = matEff.getParameter<bool >("use_hardcoded_geometry");

  bool doPairProduction     = matEff.getParameter<bool>("PairProduction");
  bool doBremsstrahlung     = matEff.getParameter<bool>("Bremsstrahlung");
  bool doEnergyLoss         = matEff.getParameter<bool>("EnergyLoss");
  bool doMultipleScattering = matEff.getParameter<bool>("MultipleScattering");
  bool doNuclearInteraction = matEff.getParameter<bool>("NuclearInteraction");
  bool doG4NuclInteraction  = matEff.getParameter<bool>("G4NuclearInteraction");
  bool doMuonBremsstrahlung = matEff.getParameter<bool>("MuonBremsstrahlung");

  double A = matEff.getParameter<double>("A");
  double Z = matEff.getParameter<double>("Z");
  double density = matEff.getParameter<double>("Density");
  double radLen = matEff.getParameter<double>("RadiationLength");
  
  // Set the minimal pT before giving up the dE/dx treatment

  if ( doPairProduction ) { 

    double photonEnergy = matEff.getParameter<double>("photonEnergy");
    PairProduction = new PairProductionSimulator(photonEnergy);
  }

  if ( doBremsstrahlung ) { 

    double bremEnergy = matEff.getParameter<double>("bremEnergy");
    double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
    Bremsstrahlung = new BremsstrahlungSimulator(bremEnergy,
                         bremEnergyFraction);
  }
  //muon Brem+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if ( doMuonBremsstrahlung ) {

    double bremEnergy = matEff.getParameter<double>("bremEnergy");
    double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
    MuonBremsstrahlung = new MuonBremsstrahlungSimulator(A,Z,density,radLen,bremEnergy,
                                                 bremEnergyFraction);

  }

  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if ( doEnergyLoss ) { 

    pTmin = matEff.getParameter<double>("pTmin");
    EnergyLoss = new EnergyLossSimulator(A,Z,density,radLen);

  }

  if ( doMultipleScattering ) { 
    
    MultipleScattering = new MultipleScatteringSimulator(A,Z,density,radLen);

  }


  if ( doNuclearInteraction ) { 

    // The energies simulated
    std::vector<double> hadronEnergies 
      = matEff.getUntrackedParameter<std::vector<double> >("hadronEnergies");

    // The particle types simulated
    std::vector<int> hadronTypes 
      = matEff.getUntrackedParameter<std::vector<int> >("hadronTypes");

    // The corresponding particle names
    std::vector<std::string> hadronNames 
      = matEff.getUntrackedParameter<std::vector<std::string> >("hadronNames");

    // The corresponding particle masses
    std::vector<double> hadronMasses 
      = matEff.getUntrackedParameter<std::vector<double> >("hadronMasses");

    // The smallest momentum for inelastic interactions
    std::vector<double> hadronPMin 
      = matEff.getUntrackedParameter<std::vector<double> >("hadronMinP");

    // The interaction length / radiation length ratio for each particle type
    std::vector<double> lengthRatio 
      = matEff.getParameter<std::vector<double> >("lengthRatio");
    //    std::map<int,double> lengthRatio;
    //    for ( unsigned i=0; i<theLengthRatio.size(); ++i )
    //      lengthRatio[ hadronTypes[i] ] = theLengthRatio[i];

    // A global fudge factor for TEC layers (which apparently do not react to 
    // hadrons the same way as all other layers...
    theTECFudgeFactor = matEff.getParameter<double>("fudgeFactor");

    // The evolution of the interaction lengths with energy
    std::vector<double> theRatios  
      = matEff.getUntrackedParameter<std::vector<double> >("ratios");
    //std::map<int,std::vector<double> > ratios;
    //for ( unsigned i=0; i<hadronTypes.size(); ++i ) { 
    //  for ( unsigned j=0; j<hadronEnergies.size(); ++j ) { 
    //  ratios[ hadronTypes[i] ].push_back(theRatios[ i*hadronEnergies.size() + j ]);
    //  }
    //}
    std::vector< std::vector<double> > ratios;
    ratios.resize(hadronTypes.size());
    for ( unsigned i=0; i<hadronTypes.size(); ++i ) { 
      for ( unsigned j=0; j<hadronEnergies.size(); ++j ) { 
    ratios[i].push_back(theRatios[ i*hadronEnergies.size() + j ]);
      }
    }

    // The smallest momentum for elastic interactions
    double pionEnergy 
      = matEff.getParameter<double>("pionEnergy");

    // The algorithm to compute the distance between primary and secondaries
    // when a nuclear interaction occurs
    unsigned distAlgo 
      = matEff.getParameter<unsigned>("distAlgo");
    double distCut 
      = matEff.getParameter<double>("distCut");

    // The file to read the starting interaction in each files
    // (random reproducibility in case of a crash)
    std::string inputFile 
      = matEff.getUntrackedParameter<std::string>("inputFile");

    // Build the ID map (i.e., what is to be considered as a proton, etc...)
    std::map<int,int> idMap;
    // Protons
    std::vector<int> idProtons 
      = matEff.getUntrackedParameter<std::vector<int> >("protons");
    for ( unsigned i=0; i<idProtons.size(); ++i ) 
      idMap[idProtons[i]] = 2212;
    // Anti-Protons
    std::vector<int> idAntiProtons 
      = matEff.getUntrackedParameter<std::vector<int> >("antiprotons");
    for ( unsigned i=0; i<idAntiProtons.size(); ++i ) 
      idMap[idAntiProtons[i]] = -2212;
    // Neutrons
    std::vector<int> idNeutrons 
      = matEff.getUntrackedParameter<std::vector<int> >("neutrons");
    for ( unsigned i=0; i<idNeutrons.size(); ++i ) 
      idMap[idNeutrons[i]] = 2112;
    // Anti-Neutrons
    std::vector<int> idAntiNeutrons 
      = matEff.getUntrackedParameter<std::vector<int> >("antineutrons");
    for ( unsigned i=0; i<idAntiNeutrons.size(); ++i ) 
      idMap[idAntiNeutrons[i]] = -2112;
    // K0L's
    std::vector<int> idK0Ls 
      = matEff.getUntrackedParameter<std::vector<int> >("K0Ls");
    for ( unsigned i=0; i<idK0Ls.size(); ++i ) 
      idMap[idK0Ls[i]] = 130;
    // K+'s
    std::vector<int> idKplusses 
      = matEff.getUntrackedParameter<std::vector<int> >("Kplusses");
    for ( unsigned i=0; i<idKplusses.size(); ++i ) 
      idMap[idKplusses[i]] = 321;
    // K-'s
    std::vector<int> idKminusses 
      = matEff.getUntrackedParameter<std::vector<int> >("Kminusses");
    for ( unsigned i=0; i<idKminusses.size(); ++i ) 
      idMap[idKminusses[i]] = -321;
    // pi+'s
    std::vector<int> idPiplusses 
      = matEff.getUntrackedParameter<std::vector<int> >("Piplusses");
    for ( unsigned i=0; i<idPiplusses.size(); ++i ) 
      idMap[idPiplusses[i]] = 211;
    // pi-'s
    std::vector<int> idPiminusses 
      = matEff.getUntrackedParameter<std::vector<int> >("Piminusses");
    for ( unsigned i=0; i<idPiminusses.size(); ++i ) 
      idMap[idPiminusses[i]] = -211;

    // Construction
    if ( doG4NuclInteraction ) { 
      double elimit = matEff.getParameter<double>("EkinBertiniGeV")*CLHEP::GeV;
      double eth = matEff.getParameter<double>("EkinLimitGeV")*CLHEP::GeV;
      NuclearInteraction = 
    new NuclearInteractionFTFSimulator(distAlgo, distCut, elimit, eth); 
    } else {
      NuclearInteraction = 
    new NuclearInteractionSimulator(hadronEnergies, hadronTypes, hadronNames, 
                    hadronMasses, hadronPMin, pionEnergy, 
                    lengthRatio, ratios, idMap, 
                    inputFile, distAlgo, distCut);
    }
  }
}

MaterialEffects::~MaterialEffects() {

  if ( PairProduction ) delete PairProduction;
  if ( Bremsstrahlung ) delete Bremsstrahlung;
  if ( EnergyLoss ) delete EnergyLoss;
  if ( MultipleScattering ) delete MultipleScattering;
  if ( NuclearInteraction ) delete NuclearInteraction;
  //Muon Brem
  if ( MuonBremsstrahlung ) delete MuonBremsstrahlung;
}

void MaterialEffects::interact(FSimEvent& mySimEvent, 
                   const TrackerLayer& layer,
                   ParticlePropagator& myTrack,
                   unsigned itrack,
                               RandomEngineAndDistribution const* random) {

  MaterialEffectsSimulator::RHEP_const_iter DaughterIter;
  double radlen;
  theEnergyLoss = 0;
  theNormalVector = normalVector(layer,myTrack);
  radlen = radLengths(layer,myTrack);

  //std::cout << "### MaterialEffects: for Track= " <<  itrack << " in layer #"
  //        << layer.layerNumber() <<  std::endl;
  //std::cout << myTrack << std::endl;

//-------------------
//  Photon Conversion
//-------------------

  if ( PairProduction && myTrack.particle().pid()==22 ) {
    
    //
    PairProduction->updateState(myTrack, radlen, random);

    if ( PairProduction->nDaughters() ) {   
      //add a vertex to the mother particle
      int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(),itrack,
                        FSimVertexType::PAIR_VERTEX);
      
      // Check if it is a valid vertex first:
      if (ivertex>=0) {
    // This was a photon that converted
    for ( DaughterIter = PairProduction->beginDaughters();
          DaughterIter != PairProduction->endDaughters(); 
          ++DaughterIter) {
      
      mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
      
    }
    // The photon converted. Return.
    return;
      }
      else {
    edm::LogWarning("MaterialEffects") <<  " WARNING: A non valid vertex was found in photon conv. -> " << ivertex << std::endl;    
      }

    }

  }

  if ( myTrack.particle().pid() == 22 ) return;

//------------------------
//   Nuclear interactions
//------------------------ 

  if ( NuclearInteraction && abs(myTrack.particle().pid()) > 100 
                          && abs(myTrack.particle().pid()) < 1000000) { 

    // Simulate a nuclear interaction
    double factor = 1.0;
    if(use_hardcoded){
      if (layer.layerNumber() >= 19 && layer.layerNumber() <= 27 ) 
    factor = theTECFudgeFactor;
    }
    NuclearInteraction->updateState(myTrack, radlen*factor, random);

    //std::cout << "MaterialEffects: nDaughters= " 
    //        << NuclearInteraction->nDaughters() << std::endl;
    if ( NuclearInteraction->nDaughters() ) { 

      //add a end vertex to the mother particle
      int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(),itrack,
                        FSimVertexType::NUCL_VERTEX);
      //std::cout << "ivertex= " << ivertex << " nDaughters= " 
      //    << NuclearInteraction->nDaughters() << std::endl;
      // Check if it is a valid vertex first:
      if (ivertex>=0) {
    // This was a hadron that interacted inelastically
    int idaugh = 0;
    for ( DaughterIter = NuclearInteraction->beginDaughters();
          DaughterIter != NuclearInteraction->endDaughters(); 
          ++DaughterIter) {
      
      // The daughter in the event
      int daughId = mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
      
      // Store the closest daughter in the mother info (for later tracking purposes)
      if ( NuclearInteraction->closestDaughterId() == idaugh ) {
        if ( mySimEvent.track(itrack).vertex().position().Pt() < 4.0 ) 
          mySimEvent.track(itrack).setClosestDaughterId(daughId);
      }
      ++idaugh;
    }
    // The hadron is destroyed. Return.
    return;
      }
      else {
    edm::LogWarning("MaterialEffects") <<  " WARNING: A non valid vertex was found in nucl. int. -> " << ivertex << std::endl;    
      }

    }
    
  }

  if ( myTrack.particle().charge() == 0 ) return;

  if ( !MuonBremsstrahlung && !Bremsstrahlung && !EnergyLoss && !MultipleScattering ) return;

//----------------
//  Bremsstrahlung
//----------------

  if ( Bremsstrahlung && abs(myTrack.particle().pid())==11 ) {
        
    Bremsstrahlung->updateState(myTrack,radlen, random);

    if ( Bremsstrahlung->nDaughters() ) {
      
      // Add a vertex, but do not attach it to the electron, because it 
      // continues its way...
      int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(),itrack,
                        FSimVertexType::BREM_VERTEX);

      // Check if it is a valid vertex first:
      if (ivertex>=0) {
    for ( DaughterIter = Bremsstrahlung->beginDaughters();
          DaughterIter != Bremsstrahlung->endDaughters(); 
          ++DaughterIter) {
      mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
    }
      }
      else {
    edm::LogWarning("MaterialEffects") <<  " WARNING: A non valid vertex was found in brem -> " << ivertex << std::endl;    
      }
      
    }
    
  } 

//---------------------------
//  Muon_Bremsstrahlung
//--------------------------

  if (  MuonBremsstrahlung && abs(myTrack.particle().pid())==13 ) {
       
    MuonBremsstrahlung->updateState(myTrack, radlen, random);

    if ( MuonBremsstrahlung->nDaughters() ) {

      // Add a vertex, but do not attach it to the muon, because it 
      // continues its way...
      int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(),itrack,
                                            FSimVertexType::BREM_VERTEX);
 
     // Check if it is a valid vertex first:
      if (ivertex>=0) {
    for ( DaughterIter = MuonBremsstrahlung->beginDaughters();
          DaughterIter != MuonBremsstrahlung->endDaughters();
          ++DaughterIter) {
      mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
    }
      }
      else {
    edm::LogWarning("MaterialEffects") <<  " WARNING: A non valid vertex was found in muon brem -> " << ivertex << std::endl;    
      }

    }

  }


  if ( EnergyLoss )
  {
    theEnergyLoss = myTrack.particle().E();
    EnergyLoss->updateState(myTrack, radlen, random);
    theEnergyLoss -= myTrack.particle().E();
  }
  


  if ( MultipleScattering && myTrack.particle().Pt() > pTmin ) {
    //    MultipleScattering->setNormalVector(normalVector(layer,myTrack));
    MultipleScattering->setNormalVector(theNormalVector);
    MultipleScattering->updateState(myTrack,radlen, random);
  }
    
}

double
MaterialEffects::radLengths(const TrackerLayer& layer,
                ParticlePropagator& myTrack) {

  // Thickness of layer
  theThickness = layer.surface().mediumProperties().radLen();

  GlobalVector P(myTrack.particle().Px(),myTrack.particle().Py(),myTrack.particle().Pz());
  
  // Effective length of track inside layer (considering crossing angle)
  //  double radlen = theThickness / fabs(P.dot(theNormalVector)/(P.mag()*theNormalVector.mag()));
  double radlen = theThickness / fabs(P.dot(theNormalVector)) * P.mag();

  // This is a series of fudge factors (from the geometry description), 
  // to describe the layer inhomogeneities (services, cables, supports...)
  double rad = myTrack.particle().R();
  double zed = fabs(myTrack.particle().Z());

  double factor = 1;

  // Are there fudge factors for this layer
  if ( layer.fudgeNumber() ) 

    // If yes, loop on them
    for ( unsigned int iLayer=0; iLayer < layer.fudgeNumber(); ++iLayer ) { 

      // Apply to R if forward layer, to Z if barrel layer
      if ( (  layer.forward() && layer.fudgeMin(iLayer) < rad && rad < layer.fudgeMax(iLayer) )  ||
       ( !layer.forward() && layer.fudgeMin(iLayer) < zed && zed < layer.fudgeMax(iLayer) ) ) {
    factor = layer.fudgeFactor(iLayer);
    break;
      }
    
  }

  theThickness *= factor;

  return radlen * factor;

}

GlobalVector
MaterialEffects::normalVector(const TrackerLayer& layer,
                  ParticlePropagator& myTrack ) const {
  return layer.forward() ?  
    layer.disk()->normalVector() :
    GlobalVector(myTrack.particle().X(),myTrack.particle().Y(),0.)/myTrack.particle().R();
}

void 
MaterialEffects::save() { 

  // Save current nuclear interactions in the event libraries.
  if ( NuclearInteraction ) NuclearInteraction->save();

}