Skip to content

Commit

Permalink
added alternative nuclear interation using G4 FTF model
Browse files Browse the repository at this point in the history
  • Loading branch information
civanch committed Jan 27, 2015
1 parent e5c03a7 commit 0181a82
Show file tree
Hide file tree
Showing 8 changed files with 436 additions and 21 deletions.
2 changes: 2 additions & 0 deletions FastSimulation/MaterialEffects/BuildFile.xml
Expand Up @@ -3,6 +3,8 @@
<use name="FastSimulation/Event"/>
<use name="FastSimulation/Utilities"/>
<use name="rootcore"/>
<use name="clhep"/>
<use name="geant4core"/>
<export>
<lib name="1"/>
</export>
31 changes: 31 additions & 0 deletions FastSimulation/MaterialEffects/interface/CMSDummyDeexitation.h
@@ -0,0 +1,31 @@
#ifndef FastSimulation_MaterialEffects_CMSDummyDeexcitation_H
#define FastSimulation_MaterialEffects_CMSDummyDeexcitation_H

/**
* This class is needed as a dummy interface to Geant4
* nuclear de-excitation module; no secondary produced
*
* \author Vladimir Ivanchenko
* $Date: 20-Jan-2015
*/

#include "G4VPreCompoundModel.hh"
#include "G4ReactionProductVector.hh"

class G4Fragment;

class CMSDummyDeexcitation : public G4VPreCompoundModel
{
public:

CMSDummyDeexcitation() {};

virtual ~CMSDummyDeexcitation() {};

// virtual G4HadFinalState * ApplyYourself(const G4HadProjectile & thePrimary,
// G4Nucleus & theNucleus);

G4ReactionProductVector* DeExcite(G4Fragment&) { return new G4ReactionProductVector(); };

};
#endif
9 changes: 5 additions & 4 deletions FastSimulation/MaterialEffects/interface/MaterialEffects.h
Expand Up @@ -38,7 +38,8 @@ class ParticlePropagator;
class PairProductionSimulator;
class BremsstrahlungSimulator;
class EnergyLossSimulator;
class NuclearInteractionSimulator;
//class NuclearInteractionSimulator;
class MaterialEffectsSimulator;
class MultipleScatteringSimulator;
class MuonBremsstrahlungSimulator;
class RandomEngineAndDistribution;
Expand Down Expand Up @@ -85,7 +86,7 @@ class MaterialEffects
return EnergyLoss;
}

/// Return the Muon Bremsstrahlung engine
/// Return the Muon Bremsstrahlung engine
inline MuonBremsstrahlungSimulator* muonBremsstrahlungSimulator() const {
return MuonBremsstrahlung;
}
Expand All @@ -104,11 +105,11 @@ class MaterialEffects

PairProductionSimulator* PairProduction;
BremsstrahlungSimulator* Bremsstrahlung;
////// Muon Brem
////// Muon Brem
MuonBremsstrahlungSimulator* MuonBremsstrahlung;
MultipleScatteringSimulator* MultipleScattering;
EnergyLossSimulator* EnergyLoss;
NuclearInteractionSimulator* NuclearInteraction;
MaterialEffectsSimulator* NuclearInteraction;

// Cuts for material effects
double pTmin;
Expand Down
Expand Up @@ -52,7 +52,6 @@ class MaterialEffectsSimulator
///Electron mass in GeV/c2
inline double eMass() const { return 0.000510998902; }


/// Compute the material effect (calls the sub class)
void updateState(ParticlePropagator& myTrack, double radlen, RandomEngineAndDistribution const*);

Expand All @@ -74,6 +73,9 @@ class MaterialEffectsSimulator
/// The id of the closest charged daughter (filled for nuclear interactions only)
inline int closestDaughterId() { return theClosestChargedDaughterId; }

/// Used by NuclearInteractionSimulator to save last sampled event
virtual void save() {};

private:

/// Overloaded in all material effects updtators
Expand All @@ -82,7 +84,6 @@ class MaterialEffectsSimulator
/// Returns the fraction of radiation lengths traversed
inline double radiationLength() const {return radLengths;}


protected:

std::vector<RawParticle> _theUpdatedState;
Expand Down
@@ -0,0 +1,83 @@
#ifndef FastSimulation_MaterialEffects_NuclearInteractionFTFSimulator_H
#define FastSimulation_MaterialEffects_NuclearInteractionFTFSimulator_H

/**
* This class computes the probability for hadrons to interact with a
* nucleon of the tracker material (inelastically) and then sample
* nuclear interaction using FTF model of Geant4
* The fraction of interaction lengths traversed by the particle in this
* tracker layer is determined in MaterialEffectsSimulator as a fraction
* the radiation lengths.
*
* \author Vladimir Ivanchenko
* $Date: 20-Jan-2015
*/

#include "FastSimulation/MaterialEffects/interface/MaterialEffectsSimulator.h"

#include "G4Nucleus.hh"
#include "G4HadProjectile.hh"
#include "G4LorentzVector.hh"
#include "G4ThreeVector.hh"

#include <vector>

class ParticlePropagator;
class RandomEngineAndDistribution;
class G4ParticleDefinition;
class G4Track;
class G4Step;
class G4TheoFSGenerator;
class G4FTFModel;
class G4ExcitedStringDecay;
class G4LundStringFragmentation;
class G4GeneratorPrecompoundInterface;

