#include <MaterialEffects.h>

Public Member Functions
double	energyLoss () const
	Return the energy loss by ionization in the current layer. More...

EnergyLossSimulator *	energyLossSimulator () const
	Return the Energy Loss engine. More...

void	interact (FSimEvent &simEvent, const TrackerLayer &layer, ParticlePropagator &PP, unsigned i, RandomEngineAndDistribution const *)

	MaterialEffects (const edm::ParameterSet &matEff)
	Constructor. More...

MultipleScatteringSimulator *	multipleScatteringSimulator () const
	Return the Multiple Scattering engine. More...

MuonBremsstrahlungSimulator *	muonBremsstrahlungSimulator () const
	Return the Muon Bremsstrahlung engine. More...

void	save ()
	Save nuclear interaction information. More...

double	thickness () const
	Return the thickness of the current layer. More...

	~MaterialEffects ()
	Default destructor. More...

Private Member Functions
GlobalVector	normalVector (const TrackerLayer &layer, ParticlePropagator &myTrack) const
	The vector normal to the surface traversed. More...

double	radLengths (const TrackerLayer &layer, ParticlePropagator &myTrack)
	The number of radiation lengths traversed. More...

Private Attributes
BremsstrahlungSimulator *	Bremsstrahlung

EnergyLossSimulator *	EnergyLoss

MultipleScatteringSimulator *	MultipleScattering

MuonBremsstrahlungSimulator *	MuonBremsstrahlung

MaterialEffectsSimulator *	NuclearInteraction

PairProductionSimulator *	PairProduction

double	pTmin

double	theEnergyLoss

GlobalVector	theNormalVector

double	theTECFudgeFactor

double	theThickness

bool	use_hardcoded

Detailed Description

Definition at line 51 of file MaterialEffects.h.

Constructor & Destructor Documentation

MaterialEffects::MaterialEffects ( const edm::ParameterSet & matEff )

Constructor.

Definition at line 26 of file MaterialEffects.cc.

     : 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 ( )

Default destructor.

Definition at line 197 of file MaterialEffects.cc.

References Bremsstrahlung, EnergyLoss, MultipleScattering, MuonBremsstrahlung, NuclearInteraction, and PairProduction.

                                   {
   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;
 }

Member Function Documentation

double MaterialEffects::energyLoss ( ) const

inline

Return the energy loss by ionization in the current layer.

Definition at line 74 of file MaterialEffects.h.

Referenced by TrajectoryManager::createPSimHits().

74 { return theEnergyLoss; }

MaterialEffects::theEnergyLoss

double theEnergyLoss

Definition: MaterialEffects.h:105

EnergyLossSimulator* MaterialEffects::energyLossSimulator ( ) const

inline

Return the Energy Loss engine.

Definition at line 80 of file MaterialEffects.h.

References MaterialEffects_cfi::EnergyLoss.

Referenced by MuonSimHitProducer::applyMaterialEffects(), and CalorimetryManager::MuonMipSimulation().

80 { return EnergyLoss; }

MaterialEffects::EnergyLoss

EnergyLossSimulator * EnergyLoss

Definition: MaterialEffects.h:98

void MaterialEffects::interact	(	FSimEvent &	simEvent,
		const TrackerLayer &	layer,
		ParticlePropagator &	PP,
		unsigned	i,
		RandomEngineAndDistribution const *	random
	)

Steer the various interaction processes in the Tracker Material and update the FSimEvent

Energy loss

Multiple scattering

Definition at line 213 of file MaterialEffects.cc.

References funct::abs(), FBaseSimEvent::addSimTrack(), FBaseSimEvent::addSimVertex(), MaterialEffectsSimulator::beginDaughters(), FSimVertexType::BREM_VERTEX, Bremsstrahlung, RawParticle::charge(), MaterialEffectsSimulator::closestDaughterId(), RawParticle::E(), MaterialEffectsSimulator::endDaughters(), EnergyLoss, DQMScaleToClient_cfi::factor, TrackerLayer::layerNumber(), MultipleScattering, MuonBremsstrahlung, MaterialEffectsSimulator::nDaughters(), normalVector(), FSimVertexType::NUCL_VERTEX, NuclearInteraction, FSimVertexType::PAIR_VERTEX, PairProduction, BaseParticlePropagator::particle(), RawParticle::pid(), FSimVertex::position(), RawParticle::Pt(), pTmin, radLengths(), FSimTrack::setClosestDaughterId(), MaterialEffectsSimulator::setNormalVector(), theEnergyLoss, theNormalVector, theTECFudgeFactor, FBaseSimEvent::track(), MaterialEffectsSimulator::updateState(), use_hardcoded, FSimTrack::vertex(), and RawParticle::vertex().

Referenced by TrajectoryManager::reconstruct().

                                                                           {
   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);
   }
 }

MultipleScatteringSimulator* MaterialEffects::multipleScatteringSimulator ( ) const

inline

Return the Multiple Scattering engine.

Definition at line 77 of file MaterialEffects.h.

References MaterialEffects_cfi::MultipleScattering.

Referenced by MuonSimHitProducer::applyMaterialEffects().

77 { return MultipleScattering; }

MaterialEffects::MultipleScattering

MultipleScatteringSimulator * MultipleScattering

Definition: MaterialEffects.h:97

MuonBremsstrahlungSimulator* MaterialEffects::muonBremsstrahlungSimulator ( ) const

inline

Return the Muon Bremsstrahlung engine.

Definition at line 83 of file MaterialEffects.h.

References MaterialEffects_cfi::MuonBremsstrahlung.

Referenced by MuonSimHitProducer::applyMaterialEffects().

83 { return MuonBremsstrahlung; }

MaterialEffects::MuonBremsstrahlung

MuonBremsstrahlungSimulator * MuonBremsstrahlung

Definition: MaterialEffects.h:96

GlobalVector MaterialEffects::normalVector	(	const TrackerLayer &	layer,
		ParticlePropagator &	myTrack
	)		const

private

The vector normal to the surface traversed.

Definition at line 416 of file MaterialEffects.cc.

References TrackerLayer::disk(), TrackerLayer::forward(), BaseParticlePropagator::particle(), RawParticle::R(), RawParticle::X(), and RawParticle::Y().

Referenced by interact().

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

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

private

The number of radiation lengths traversed.

Definition at line 381 of file MaterialEffects.cc.

References DQMScaleToClient_cfi::factor, TrackerLayer::forward(), TrackerLayer::fudgeFactor(), TrackerLayer::fudgeMax(), TrackerLayer::fudgeMin(), TrackerLayer::fudgeNumber(), Surface::mediumProperties(), BaseParticlePropagator::particle(), RawParticle::Px(), RawParticle::Py(), RawParticle::Pz(), RawParticle::R(), MediumProperties::radLen(), TrackerLayer::surface(), theNormalVector, theThickness, and RawParticle::Z().

Referenced by interact().

                                                                                          {
   // 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;
 }

void MaterialEffects::save ( )

Save nuclear interaction information.

Definition at line 421 of file MaterialEffects.cc.

References NuclearInteraction, and MaterialEffectsSimulator::save().

Referenced by Vispa.Main.TabController.TabController::allowClose(), Vispa.Main.TabController.TabController::checkModificationTimestamp(), and TrajectoryManager::reconstruct().

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