67 edm::LogVerbatim(
"ForwardSim") <<
"============================================================================";
68 edm::LogVerbatim(
"ForwardSim") <<
"CastorTestAnalysis:: Initialized as observer";
77 edm::LogVerbatim(
"ForwardSim") <<
"============================================================================";
82 new TNtuple(
"NTcastorstep",
"NTcastorstep",
"evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz");
86 "NTcastorevent",
"NTcastorevent",
"evt:ihit:detector:sector:module:enem:totalenergy:hitenergy:x:y:z");
105 edm::LogVerbatim(
"ForwardSim") << std::endl <<
"CastorTestAnalysis: Starting Run";
107 edm::LogVerbatim(
"ForwardSim") <<
"CastorTestAnalysis: output step root file created";
113 edm::LogVerbatim(
"ForwardSim") <<
"CastorTestAnalysis: output event root file created";
131 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
133 G4double stepL = aStep->GetStepLength();
134 G4double stepE = aStep->GetTotalEnergyDeposit();
139 G4Track* theTrack = aStep->GetTrack();
140 G4int theTrackID = theTrack->GetTrackID();
141 G4double
theCharge = theTrack->GetDynamicParticle()->GetCharge();
143 G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
145 const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
146 G4double vpx = vert_mom.x();
147 G4double vpy = vert_mom.y();
148 G4double vpz = vert_mom.z();
149 double eta = 0.5 *
log((1. + vpz) / (1. - vpz));
150 double phi = atan2(vpy, vpx);
152 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
197 if (theTrack->GetNextVolume() != 0)
198 edm::LogVerbatim(
"ForwardSim") <<
" NextVolume: " << theTrack->GetNextVolume()->GetName();
212 <<
"CastorTest::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
215 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
218 int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID(
"CastorFI");
242 if (theCAFI->entries() > 0)
250 G4PrimaryParticle* thePrim =
nullptr;
251 G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
254 edm::LogVerbatim(
"ForwardSim") <<
"CASTORTest End Of Event ERROR: no vertex";
256 for (
int i = 0;
i < nvertex;
i++) {
257 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(
i);
258 if (avertex ==
nullptr)
259 edm::LogVerbatim(
"ForwardSim") <<
"CASTORTest End Of Event ERR: pointer to vertex = 0";
261 int npart = avertex->GetNumberOfParticle();
263 edm::LogVerbatim(
"ForwardSim") <<
"CASTORTest End Of Event ERR: no primary!";
264 if (thePrim ==
nullptr)
265 thePrim = avertex->GetPrimary(trackID);
268 double px = 0.,
py = 0., pz = 0., pInit = 0;
270 double eta = 0.,
phi = 0.;
272 if (thePrim !=
nullptr) {
273 px = thePrim->GetPx();
274 py = thePrim->GetPy();
275 pz = thePrim->GetPz();
278 edm::LogVerbatim(
"ForwardSim") <<
"CASTORTest End Of Event ERR: primary has p=0 ";
281 float costheta = pz / pInit;
293 edm::LogVerbatim(
"ForwardSim") <<
"CASTORTest End Of Event ERR: could not find primary ";
297 <<
", " <<
eta <<
", " <<
phi;
301 int iEvt = (*evt)()->GetEventID();
304 else if ((iEvt < 100) && (iEvt % 10 == 0))
306 else if ((iEvt < 1000) && (iEvt % 100 == 0))
308 else if ((iEvt < 10000) && (iEvt % 1000 == 0))
311 edm::LogVerbatim(
"ForwardSim") << std::endl <<
"===>>> Done writing user histograms ";
318 int nentries = hc->entries();
321 unsigned int volumeID = 0;
322 int det = 0,
zside, sector, zmodule;
323 std::map<int, float, std::less<int> > themap;
324 double totalEnergy = 0;
325 double hitEnergy = 0;
326 double en_in_sd = 0.;
328 for (
int ihit = 0; ihit < nentries; ihit++) {
333 for (
int ihit = 0; ihit < nentries; ihit++) {
376 edm::LogVerbatim(
"ForwardSim") <<
"CastorTestAnalysis: Ntuple step written";
384 edm::LogVerbatim(
"ForwardSim") <<
"CastorTestAnalysis: Ntuple event written";