File indexing completed on 2023-03-17 11:24:25
0001
0002
0003
0004
0005
0006
0007 #include "DataFormats/Math/interface/Point3D.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0012 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0013 #include "SimG4CMS/Forward/interface/ZdcNumberingScheme.h"
0014
0015 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0016 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0017 #include "SimG4Core/Notification/interface/EndOfRun.h"
0018 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0019 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0020 #include "SimG4Core/Notification/interface/Observer.h"
0021 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0022
0023 #include "G4SDManager.hh"
0024 #include "G4Step.hh"
0025 #include "G4Track.hh"
0026 #include "G4Event.hh"
0027 #include "G4PrimaryVertex.hh"
0028 #include "G4VProcess.hh"
0029 #include "G4HCofThisEvent.hh"
0030 #include "G4UserEventAction.hh"
0031
0032 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0033 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0034 #include <CLHEP/Random/Randomize.h>
0035
0036 #include "TROOT.h"
0037 #include "TH1.h"
0038 #include "TH2.h"
0039 #include "TProfile.h"
0040 #include "TNtuple.h"
0041 #include "TRandom.h"
0042 #include "TLorentzVector.h"
0043 #include "TUnixSystem.h"
0044 #include "TSystem.h"
0045 #include "TMath.h"
0046 #include "TF1.h"
0047 #include "TFile.h"
0048
0049 #include <cassert>
0050 #include <cmath>
0051 #include <iostream>
0052 #include <iomanip>
0053 #include <map>
0054 #include <string>
0055 #include <vector>
0056
0057 class ZdcTestAnalysis : public SimWatcher,
0058 public Observer<const BeginOfJob*>,
0059 public Observer<const BeginOfRun*>,
0060 public Observer<const EndOfRun*>,
0061 public Observer<const BeginOfEvent*>,
0062 public Observer<const EndOfEvent*>,
0063 public Observer<const G4Step*> {
0064 public:
0065 ZdcTestAnalysis(const edm::ParameterSet& p);
0066 ~ZdcTestAnalysis() override;
0067
0068 private:
0069
0070 void update(const BeginOfJob* run) override;
0071 void update(const BeginOfRun* run) override;
0072 void update(const EndOfRun* run) override;
0073 void update(const BeginOfEvent* evt) override;
0074 void update(const EndOfEvent* evt) override;
0075 void update(const G4Step* step) override;
0076
0077 void finish();
0078
0079 int verbosity;
0080 int doNTzdcstep;
0081 int doNTzdcevent;
0082 std::string stepNtFileName;
0083 std::string eventNtFileName;
0084
0085 TFile* zdcOutputEventFile;
0086 TFile* zdcOutputStepFile;
0087
0088 TNtuple* zdcstepntuple;
0089 TNtuple* zdceventntuple;
0090
0091 int eventIndex;
0092 int stepIndex;
0093
0094 Float_t zdcsteparray[18];
0095 Float_t zdceventarray[16];
0096
0097 ZdcNumberingScheme* theZdcNumScheme;
0098 };
0099
0100 enum ntzdcs_elements {
0101 ntzdcs_evt,
0102 ntzdcs_trackid,
0103 ntzdcs_charge,
0104 ntzdcs_pdgcode,
0105 ntzdcs_x,
0106 ntzdcs_y,
0107 ntzdcs_z,
0108 ntzdcs_stepl,
0109 ntzdcs_stepe,
0110 ntzdcs_eta,
0111 ntzdcs_phi,
0112 ntzdcs_vpx,
0113 ntzdcs_vpy,
0114 ntzdcs_vpz,
0115 ntzdcs_idx,
0116 ntzdcs_idl,
0117 ntzdcs_pvtype,
0118 ntzdcs_ncherphot
0119 };
0120
0121 enum ntzdce_elements {
0122 ntzdce_evt,
0123 ntzdce_ihit,
0124 ntzdce_fiberid,
0125 ntzdce_zside,
0126 ntzdce_subdet,
0127 ntzdce_layer,
0128 ntzdce_fiber,
0129 ntzdce_channel,
0130 ntzdce_enem,
0131 ntzdce_enhad,
0132 ntzdce_hitenergy,
0133 ntzdce_x,
0134 ntzdce_y,
0135 ntzdce_z,
0136 ntzdce_time,
0137 ntzdce_etot
0138 };
0139
0140 ZdcTestAnalysis::ZdcTestAnalysis(const edm::ParameterSet& p) {
0141
0142 edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("ZdcTestAnalysis");
0143 verbosity = m_Anal.getParameter<int>("Verbosity");
0144 doNTzdcstep = m_Anal.getParameter<int>("StepNtupleFlag");
0145 doNTzdcevent = m_Anal.getParameter<int>("EventNtupleFlag");
0146 stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
0147 eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
0148
0149 if (verbosity > 0)
0150 edm::LogVerbatim("ZdcAnalysis") << "\n============================================================================";
0151 edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis:: Initialized as observer";
0152 if (doNTzdcstep > 0) {
0153 edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple will be created";
0154 edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple file: " << stepNtFileName;
0155 }
0156 if (doNTzdcevent > 0) {
0157 edm::LogVerbatim("ZdcAnalysis") << " Event Ntuple will be created";
0158 edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple file: " << stepNtFileName;
0159 }
0160 edm::LogVerbatim("ZdcAnalysis") << "============================================================================"
0161 << std::endl;
0162
0163 if (doNTzdcstep > 0)
0164 zdcstepntuple =
0165 new TNtuple("NTzdcstep",
0166 "NTzdcstep",
0167 "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz:idx:idl:pvtype:ncherphot");
0168
0169 if (doNTzdcevent > 0)
0170 zdceventntuple =
0171 new TNtuple("NTzdcevent",
0172 "NTzdcevent",
0173 "evt:ihit:fiberid:zside:subdet:layer:fiber:channel:enem:enhad:hitenergy:x:y:z:time:etot");
0174
0175 theZdcNumScheme = nullptr;
0176
0177 }
0178
0179 ZdcTestAnalysis::~ZdcTestAnalysis() {
0180
0181 finish();
0182 delete theZdcNumScheme;
0183 }
0184
0185 void ZdcTestAnalysis::update(const BeginOfJob* job) {
0186
0187 edm::LogVerbatim("ZdcAnalysis") << "beggining of job";
0188 }
0189
0190
0191 void ZdcTestAnalysis::update(const BeginOfRun* run) {
0192
0193
0194 edm::LogVerbatim("ZdcAnalysis") << "\nZdcTestAnalysis: Begining of Run";
0195 if (doNTzdcstep) {
0196 edm::LogVerbatim("ZdcAnalysis") << "ZDCTestAnalysis: output step file created";
0197 TString stepfilename = stepNtFileName;
0198 zdcOutputStepFile = new TFile(stepfilename, "RECREATE");
0199 }
0200
0201 if (doNTzdcevent) {
0202 edm::LogVerbatim("ZdcAnalysis") << "ZDCTestAnalysis: output event file created";
0203 TString stepfilename = eventNtFileName;
0204 zdcOutputEventFile = new TFile(stepfilename, "RECREATE");
0205 }
0206
0207 eventIndex = 0;
0208 }
0209
0210 void ZdcTestAnalysis::update(const BeginOfEvent* evt) {
0211
0212 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest: Processing Event Number: " << eventIndex;
0213 eventIndex++;
0214 stepIndex = 0;
0215 }
0216
0217 void ZdcTestAnalysis::update(const G4Step* aStep) {
0218
0219 stepIndex++;
0220
0221 if (doNTzdcstep) {
0222 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0223
0224 G4double stepL = aStep->GetStepLength();
0225 G4double stepE = aStep->GetTotalEnergyDeposit();
0226
0227 if (verbosity >= 2)
0228 edm::LogVerbatim("ZdcAnalysis") << "Step " << stepL << "," << stepE;
0229
0230 G4Track* theTrack = aStep->GetTrack();
0231 G4int theTrackID = theTrack->GetTrackID();
0232 G4double theCharge = theTrack->GetDynamicParticle()->GetCharge();
0233 G4String particleType = theTrack->GetDefinition()->GetParticleName();
0234 G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
0235
0236 const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
0237 G4double vpx = vert_mom.x();
0238 G4double vpy = vert_mom.y();
0239 G4double vpz = vert_mom.z();
0240 double eta = 0.5 * log((1. + vpz) / (1. - vpz));
0241 double phi = atan2(vpy, vpx);
0242
0243 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0244 G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0245
0246 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0247 int idx = touch->GetReplicaNumber(0);
0248 int idLayer = -1;
0249 int thePVtype = -1;
0250
0251 int historyDepth = touch->GetHistoryDepth();
0252
0253 if (historyDepth > 0) {
0254 std::vector<int> theReplicaNumbers(historyDepth);
0255 std::vector<G4VPhysicalVolume*> thePhysicalVolumes(historyDepth);
0256 std::vector<G4String> thePVnames(historyDepth);
0257 std::vector<G4LogicalVolume*> theLogicalVolumes(historyDepth);
0258 std::vector<G4String> theLVnames(historyDepth);
0259 std::vector<G4Material*> theMaterials(historyDepth);
0260 std::vector<G4String> theMaterialNames(historyDepth);
0261
0262 for (int jj = 0; jj < historyDepth; jj++) {
0263 theReplicaNumbers[jj] = touch->GetReplicaNumber(jj);
0264 thePhysicalVolumes[jj] = touch->GetVolume(jj);
0265 thePVnames[jj] = thePhysicalVolumes[jj]->GetName();
0266 theLogicalVolumes[jj] = thePhysicalVolumes[jj]->GetLogicalVolume();
0267 theLVnames[jj] = theLogicalVolumes[jj]->GetName();
0268 theMaterials[jj] = theLogicalVolumes[jj]->GetMaterial();
0269 theMaterialNames[jj] = theMaterials[jj]->GetName();
0270 if (verbosity >= 2)
0271 edm::LogVerbatim("ZdcAnalysis") << " GHD " << jj << ": " << theReplicaNumbers[jj] << "," << thePVnames[jj]
0272 << "," << theLVnames[jj] << "," << theMaterialNames[jj];
0273 }
0274
0275 idLayer = theReplicaNumbers[1];
0276 if (thePVnames[0] == "ZDC_EMLayer")
0277 thePVtype = 1;
0278 else if (thePVnames[0] == "ZDC_EMAbsorber")
0279 thePVtype = 2;
0280 else if (thePVnames[0] == "ZDC_EMFiber")
0281 thePVtype = 3;
0282 else if (thePVnames[0] == "ZDC_HadLayer")
0283 thePVtype = 7;
0284 else if (thePVnames[0] == "ZDC_HadAbsorber")
0285 thePVtype = 8;
0286 else if (thePVnames[0] == "ZDC_HadFiber")
0287 thePVtype = 9;
0288 else if (thePVnames[0] == "ZDC_PhobosLayer")
0289 thePVtype = 11;
0290 else if (thePVnames[0] == "ZDC_PhobosAbsorber")
0291 thePVtype = 12;
0292 else if (thePVnames[0] == "ZDC_PhobosFiber")
0293 thePVtype = 13;
0294 else {
0295 thePVtype = 0;
0296 if (verbosity >= 2)
0297 edm::LogVerbatim("ZdcAnalysis") << " pvtype=0 hd=" << historyDepth << " " << theReplicaNumbers[0] << ","
0298 << thePVnames[0] << "," << theLVnames[0] << "," << theMaterialNames[0];
0299 }
0300 } else if (historyDepth == 0) {
0301 int theReplicaNumber = touch->GetReplicaNumber(0);
0302 G4VPhysicalVolume* thePhysicalVolume = touch->GetVolume(0);
0303 const G4String& thePVname = thePhysicalVolume->GetName();
0304 G4LogicalVolume* theLogicalVolume = thePhysicalVolume->GetLogicalVolume();
0305 const G4String& theLVname = theLogicalVolume->GetName();
0306 G4Material* theMaterial = theLogicalVolume->GetMaterial();
0307 const G4String& theMaterialName = theMaterial->GetName();
0308 if (verbosity >= 2)
0309 edm::LogVerbatim("ZdcAnalysis") << " hd=0 " << theReplicaNumber << "," << thePVname << "," << theLVname << ","
0310 << theMaterialName;
0311 } else {
0312 edm::LogVerbatim("ZdcAnalysis") << " hd<0: hd=" << historyDepth;
0313 }
0314
0315 double NCherPhot = -1;
0316 zdcsteparray[ntzdcs_evt] = (float)eventIndex;
0317 zdcsteparray[ntzdcs_trackid] = (float)theTrackID;
0318 zdcsteparray[ntzdcs_charge] = theCharge;
0319 zdcsteparray[ntzdcs_pdgcode] = (float)pdgcode;
0320 zdcsteparray[ntzdcs_x] = hitPoint.x();
0321 zdcsteparray[ntzdcs_y] = hitPoint.y();
0322 zdcsteparray[ntzdcs_z] = hitPoint.z();
0323 zdcsteparray[ntzdcs_stepl] = stepL;
0324 zdcsteparray[ntzdcs_stepe] = stepE;
0325 zdcsteparray[ntzdcs_eta] = eta;
0326 zdcsteparray[ntzdcs_phi] = phi;
0327 zdcsteparray[ntzdcs_vpx] = vpx;
0328 zdcsteparray[ntzdcs_vpy] = vpy;
0329 zdcsteparray[ntzdcs_vpz] = vpz;
0330 zdcsteparray[ntzdcs_idx] = (float)idx;
0331 zdcsteparray[ntzdcs_idl] = (float)idLayer;
0332 zdcsteparray[ntzdcs_pvtype] = thePVtype;
0333 zdcsteparray[ntzdcs_ncherphot] = NCherPhot;
0334 zdcstepntuple->Fill(zdcsteparray);
0335 }
0336 }
0337
0338 void ZdcTestAnalysis::update(const EndOfEvent* evt) {
0339
0340
0341
0342 edm::LogVerbatim("ZdcAnalysis") << "\nZdcTest::upDate(const EndOfEvent * evt) - event #" << (*evt)()->GetEventID()
0343 << "\n # of aSteps followed in event = " << stepIndex;
0344
0345
0346 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0347 edm::LogVerbatim("ZdcAnalysis") << " accessed all HC";
0348
0349 int theZDCHCid = G4SDManager::GetSDMpointer()->GetCollectionID("ZDCHITS");
0350 edm::LogVerbatim("ZdcAnalysis") << " - theZDCHCid = " << theZDCHCid;
0351
0352 CaloG4HitCollection* theZDCHC = (CaloG4HitCollection*)allHC->GetHC(theZDCHCid);
0353 edm::LogVerbatim("ZdcAnalysis") << " - theZDCHC = " << theZDCHC;
0354
0355 if (!theZdcNumScheme) {
0356 theZdcNumScheme = new ZdcNumberingScheme(1);
0357 }
0358
0359
0360 int maxTime = 0;
0361 int fiberID = 0;
0362 unsigned int unsignedfiberID = 0;
0363 std::map<int, float, std::less<int> > energyInFibers;
0364 std::map<int, float, std::less<int> > primaries;
0365 float totalEnergy = 0;
0366 int nentries = theZDCHC->entries();
0367 edm::LogVerbatim("ZdcAnalysis") << " theZDCHC has " << nentries << " entries";
0368
0369 if (doNTzdcevent) {
0370 if (nentries > 0) {
0371 for (int ihit = 0; ihit < nentries; ihit++) {
0372 CaloG4Hit* caloHit = (*theZDCHC)[ihit];
0373 totalEnergy += caloHit->getEnergyDeposit();
0374 }
0375
0376 for (int ihit = 0; ihit < nentries; ihit++) {
0377 CaloG4Hit* aHit = (*theZDCHC)[ihit];
0378 fiberID = aHit->getUnitID();
0379 unsignedfiberID = aHit->getUnitID();
0380 double enEm = aHit->getEM();
0381 double enHad = aHit->getHadr();
0382 math::XYZPoint hitPoint = aHit->getPosition();
0383 double hitEnergy = aHit->getEnergyDeposit();
0384 if (verbosity >= 1)
0385 edm::LogVerbatim("ZdcAnalysis")
0386 << " entry #" << ihit << ": fiberID=0x" << std::hex << fiberID << std::dec << "; enEm=" << enEm
0387 << "; enHad=" << enHad << "; hitEnergy=" << hitEnergy << "z=" << hitPoint.z();
0388 energyInFibers[fiberID] += enEm + enHad;
0389 primaries[aHit->getTrackID()] += enEm + enHad;
0390 float time = aHit->getTimeSliceID();
0391 if (time > maxTime)
0392 maxTime = (int)time;
0393
0394 int thesubdet, thelayer, thefiber, thechannel, thez;
0395 theZdcNumScheme->unpackZdcIndex(fiberID, thesubdet, thelayer, thefiber, thechannel, thez);
0396 int unsignedsubdet, unsignedlayer, unsignedfiber, unsignedchannel, unsignedz;
0397 theZdcNumScheme->unpackZdcIndex(
0398 unsignedfiberID, unsignedsubdet, unsignedlayer, unsignedfiber, unsignedchannel, unsignedz);
0399
0400
0401
0402
0403
0404
0405 zdceventarray[ntzdce_evt] = (float)eventIndex;
0406 zdceventarray[ntzdce_ihit] = (float)ihit;
0407 zdceventarray[ntzdce_fiberid] = (float)fiberID;
0408 zdceventarray[ntzdce_zside] = (float)thez;
0409 zdceventarray[ntzdce_subdet] = (float)thesubdet;
0410 zdceventarray[ntzdce_layer] = (float)thelayer;
0411 zdceventarray[ntzdce_fiber] = (float)thefiber;
0412 zdceventarray[ntzdce_channel] = (float)thechannel;
0413 zdceventarray[ntzdce_enem] = enEm;
0414 zdceventarray[ntzdce_enhad] = enHad;
0415 zdceventarray[ntzdce_hitenergy] = hitEnergy;
0416 zdceventarray[ntzdce_x] = hitPoint.x();
0417 zdceventarray[ntzdce_y] = hitPoint.y();
0418 zdceventarray[ntzdce_z] = hitPoint.z();
0419 zdceventarray[ntzdce_time] = time;
0420 zdceventarray[ntzdce_etot] = totalEnergy;
0421 zdceventntuple->Fill(zdceventarray);
0422 }
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432 int trackID = 0;
0433 G4PrimaryParticle* thePrim = nullptr;
0434 G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0435 edm::LogVerbatim("ZdcAnalysis") << "Event has " << nvertex << " vertex";
0436 if (nvertex == 0)
0437 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERROR: no vertex";
0438
0439 for (int i = 0; i < nvertex; i++) {
0440 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0441 if (avertex == nullptr) {
0442 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: pointer to vertex = 0";
0443 } else {
0444 edm::LogVerbatim("ZdcAnalysis") << "Vertex number :" << i;
0445 int npart = avertex->GetNumberOfParticle();
0446 if (npart == 0)
0447 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: no primary!";
0448 if (thePrim == nullptr)
0449 thePrim = avertex->GetPrimary(trackID);
0450 }
0451 }
0452
0453 double px = 0., py = 0., pz = 0.;
0454 double pInit = 0.;
0455
0456 if (thePrim != nullptr) {
0457 px = thePrim->GetPx();
0458 py = thePrim->GetPy();
0459 pz = thePrim->GetPz();
0460 pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0461 if (pInit == 0) {
0462 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: primary has p=0 ";
0463 }
0464 } else {
0465 edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: could not find primary ";
0466 }
0467
0468 }
0469
0470 }
0471
0472 int iEvt = (*evt)()->GetEventID();
0473 if (iEvt < 10)
0474 edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0475 else if ((iEvt < 100) && (iEvt % 10 == 0))
0476 edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0477 else if ((iEvt < 1000) && (iEvt % 100 == 0))
0478 edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0479 else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0480 edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0481 }
0482
0483 void ZdcTestAnalysis::update(const EndOfRun* run) {}
0484
0485 void ZdcTestAnalysis::finish() {
0486 if (doNTzdcstep) {
0487 zdcOutputStepFile->cd();
0488 zdcstepntuple->Write();
0489 edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple step written for event: " << eventIndex;
0490 zdcOutputStepFile->Close();
0491 edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Step file closed";
0492 }
0493
0494 if (doNTzdcevent) {
0495 zdcOutputEventFile->cd();
0496 zdceventntuple->Write("", TObject::kOverwrite);
0497 edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple event written for event: " << eventIndex;
0498 zdcOutputEventFile->Close();
0499 edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Event file closed";
0500 }
0501 }
0502
0503 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0504 #include "FWCore/PluginManager/interface/ModuleDef.h"
0505
0506 DEFINE_SIMWATCHER(ZdcTestAnalysis);