#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.

References HLT_25ns14e33_v1_cff::A, Bremsstrahlung, HLT_25ns14e33_v1_cff::distCut, EnergyLoss, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), i, analyzePatCleaning_cfg::inputFile, j, MultipleScattering, MuonBremsstrahlung, NuclearInteraction, PairProduction, pTmin, plotFactory::ratios, AlCaHLTBitMon_QueryRunRegistry::string, theTECFudgeFactor, use_hardcoded, and Gflash::Z.

   : PairProduction(0), Bremsstrahlung(0),MuonBremsstrahlung(0),
     MultipleScattering(0), EnergyLoss(0), 
     NuclearInteraction(0),
     pTmin(999.), use_hardcoded(1)
 {
   // 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 ) { 
       NuclearInteraction = new NuclearInteractionFTFSimulator(distAlgo, distCut); 
     } else {
       NuclearInteraction = 
     new NuclearInteractionSimulator(hadronEnergies, hadronTypes, hadronNames, 
                     hadronMasses, hadronPMin, pionEnergy, 
                     lengthRatio, ratios, idMap, 
                     inputFile, distAlgo, distCut);
     }
   }
 }

MaterialEffects::~MaterialEffects ( )

Default destructor.

Definition at line 217 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 77 of file MaterialEffects.h.

References theEnergyLoss.

Referenced by TrajectoryManager::createPSimHits().

77 { return theEnergyLoss; }

MaterialEffects::theEnergyLoss

double theEnergyLoss

Definition: MaterialEffects.h:118

EnergyLossSimulator* MaterialEffects::energyLossSimulator ( ) const

inline

Return the Energy Loss engine.

Definition at line 85 of file MaterialEffects.h.

References EnergyLoss.

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

                                                           { 
     return EnergyLoss;
   }

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 228 of file MaterialEffects.cc.

References funct::abs(), FBaseSimEvent::addSimTrack(), FBaseSimEvent::addSimVertex(), MaterialEffectsSimulator::beginDaughters(), FSimVertexType::BREM_VERTEX, Bremsstrahlung, RawParticle::charge(), MaterialEffectsSimulator::closestDaughterId(), MaterialEffectsSimulator::endDaughters(), EnergyLoss, TrackerLayer::layerNumber(), MultipleScattering, MuonBremsstrahlung, MaterialEffectsSimulator::nDaughters(), normalVector(), FSimVertexType::NUCL_VERTEX, NuclearInteraction, FSimVertexType::PAIR_VERTEX, PairProduction, RawParticle::pid(), FSimVertex::position(), 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.pid()==22 ) {
     
     //
     PairProduction->updateState(myTrack, radlen, random);
 
     if ( PairProduction->nDaughters() ) {   
       //add a vertex to the mother particle
       int ivertex = mySimEvent.addSimVertex(myTrack.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.pid() == 22 ) return;
 
 //------------------------
 //   Nuclear interactions
 //------------------------ 
 
   if ( NuclearInteraction && abs(myTrack.pid()) > 100 
                           && abs(myTrack.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.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.charge() == 0 ) return;
 
   if ( !MuonBremsstrahlung && !Bremsstrahlung && !EnergyLoss && !MultipleScattering ) return;
 
 //----------------
 //  Bremsstrahlung
 //----------------
 
   if ( Bremsstrahlung && abs(myTrack.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.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.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.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.E();
     EnergyLoss->updateState(myTrack, radlen, random);
     theEnergyLoss -= myTrack.E();
   }
   
 
 
   if ( MultipleScattering && myTrack.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 80 of file MaterialEffects.h.

References MultipleScattering.

Referenced by MuonSimHitProducer::applyMaterialEffects().

                                                                           { 
     return MultipleScattering;
   }

MuonBremsstrahlungSimulator* MaterialEffects::muonBremsstrahlungSimulator ( ) const

inline

Return the Muon Bremsstrahlung engine.

Definition at line 90 of file MaterialEffects.h.

References MuonBremsstrahlung.

Referenced by MuonSimHitProducer::applyMaterialEffects().

                                                                           {
     return MuonBremsstrahlung;
   }

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

private

The vector normal to the surface traversed.

Definition at line 466 of file MaterialEffects.cc.

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

Referenced by interact().

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

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

private

The number of radiation lengths traversed.

Definition at line 425 of file MaterialEffects.cc.

References TrackerLayer::forward(), TrackerLayer::fudgeFactor(), TrackerLayer::fudgeMax(), TrackerLayer::fudgeMin(), TrackerLayer::fudgeNumber(), Surface::mediumProperties(), P, 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.Px(),myTrack.Py(),myTrack.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.R();
   double zed = fabs(myTrack.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 474 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();
 
 }