Forward/plugins/ZdcTestAnalysis.cc

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