Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-13 23:27:40

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: ZdcTestAnalysis.cc
0003 // Date: 03.06 Edmundo Garcia
0004 // Description: simulation analysis steering code
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   // observer classes
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   //constructor
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   //theZdcSD = new ZdcSD("ZDCHITSB", new ZdcNumberingScheme());
0177 }
0178 
0179 ZdcTestAnalysis::~ZdcTestAnalysis() {
0180   // destructor
0181   finish();
0182   delete theZdcNumScheme;
0183 }
0184 
0185 void ZdcTestAnalysis::update(const BeginOfJob* job) {
0186   //job
0187   edm::LogVerbatim("ZdcAnalysis") << "beggining of job";
0188 }
0189 
0190 //==================================================================== per RUN
0191 void ZdcTestAnalysis::update(const BeginOfRun* run) {
0192   //run
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   //event
0212   edm::LogVerbatim("ZdcAnalysis") << "ZdcTest: Processing Event Number: " << eventIndex;
0213   eventIndex++;
0214   stepIndex = 0;
0215 }
0216 
0217 void ZdcTestAnalysis::update(const G4Step* aStep) {
0218   //step;
0219   stepIndex++;
0220 
0221   if (doNTzdcstep) {
0222     G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0223     // G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
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   //end of event
0340 
0341   // Look for the Hit Collection
0342   edm::LogVerbatim("ZdcAnalysis") << "\nZdcTest::upDate(const EndOfEvent * evt) - event #" << (*evt)()->GetEventID()
0343                                   << "\n  # of aSteps followed in event = " << stepIndex;
0344 
0345   // access to the G4 hit collections
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   float ETot = 0., SEnergy = 0.;
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         // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
0401         // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
0402         // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
0403         // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
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       for (std::map<int, float, std::less<int> >::iterator is = energyInFibers.begin(); is != energyInFibers.end();
0425            is++) {
0426         ETot = (*is).second;
0427         SEnergy += ETot;
0428       }
0429 
0430       // Find Primary info:
0431       int trackID = 0;
0432       G4PrimaryParticle* thePrim = nullptr;
0433       G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0434       edm::LogVerbatim("ZdcAnalysis") << "Event has " << nvertex << " vertex";
0435       if (nvertex == 0)
0436         edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event  ERROR: no vertex";
0437 
0438       for (int i = 0; i < nvertex; i++) {
0439         G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0440         if (avertex == nullptr) {
0441           edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: pointer to vertex = 0";
0442         } else {
0443           edm::LogVerbatim("ZdcAnalysis") << "Vertex number :" << i;
0444           int npart = avertex->GetNumberOfParticle();
0445           if (npart == 0)
0446             edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: no primary!";
0447           if (thePrim == nullptr)
0448             thePrim = avertex->GetPrimary(trackID);
0449         }
0450       }
0451 
0452       double px = 0., py = 0., pz = 0.;
0453       double pInit = 0.;
0454 
0455       if (thePrim != nullptr) {
0456         px = thePrim->GetPx();
0457         py = thePrim->GetPy();
0458         pz = thePrim->GetPz();
0459         pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0460         if (pInit == 0) {
0461           edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event  ERR: primary has p=0 ";
0462         }
0463       } else {
0464         edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: could not find primary ";
0465       }
0466 
0467     }  // nentries > 0
0468 
0469   }  // createNTzdcevent
0470 
0471   int iEvt = (*evt)()->GetEventID();
0472   if (iEvt < 10)
0473     edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0474   else if ((iEvt < 100) && (iEvt % 10 == 0))
0475     edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0476   else if ((iEvt < 1000) && (iEvt % 100 == 0))
0477     edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0478   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0479     edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
0480 }
0481 
0482 void ZdcTestAnalysis::update(const EndOfRun* run) {}
0483 
0484 void ZdcTestAnalysis::finish() {
0485   if (doNTzdcstep) {
0486     zdcOutputStepFile->cd();
0487     zdcstepntuple->Write();
0488     edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple step  written for event: " << eventIndex;
0489     zdcOutputStepFile->Close();
0490     edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Step file closed";
0491   }
0492 
0493   if (doNTzdcevent) {
0494     zdcOutputEventFile->cd();
0495     zdceventntuple->Write("", TObject::kOverwrite);
0496     edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple event written for event: " << eventIndex;
0497     zdcOutputEventFile->Close();
0498     edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Event file closed";
0499   }
0500 }
0501 
0502 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0503 #include "FWCore/PluginManager/interface/ModuleDef.h"
0504 
0505 DEFINE_SIMWATCHER(ZdcTestAnalysis);