Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-09 22:38:41

0001 #include "DBremWatcher.h"
0002 
0003 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0004 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0005 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0006 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0007 #include "SimG4Core/Notification/interface/TrackInformation.h"
0008 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0009 #include "SimG4Core/Notification/interface/Observer.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "G4Region.hh"
0013 #include "G4RegionStore.hh"
0014 #include "G4LogicalVolumeStore.hh"
0015 #include "G4ProductionCuts.hh"
0016 #include "G4ProcessTable.hh"
0017 #include "G4ProcessManager.hh"
0018 #include "G4MuonMinus.hh"
0019 #include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
0020 #include "G4Track.hh"
0021 #include "G4PhysicalConstants.hh"
0022 #include "G4SystemOfUnits.hh"
0023 #include "G4VProcess.hh"
0024 #include "G4VParticleChange.hh"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include <algorithm>
0027 
0028 DBremWatcher::DBremWatcher(edm::ParameterSet const& p) {
0029   edm::ParameterSet ps = p.getParameter<edm::ParameterSet>("DBremWatcher");
0030   biasFactor = ps.getUntrackedParameter<double>("DBremBiasFactor", 1);
0031   m_weight = 0;
0032   foundbrem = false;
0033   finaltraj = G4ThreeVector(0, 0, 0);
0034   aPrimeTraj = G4ThreeVector(0, 0, 0);
0035   VertexPos = G4ThreeVector(0, 0, 0);
0036   f_energy = 0;
0037   produces<float>("DBremEventWeight");
0038   produces<float>("DBremLocationX");
0039   produces<float>("DBremLocationY");
0040   produces<float>("DBremLocationZ");
0041   //produces<std::string>("DBremMaterial");
0042   produces<float>("DBremAngle");
0043   produces<float>("DBremInitialEnergy");
0044   produces<float>("DBremFinalEnergy");
0045   produces<float>("BiasFactor");
0046   pdgs_ = ps.getUntrackedParameter<std::vector<int>>("PDGCodes");
0047   edm::LogInfo("DBremWatcher") << "DBremWatcher:: Save Sim Track if PDG code "
0048                                << "is one from the list of " << pdgs_.size() << " items";
0049   for (unsigned int k = 0; k < pdgs_.size(); ++k)
0050     edm::LogInfo("DBremWatcher") << "[" << k << "] " << pdgs_[k];
0051 }
0052 
0053 DBremWatcher::~DBremWatcher() {}
0054 
0055 void DBremWatcher::update(const BeginOfTrack* trk) {
0056   G4Track* theTrack = (G4Track*)((*trk)());
0057   TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0058   if (trkInfo) {
0059     int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0060     G4ThreeVector Vpos = theTrack->GetVertexPosition();
0061     const G4VProcess* TrPro = theTrack->GetCreatorProcess();
0062     if (TrPro != nullptr) {
0063       if ((theTrack->GetCreatorProcess()->GetProcessName()) == "muDBrem") {
0064         if (std::find(pdgs_.begin(), pdgs_.end(), pdg) == pdgs_.end()) {
0065           //Found the deflected muon track
0066           trkInfo->storeTrack(true);
0067           if (!theTrack->IsGoodForTracking()) {
0068             theTrack->SetGoodForTrackingFlag(true);
0069           }
0070           f_energy = theTrack->GetTotalEnergy();
0071           foundbrem = true;
0072           finaltraj = theTrack->GetMomentum();
0073         } else {
0074           m_weight = theTrack->GetWeight();
0075         }
0076       }
0077     }
0078     if (std::find(pdgs_.begin(), pdgs_.end(), pdg) != pdgs_.end()) {
0079       //Found an A'
0080       trkInfo->storeTrack(true);
0081       VertexPos = Vpos;
0082       aPrimeTraj = theTrack->GetMomentum();
0083       LogDebug("DBremWatcher") << "Save SimTrack the Track " << theTrack->GetTrackID() << " Type "
0084                                << theTrack->GetDefinition()->GetParticleName() << " Momentum "
0085                                << theTrack->GetMomentum() / MeV << " MeV/c";
0086     }
0087   }
0088 }
0089 
0090 void DBremWatcher::update(const BeginOfRun* run) {}
0091 
0092 void DBremWatcher::update(const BeginOfEvent* event) {
0093   G4String pname = "muDBrem";
0094   G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
0095   G4bool state = true;
0096   ptable->SetProcessActivation(pname, state);
0097   foundbrem = false;
0098 }
0099 
0100 void DBremWatcher::update(const EndOfEvent* event) {}
0101 
0102 void DBremWatcher::update(const EndOfTrack* trk) {
0103   G4Track* theTrack = (G4Track*)((*trk)());
0104   TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0105   const G4VProcess* TrPro = theTrack->GetCreatorProcess();
0106   if (trkInfo && TrPro != nullptr) {
0107     int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0108 
0109     if (std::find(pdgs_.begin(), pdgs_.end(), pdg) == pdgs_.end() &&
0110         (theTrack->GetCreatorProcess()->GetProcessName()) == "muDBrem") {
0111       trkInfo->storeTrack(true);
0112     }
0113   }
0114 }
0115 
0116 void DBremWatcher::produce(edm::Event& fEvent, const edm::EventSetup&) {
0117   if (foundbrem) {
0118     std::unique_ptr<float> weight = std::make_unique<float>(m_weight);
0119     fEvent.put(std::move(weight), "DBremEventWeight");
0120     std::unique_ptr<float> vtxposx = std::make_unique<float>(VertexPos.x());
0121     std::unique_ptr<float> vtxposy = std::make_unique<float>(VertexPos.y());
0122     std::unique_ptr<float> vtxposz = std::make_unique<float>(VertexPos.z());
0123     fEvent.put(std::move(vtxposx), "DBremLocationX");
0124     fEvent.put(std::move(vtxposy), "DBremLocationY");
0125     fEvent.put(std::move(vtxposz), "DBremLocationZ");
0126     std::unique_ptr<float> finalE = std::make_unique<float>(f_energy / GeV);
0127     fEvent.put(std::move(finalE), "DBremFinalEnergy");
0128     float deflectionAngle = -1;
0129     float initialEnergy = sqrt(pow(finaltraj.x() + aPrimeTraj.x(), 2) + pow(finaltraj.y() + aPrimeTraj.y(), 2) +
0130                                pow(finaltraj.z() + aPrimeTraj.z(), 2) + pow(0.1056 * GeV, 2));
0131     G4ThreeVector mother(
0132         finaltraj.x() + aPrimeTraj.x(), finaltraj.y() + aPrimeTraj.y(), finaltraj.z() + aPrimeTraj.z());
0133     deflectionAngle = mother.angle(finaltraj);
0134     std::unique_ptr<float> dAngle = std::make_unique<float>(deflectionAngle);
0135     std::unique_ptr<float> initialE = std::make_unique<float>(initialEnergy / GeV);
0136     fEvent.put(std::move(dAngle), "DBremAngle");
0137     fEvent.put(std::move(initialE), "DBremInitialEnergy");
0138     std::unique_ptr<float> bias = std::make_unique<float>(biasFactor);
0139     fEvent.put(std::move(bias), "BiasFactor");
0140   } else {
0141     std::unique_ptr<float> weight = std::make_unique<float>(0.);
0142     fEvent.put(std::move(weight), "DBremEventWeight");
0143   }
0144 }
0145 
0146 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0147 
0148 DEFINE_SIMWATCHER(DBremWatcher);