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