class NuclearInteractionFTFSimulator : public MaterialEffectsSimulator
{
public:

/// Constructor
NuclearInteractionFTFSimulator(unsigned int distAlgo, double distCut);

/// Default Destructor
~NuclearInteractionFTFSimulator();

private:

/// Generate a nuclear interaction according to the probability that it happens
void compute(ParticlePropagator& Particle, RandomEngineAndDistribution const*);

void saveDaughter(ParticlePropagator& Particle, const G4LorentzVector& lv, int pdgid);

double distanceToPrimary(const RawParticle& Particle,
const RawParticle& aDaughter) const;

std::vector<const G4ParticleDefinition*> theG4Hadron;
std::vector<double> theNuclIntLength;
std::vector<int> theId;

G4TheoFSGenerator* theHadronicModel;
G4FTFModel* theStringModel;
G4ExcitedStringDecay* theStringDecay;
G4LundStringFragmentation* theLund;
G4GeneratorPrecompoundInterface* theCascade;

G4Step* dummyStep;
G4Track* currTrack;
const G4ParticleDefinition* currParticle;
G4Nucleus targetNucleus;
G4HadProjectile theProjectile;
G4LorentzVector curr4Mom;
G4ThreeVector vectProj;
G4ThreeVector theBoost;

double theEnergyLimit;

int numHadrons;

unsigned int theDistAlgo;
double theDistCut;
double distMin;
};
#endif
9 changes: 7 additions & 2 deletions FastSimulation/MaterialEffects/python/MaterialEffects_cfi.py
Expand Up @@ -41,8 +41,10 @@
# Smallest pT for the Mutliple Scattering
pTmin = cms.double(0.2),
# Enable Nuclear Interactions
NuclearInteraction = cms.bool(False), # temporary 12.01.14 until bug-fix for nuclear interactions
# The energies of the pions used in the above files (same order)
NuclearInteraction = cms.bool(False), # temporary 12.01.14 until bug-fix for nuclear interactions
#
G4NuclearInteraction = cms.bool(False),
# The energies of the pions used in the above files (same order)
pionEnergies = cms.untracked.vdouble(
1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 12.0, 15.0, 20.0,
30.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0
Expand Down Expand Up @@ -180,6 +182,7 @@
# Smallest pT for the Mutliple Scattering
pTmin = cms.double(0.3),
# Enable Nuclear Interactions
G4NuclearInteraction = cms.bool(False),
NuclearInteraction = cms.bool(False)

)
Expand Down Expand Up @@ -220,6 +223,7 @@
# Smallest pT for the Mutliple Scattering
pTmin = cms.double(0.3),
# Enable Nuclear Interactions
G4NuclearInteraction = cms.bool(False),
NuclearInteraction = cms.bool(False)
)
)
Expand Down Expand Up @@ -259,6 +263,7 @@
# Smallest pT for the Mutliple Scattering
pTmin = cms.double(0.3),
# Enable Nuclear Interactions
G4NuclearInteraction = cms.bool(False),
NuclearInteraction = cms.bool(False)

)
Expand Down
39 changes: 26 additions & 13 deletions FastSimulation/MaterialEffects/src/MaterialEffects.cc
Expand Up @@ -16,6 +16,7 @@
#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>
Expand All @@ -37,6 +38,7 @@ MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
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");
Expand All @@ -59,8 +61,8 @@ MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
Bremsstrahlung = new BremsstrahlungSimulator(bremEnergy,
bremEnergyFraction);
}
//muon Brem+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if ( doMuonBremsstrahlung ) {
//muon Brem+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if ( doMuonBremsstrahlung ) {

double bremEnergy = matEff.getParameter<double>("bremEnergy");
double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
Expand All @@ -69,8 +71,7 @@ MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)

}


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

if ( doEnergyLoss ) {

Expand All @@ -85,6 +86,7 @@ MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)

}


if ( doNuclearInteraction ) {

// The energies simulated
Expand Down Expand Up @@ -200,11 +202,15 @@ MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
idMap[idPiminusses[i]] = -211;

// Construction
NuclearInteraction =
new NuclearInteractionSimulator(pionEnergies, pionTypes, pionNames,
pionMasses, pionPMin, pionEnergy,
lengthRatio, ratios, idMap,
inputFile, distAlgo, distCut);
if ( doG4NuclInteraction ) {
NuclearInteraction = new NuclearInteractionFTFSimulator(distAlgo, distCut);
} else {
NuclearInteraction =
new NuclearInteractionSimulator(pionEnergies, pionTypes, pionNames,
pionMasses, pionPMin, pionEnergy,
lengthRatio, ratios, idMap,
inputFile, distAlgo, distCut);
}
}
}

Expand All @@ -215,7 +221,7 @@ MaterialEffects::~MaterialEffects() {
if ( EnergyLoss ) delete EnergyLoss;
if ( MultipleScattering ) delete MultipleScattering;
if ( NuclearInteraction ) delete NuclearInteraction;
//Muon Brem
//Muon Brem
if ( MuonBremsstrahlung ) delete MuonBremsstrahlung;
}

Expand All @@ -231,6 +237,10 @@ void MaterialEffects::interact(FSimEvent& mySimEvent,
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
//-------------------
Expand Down Expand Up @@ -283,12 +293,15 @@ void MaterialEffects::interact(FSimEvent& mySimEvent,
}
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
Expand All @@ -301,11 +314,11 @@ void MaterialEffects::interact(FSimEvent& mySimEvent,
int daughId = mySimEvent.addSimTrack(&(*DaughterIter), ivertex);

// Store the closest daughter in the mother info (for later tracking purposes)
if ( NuclearInteraction->closestDaughterId() == idaugh++ ) {
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;
Expand Down

0 comments on commit 0181a82

Please sign in to comment.