MaterialEffects/src/MaterialEffects.cc

0001
0002 //Framework Headers
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 //TrackingTools Headers
0007
0008 // Famos Headers
0009 #include "FastSimulation/Event/interface/FSimEvent.h"
0010 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0011 #include "FastSimulation/TrackerSetup/interface/TrackerLayer.h"
0012 #include "FastSimulation/MaterialEffects/interface/MaterialEffects.h"
0013
0014 #include "FastSimulation/MaterialEffects/interface/PairProductionSimulator.h"
0015 #include "FastSimulation/MaterialEffects/interface/MultipleScatteringSimulator.h"
0016 #include "FastSimulation/MaterialEffects/interface/BremsstrahlungSimulator.h"
0017 #include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
0018 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionSimulator.h"
0019 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionFTFSimulator.h"
0020 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0021
0022 #include <list>
0023 #include <map>
0024 #include <string>
0025
0026 MaterialEffects::MaterialEffects(const edm::ParameterSet& matEff)
0027     : PairProduction(nullptr),
0028       Bremsstrahlung(nullptr),
0029       MuonBremsstrahlung(nullptr),
0030       MultipleScattering(nullptr),
0031       EnergyLoss(nullptr),
0032       NuclearInteraction(nullptr),
0033       pTmin(999.),
0034       use_hardcoded(true) {
0035   // Set the minimal photon energy for a Brem from e+/-
0036
0037   use_hardcoded = matEff.getParameter<bool>("use_hardcoded_geometry");
0038
0039   bool doPairProduction = matEff.getParameter<bool>("PairProduction");
0040   bool doBremsstrahlung = matEff.getParameter<bool>("Bremsstrahlung");
0041   bool doEnergyLoss = matEff.getParameter<bool>("EnergyLoss");
0042   bool doMultipleScattering = matEff.getParameter<bool>("MultipleScattering");
0043   bool doNuclearInteraction = matEff.getParameter<bool>("NuclearInteraction");
0044   bool doG4NuclInteraction = matEff.getParameter<bool>("G4NuclearInteraction");
0045   bool doMuonBremsstrahlung = matEff.getParameter<bool>("MuonBremsstrahlung");
0046
0047   double A = matEff.getParameter<double>("A");
0048   double Z = matEff.getParameter<double>("Z");
0049   double density = matEff.getParameter<double>("Density");
0050   double radLen = matEff.getParameter<double>("RadiationLength");
0051
0052   // Set the minimal pT before giving up the dE/dx treatment
0053
0054   if (doPairProduction) {
0055     double photonEnergy = matEff.getParameter<double>("photonEnergy");
0056     PairProduction = new PairProductionSimulator(photonEnergy);
0057   }
0058
0059   if (doBremsstrahlung) {
0060     double bremEnergy = matEff.getParameter<double>("bremEnergy");
0061     double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
0062     Bremsstrahlung = new BremsstrahlungSimulator(bremEnergy, bremEnergyFraction);
0063   }
0064   //muon Brem+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0065   if (doMuonBremsstrahlung) {
0066     double bremEnergy = matEff.getParameter<double>("bremEnergy");
0067     double bremEnergyFraction = matEff.getParameter<double>("bremEnergyFraction");
0068     MuonBremsstrahlung = new MuonBremsstrahlungSimulator(A, Z, density, radLen, bremEnergy, bremEnergyFraction);
0069   }
0070
0071   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0072
0073   if (doEnergyLoss) {
0074     pTmin = matEff.getParameter<double>("pTmin");
0075     EnergyLoss = new EnergyLossSimulator(A, Z, density, radLen);
0076   }
0077
0078   if (doMultipleScattering) {
0079     MultipleScattering = new MultipleScatteringSimulator(A, Z, density, radLen);
0080   }
0081
0082   if (doNuclearInteraction) {
0083     // The energies simulated
0084     std::vector<double> hadronEnergies = matEff.getUntrackedParameter<std::vector<double> >("hadronEnergies");
0085
0086     // The particle types simulated
0087     std::vector<int> hadronTypes = matEff.getUntrackedParameter<std::vector<int> >("hadronTypes");
0088
0089     // The corresponding particle names
0090     std::vector<std::string> hadronNames = matEff.getUntrackedParameter<std::vector<std::string> >("hadronNames");
0091
0092     // The corresponding particle masses
0093     std::vector<double> hadronMasses = matEff.getUntrackedParameter<std::vector<double> >("hadronMasses");
0094
0095     // The smallest momentum for inelastic interactions
0096     std::vector<double> hadronPMin = matEff.getUntrackedParameter<std::vector<double> >("hadronMinP");
0097
0098     // The interaction length / radiation length ratio for each particle type
0099     std::vector<double> lengthRatio = matEff.getParameter<std::vector<double> >("lengthRatio");
0100     //    std::map<int,double> lengthRatio;
0101     //    for ( unsigned i=0; i<theLengthRatio.size(); ++i )
0102     //      lengthRatio[ hadronTypes[i] ] = theLengthRatio[i];
0103
0104     // A global fudge factor for TEC layers (which apparently do not react to
0105     // hadrons the same way as all other layers...
0106     theTECFudgeFactor = matEff.getParameter<double>("fudgeFactor");
0107
0108     // The evolution of the interaction lengths with energy
0109     std::vector<double> theRatios = matEff.getUntrackedParameter<std::vector<double> >("ratios");
0110     //std::map<int,std::vector<double> > ratios;
0111     //for ( unsigned i=0; i<hadronTypes.size(); ++i ) {
0112     //  for ( unsigned j=0; j<hadronEnergies.size(); ++j ) {
0113     //  ratios[ hadronTypes[i] ].push_back(theRatios[ i*hadronEnergies.size() + j ]);
0114     //  }
0115     //}
0116     std::vector<std::vector<double> > ratios;
0117     ratios.resize(hadronTypes.size());
0118     for (unsigned i = 0; i < hadronTypes.size(); ++i) {
0119       for (unsigned j = 0; j < hadronEnergies.size(); ++j) {
0120         ratios[i].push_back(theRatios[i * hadronEnergies.size() + j]);
0121       }
0122     }
0123
0124     // The smallest momentum for elastic interactions
0125     double pionEnergy = matEff.getParameter<double>("pionEnergy");
0126
0127     // The algorithm to compute the distance between primary and secondaries
0128     // when a nuclear interaction occurs
0129     unsigned distAlgo = matEff.getParameter<unsigned>("distAlgo");
0130     double distCut = matEff.getParameter<double>("distCut");
0131
0132     // The file to read the starting interaction in each files
0133     // (random reproducibility in case of a crash)
0134     std::string inputFile = matEff.getUntrackedParameter<std::string>("inputFile");
0135
0136     // Build the ID map (i.e., what is to be considered as a proton, etc...)
0137     std::map<int, int> idMap;
0138     // Protons
0139     std::vector<int> idProtons = matEff.getUntrackedParameter<std::vector<int> >("protons");
0140     for (unsigned i = 0; i < idProtons.size(); ++i)
0141       idMap[idProtons[i]] = 2212;
0142     // Anti-Protons
0143     std::vector<int> idAntiProtons = matEff.getUntrackedParameter<std::vector<int> >("antiprotons");
0144     for (unsigned i = 0; i < idAntiProtons.size(); ++i)
0145       idMap[idAntiProtons[i]] = -2212;
0146     // Neutrons
0147     std::vector<int> idNeutrons = matEff.getUntrackedParameter<std::vector<int> >("neutrons");
0148     for (unsigned i = 0; i < idNeutrons.size(); ++i)
0149       idMap[idNeutrons[i]] = 2112;
0150     // Anti-Neutrons
0151     std::vector<int> idAntiNeutrons = matEff.getUntrackedParameter<std::vector<int> >("antineutrons");
0152     for (unsigned i = 0; i < idAntiNeutrons.size(); ++i)
0153       idMap[idAntiNeutrons[i]] = -2112;
0154     // K0L's
0155     std::vector<int> idK0Ls = matEff.getUntrackedParameter<std::vector<int> >("K0Ls");
0156     for (unsigned i = 0; i < idK0Ls.size(); ++i)
0157       idMap[idK0Ls[i]] = 130;
0158     // K+'s
0159     std::vector<int> idKplusses = matEff.getUntrackedParameter<std::vector<int> >("Kplusses");
0160     for (unsigned i = 0; i < idKplusses.size(); ++i)
0161       idMap[idKplusses[i]] = 321;
0162     // K-'s
0163     std::vector<int> idKminusses = matEff.getUntrackedParameter<std::vector<int> >("Kminusses");
0164     for (unsigned i = 0; i < idKminusses.size(); ++i)
0165       idMap[idKminusses[i]] = -321;
0166     // pi+'s
0167     std::vector<int> idPiplusses = matEff.getUntrackedParameter<std::vector<int> >("Piplusses");
0168     for (unsigned i = 0; i < idPiplusses.size(); ++i)
0169       idMap[idPiplusses[i]] = 211;
0170     // pi-'s
0171     std::vector<int> idPiminusses = matEff.getUntrackedParameter<std::vector<int> >("Piminusses");
0172     for (unsigned i = 0; i < idPiminusses.size(); ++i)
0173       idMap[idPiminusses[i]] = -211;
0174
0175     // Construction
0176     if (doG4NuclInteraction) {
0177       double elimit = matEff.getParameter<double>("EkinBertiniGeV") * CLHEP::GeV;
0178       double eth = matEff.getParameter<double>("EkinLimitGeV") * CLHEP::GeV;
0179       NuclearInteraction = new NuclearInteractionFTFSimulator(distAlgo, distCut, elimit, eth);
0180     } else {
0181       NuclearInteraction = new NuclearInteractionSimulator(hadronEnergies,
0182                                                            hadronTypes,
0183                                                            hadronNames,
0184                                                            hadronMasses,
0185                                                            hadronPMin,
0186                                                            pionEnergy,
0187                                                            lengthRatio,
0188                                                            ratios,
0189                                                            idMap,
0190                                                            inputFile,
0191                                                            distAlgo,
0192                                                            distCut);
0193     }
0194   }
0195 }
0196
0197 MaterialEffects::~MaterialEffects() {
0198   if (PairProduction)
0199     delete PairProduction;
0200   if (Bremsstrahlung)
0201     delete Bremsstrahlung;
0202   if (EnergyLoss)
0203     delete EnergyLoss;
0204   if (MultipleScattering)
0205     delete MultipleScattering;
0206   if (NuclearInteraction)
0207     delete NuclearInteraction;
0208   //Muon Brem
0209   if (MuonBremsstrahlung)
0210     delete MuonBremsstrahlung;
0211 }
0212
0213 void MaterialEffects::interact(FSimEvent& mySimEvent,
0214                                const TrackerLayer& layer,
0215                                ParticlePropagator& myTrack,
0216                                unsigned itrack,
0217                                RandomEngineAndDistribution const* random) {
0218   MaterialEffectsSimulator::RHEP_const_iter DaughterIter;
0219   double radlen;
0220   theEnergyLoss = 0;
0221   theNormalVector = normalVector(layer, myTrack);
0222   radlen = radLengths(layer, myTrack);
0223
0224   //std::cout << "### MaterialEffects: for Track= " <<  itrack << " in layer #"
0225   //        << layer.layerNumber() <<  std::endl;
0226   //std::cout << myTrack << std::endl;
0227
0228   //-------------------
0229   //  Photon Conversion
0230   //-------------------
0231
0232   if (PairProduction && myTrack.particle().pid() == 22) {
0233     //
0234     PairProduction->updateState(myTrack, radlen, random);
0235
0236     if (PairProduction->nDaughters()) {
0237       //add a vertex to the mother particle
0238       int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::PAIR_VERTEX);
0239
0240       // Check if it is a valid vertex first:
0241       if (ivertex >= 0) {
0242         // This was a photon that converted
0243         for (DaughterIter = PairProduction->beginDaughters(); DaughterIter != PairProduction->endDaughters();
0244              ++DaughterIter) {
0245           mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0246         }
0247         // The photon converted. Return.
0248         return;
0249       } else {
0250         edm::LogWarning("MaterialEffects")
0251             << " WARNING: A non valid vertex was found in photon conv. -> " << ivertex << std::endl;
0252       }
0253     }
0254   }
0255
0256   if (myTrack.particle().pid() == 22)
0257     return;
0258
0259   //------------------------
0260   //   Nuclear interactions
0261   //------------------------
0262
0263   if (NuclearInteraction && abs(myTrack.particle().pid()) > 100 && abs(myTrack.particle().pid()) < 1000000) {
0264     // Simulate a nuclear interaction
0265     double factor = 1.0;
0266     if (use_hardcoded) {
0267       if (layer.layerNumber() >= 19 && layer.layerNumber() <= 27)
0268         factor = theTECFudgeFactor;
0269     }
0270     NuclearInteraction->updateState(myTrack, radlen * factor, random);
0271
0272     //std::cout << "MaterialEffects: nDaughters= "
0273     //        << NuclearInteraction->nDaughters() << std::endl;
0274     if (NuclearInteraction->nDaughters()) {
0275       //add a end vertex to the mother particle
0276       int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::NUCL_VERTEX);
0277       //std::cout << "ivertex= " << ivertex << " nDaughters= "
0278       //    << NuclearInteraction->nDaughters() << std::endl;
0279       // Check if it is a valid vertex first:
0280       if (ivertex >= 0) {
0281         // This was a hadron that interacted inelastically
0282         int idaugh = 0;
0283         for (DaughterIter = NuclearInteraction->beginDaughters(); DaughterIter != NuclearInteraction->endDaughters();
0284              ++DaughterIter) {
0285           // The daughter in the event
0286           int daughId = mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0287
0288           // Store the closest daughter in the mother info (for later tracking purposes)
0289           if (NuclearInteraction->closestDaughterId() == idaugh) {
0290             if (mySimEvent.track(itrack).vertex().position().Pt() < 4.0)
0291               mySimEvent.track(itrack).setClosestDaughterId(daughId);
0292           }
0293           ++idaugh;
0294         }
0295         // The hadron is destroyed. Return.
0296         return;
0297       } else {
0298         edm::LogWarning("MaterialEffects")
0299             << " WARNING: A non valid vertex was found in nucl. int. -> " << ivertex << std::endl;
0300       }
0301     }
0302   }
0303
0304   if (myTrack.particle().charge() == 0)
0305     return;
0306
0307   if (!MuonBremsstrahlung && !Bremsstrahlung && !EnergyLoss && !MultipleScattering)
0308     return;
0309
0310   //----------------
0311   //  Bremsstrahlung
0312   //----------------
0313
0314   if (Bremsstrahlung && abs(myTrack.particle().pid()) == 11) {
0315     Bremsstrahlung->updateState(myTrack, radlen, random);
0316
0317     if (Bremsstrahlung->nDaughters()) {
0318       // Add a vertex, but do not attach it to the electron, because it
0319       // continues its way...
0320       int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::BREM_VERTEX);
0321
0322       // Check if it is a valid vertex first:
0323       if (ivertex >= 0) {
0324         for (DaughterIter = Bremsstrahlung->beginDaughters(); DaughterIter != Bremsstrahlung->endDaughters();
0325              ++DaughterIter) {
0326           mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0327         }
0328       } else {
0329         edm::LogWarning("MaterialEffects")
0330             << " WARNING: A non valid vertex was found in brem -> " << ivertex << std::endl;
0331       }
0332     }
0333   }
0334
0335   //---------------------------
0336   //  Muon_Bremsstrahlung
0337   //--------------------------
0338
0339   if (MuonBremsstrahlung && abs(myTrack.particle().pid()) == 13) {
0340     MuonBremsstrahlung->updateState(myTrack, radlen, random);
0341
0342     if (MuonBremsstrahlung->nDaughters()) {
0343       // Add a vertex, but do not attach it to the muon, because it
0344       // continues its way...
0345       int ivertex = mySimEvent.addSimVertex(myTrack.particle().vertex(), itrack, FSimVertexType::BREM_VERTEX);
0346
0347       // Check if it is a valid vertex first:
0348       if (ivertex >= 0) {
0349         for (DaughterIter = MuonBremsstrahlung->beginDaughters(); DaughterIter != MuonBremsstrahlung->endDaughters();
0350              ++DaughterIter) {
0351           mySimEvent.addSimTrack(&(*DaughterIter), ivertex);
0352         }
0353       } else {
0354         edm::LogWarning("MaterialEffects")
0355             << " WARNING: A non valid vertex was found in muon brem -> " << ivertex << std::endl;
0356       }
0357     }
0358   }
0359
0360   ////--------------
0361   ////  Energy loss
0362   ///---------------
0363
0364   if (EnergyLoss) {
0365     theEnergyLoss = myTrack.particle().E();
0366     EnergyLoss->updateState(myTrack, radlen, random);
0367     theEnergyLoss -= myTrack.particle().E();
0368   }
0369
0370   ////----------------------
0371   ////  Multiple scattering
0372   ///-----------------------
0373
0374   if (MultipleScattering && myTrack.particle().Pt() > pTmin) {
0375     //    MultipleScattering->setNormalVector(normalVector(layer,myTrack));
0376     MultipleScattering->setNormalVector(theNormalVector);
0377     MultipleScattering->updateState(myTrack, radlen, random);
0378   }
0379 }
0380
0381 double MaterialEffects::radLengths(const TrackerLayer& layer, ParticlePropagator& myTrack) {
0382   // Thickness of layer
0383   theThickness = layer.surface().mediumProperties().radLen();
0384
0385   GlobalVector P(myTrack.particle().Px(), myTrack.particle().Py(), myTrack.particle().Pz());
0386
0387   // Effective length of track inside layer (considering crossing angle)
0388   //  double radlen = theThickness / fabs(P.dot(theNormalVector)/(P.mag()*theNormalVector.mag()));
0389   double radlen = theThickness / fabs(P.dot(theNormalVector)) * P.mag();
0390
0391   // This is a series of fudge factors (from the geometry description),
0392   // to describe the layer inhomogeneities (services, cables, supports...)
0393   double rad = myTrack.particle().R();
0394   double zed = fabs(myTrack.particle().Z());
0395
0396   double factor = 1;
0397
0398   // Are there fudge factors for this layer
0399   if (layer.fudgeNumber())
0400
0401     // If yes, loop on them
0402     for (unsigned int iLayer = 0; iLayer < layer.fudgeNumber(); ++iLayer) {
0403       // Apply to R if forward layer, to Z if barrel layer
0404       if ((layer.forward() && layer.fudgeMin(iLayer) < rad && rad < layer.fudgeMax(iLayer)) ||
0405           (!layer.forward() && layer.fudgeMin(iLayer) < zed && zed < layer.fudgeMax(iLayer))) {
0406         factor = layer.fudgeFactor(iLayer);
0407         break;
0408       }
0409     }
0410
0411   theThickness *= factor;
0412
0413   return radlen * factor;
0414 }
0415
0416 GlobalVector MaterialEffects::normalVector(const TrackerLayer& layer, ParticlePropagator& myTrack) const {
0417   return layer.forward() ? layer.disk()->normalVector()
0418                          : GlobalVector(myTrack.particle().X(), myTrack.particle().Y(), 0.) / myTrack.particle().R();
0419 }
0420
0421 void MaterialEffects::save() {
0422   // Save current nuclear interactions in the event libraries.
0423   if (NuclearInteraction)
0424     NuclearInteraction->save();
0425 }