File indexing completed on 2024-05-10 02:21:24
0001 #include "SimG4Core/CheckSecondary/interface/StoreSecondary.h"
0002 #include "SimG4Core/CheckSecondary/interface/TreatSecondary.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006
0007 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0008 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0013 #include <CLHEP/Units/SystemOfUnits.h>
0014 #include "G4HCofThisEvent.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017
0018 #include <cmath>
0019 #include <iomanip>
0020 #include <iostream>
0021
0022 StoreSecondary::StoreSecondary(const edm::ParameterSet &p) {
0023 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("StoreSecondary");
0024 treatSecondary = new TreatSecondary(m_p);
0025
0026 produces<std::vector<math::XYZTLorentzVector>>("SecondaryMomenta");
0027 produces<std::vector<int>>("SecondaryParticles");
0028
0029 edm::LogInfo("CheckSecondary") << "Instantiate StoreSecondary to store "
0030 << "secondaries after 1st hadronic inelastic"
0031 << " interaction";
0032 }
0033
0034 StoreSecondary::~StoreSecondary() { delete treatSecondary; }
0035
0036 void StoreSecondary::produce(edm::Event &e, const edm::EventSetup &) {
0037 std::unique_ptr<std::vector<math::XYZTLorentzVector>> secMom(new std::vector<math::XYZTLorentzVector>);
0038 *secMom = secondaries;
0039 e.put(std::move(secMom), "SecondaryMomenta");
0040
0041 std::unique_ptr<std::vector<int>> secNumber(new std::vector<int>);
0042 *secNumber = nsecs;
0043 e.put(std::move(secNumber), "SecondaryParticles");
0044
0045 LogDebug("CheckSecondary") << "StoreSecondary:: Event " << e.id() << " with " << nsecs.size()
0046 << " hadronic collisions with "
0047 << "secondaries produced in each step";
0048 for (unsigned int i = 0; i < nsecs.size(); i++) {
0049 LogDebug("CheckSecondary") << " " << nsecs[i] << " from " << procs[i];
0050 }
0051 LogDebug("CheckSecondary") << " and " << secondaries.size() << " secondaries"
0052 << " produced in the first interactions:";
0053 for (unsigned int i = 0; i < secondaries.size(); i++)
0054 LogDebug("CheckSecondary") << "Secondary " << i << " " << secondaries[i];
0055 }
0056
0057 void StoreSecondary::update(const BeginOfEvent *) {
0058 nsecs.clear();
0059 procs.clear();
0060 secondaries.clear();
0061 }
0062
0063 void StoreSecondary::update(const BeginOfTrack *trk) {
0064 const G4Track *thTk = (*trk)();
0065 treatSecondary->initTrack(thTk);
0066 if (nsecs.empty() && thTk->GetParentID() <= 0)
0067 storeIt = true;
0068 else
0069 storeIt = false;
0070 nHad = 0;
0071 }
0072
0073 void StoreSecondary::update(const G4Step *aStep) {
0074 std::string name;
0075 int procID;
0076 bool hadrInt;
0077 double deltaE;
0078 std::vector<int> charge;
0079 std::vector<math::XYZTLorentzVector> tracks = treatSecondary->tracks(aStep, name, procID, hadrInt, deltaE, charge);
0080 if (hadrInt) {
0081 ++nHad;
0082 if (storeIt) {
0083 int sec = (int)(tracks.size());
0084 nsecs.push_back(sec);
0085 procs.push_back(name);
0086 if (nHad == 1) {
0087 for (int i = 0; i < sec; i++)
0088 secondaries.push_back(tracks[i]);
0089 }
0090 }
0091 }
0092 }
