Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:18 

0001 // -*- C++ -*-
0002 //
0003 // Package:     Forward
0004 // Class  :     CastorTestAnalysis
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author: P. Katsas
0010 //         Created: 02/2007
0011 //
0012 
0013 #include "SimG4Core/Notification/interface/Observer.h"
0014 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0015 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0016 #include "SimG4Core/Notification/interface/EndOfRun.h"
0017 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0018 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0019 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0020 
0021 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0022 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0023 #include "SimG4CMS/Forward/interface/CastorNumberingScheme.h"
0024 
0025 #include "DataFormats/Math/interface/Point3D.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 
0029 #include "G4SDManager.hh"
0030 #include "G4Step.hh"
0031 #include "G4Track.hh"
0032 #include "G4Event.hh"
0033 #include "G4PrimaryVertex.hh"
0034 #include "G4VProcess.hh"
0035 #include "G4HCofThisEvent.hh"
0036 #include "G4UserEventAction.hh"
0037 
0038 #include <CLHEP/Units/SystemOfUnits.h>
0039 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0040 #include <CLHEP/Random/Randomize.h>
0041 
0042 #include "TROOT.h"
0043 #include "TFile.h"
0044 #include "TH1.h"
0045 #include "TH2.h"
0046 #include "TProfile.h"
0047 #include "TNtuple.h"
0048 #include "TRandom.h"
0049 #include "TLorentzVector.h"
0050 #include "TUnixSystem.h"
0051 #include "TSystem.h"
0052 #include "TMath.h"
0053 #include "TF1.h"
0054 
0055 #include <cassert>
0056 #include <cmath>
0057 #include <iostream>
0058 #include <iomanip>
0059 #include <map>
0060 #include <string>
0061 #include <vector>
0062 
0063 //#define EDM_ML_DEBUG
0064 
0065 class CastorTestAnalysis : public SimWatcher,
0066                            public Observer<const BeginOfJob *>,
0067                            public Observer<const BeginOfRun *>,
0068                            public Observer<const EndOfRun *>,
0069                            public Observer<const BeginOfEvent *>,
0070                            public Observer<const EndOfEvent *>,
0071                            public Observer<const G4Step *> {
0072 public:
0073   CastorTestAnalysis(const edm::ParameterSet &p);
0074   ~CastorTestAnalysis() override;
0075 
0076 private:
0077   // observer classes
0078   void update(const BeginOfJob *run) override;
0079   void update(const BeginOfRun *run) override;
0080   void update(const EndOfRun *run) override;
0081   void update(const BeginOfEvent *evt) override;
0082   void update(const EndOfEvent *evt) override;
0083   void update(const G4Step *step) override;
0084 
0085 private:
0086   void getCastorBranchData(const CaloG4HitCollection *hc);
0087   void Finish();
0088 
0089   int verbosity;
0090   int doNTcastorstep;
0091   int doNTcastorevent;
0092   std::string stepNtFileName;
0093   std::string eventNtFileName;
0094 
0095   TFile *castorOutputEventFile;
0096   TFile *castorOutputStepFile;
0097 
0098   TNtuple *castorstepntuple;
0099   TNtuple *castoreventntuple;
0100 
0101   CastorNumberingScheme *theCastorNumScheme;
0102 
0103   int eventIndex;
0104   int stepIndex;
0105   int eventGlobalHit;
0106 
0107   Float_t castorsteparray[14];
0108   Float_t castoreventarray[11];
0109 };
0110 
0111 enum ntcastors_elements {
0112   ntcastors_evt,
0113   ntcastors_trackid,
0114   ntcastors_charge,
0115   ntcastors_pdgcode,
0116   ntcastors_x,
0117   ntcastors_y,
0118   ntcastors_z,
0119   ntcastors_stepl,
0120   ntcastors_stepe,
0121   ntcastors_eta,
0122   ntcastors_phi,
0123   ntcastors_vpx,
0124   ntcastors_vpy,
0125   ntcastors_vpz
0126 };
0127 
0128 enum ntcastore_elements {
0129   ntcastore_evt,
0130   ntcastore_ihit,
0131   ntcastore_detector,
0132   ntcastore_sector,
0133   ntcastore_module,
0134   ntcastore_enem,
0135   ntcastore_enhad,
0136   ntcastore_hitenergy,
0137   ntcastore_x,
0138   ntcastore_y,
0139   ntcastore_z
0140 };
0141 
0142 CastorTestAnalysis::CastorTestAnalysis(const edm::ParameterSet &p) {
0143   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("CastorTestAnalysis");
0144   verbosity = m_Anal.getParameter<int>("Verbosity");
0145   doNTcastorstep = m_Anal.getParameter<int>("StepNtupleFlag");
0146   doNTcastorevent = m_Anal.getParameter<int>("EventNtupleFlag");
0147   stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
0148   eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
0149 
0150   if (verbosity > 0) {
0151     edm::LogVerbatim("ForwardSim") << std::endl;
0152     edm::LogVerbatim("ForwardSim") << "============================================================================";
0153     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis:: Initialized as observer";
0154     if (doNTcastorstep > 0) {
0155       edm::LogVerbatim("ForwardSim") << " Step Ntuple will be created";
0156       edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
0157     }
0158     if (doNTcastorevent > 0) {
0159       edm::LogVerbatim("ForwardSim") << " Event Ntuple will be created";
0160       edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
0161     }
0162     edm::LogVerbatim("ForwardSim") << "============================================================================";
0163     edm::LogVerbatim("ForwardSim") << std::endl;
0164   }
0165   if (doNTcastorstep > 0)
0166     castorstepntuple =
0167         new TNtuple("NTcastorstep", "NTcastorstep", "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz");
0168 
0169   if (doNTcastorevent > 0)
0170     castoreventntuple = new TNtuple(
0171         "NTcastorevent", "NTcastorevent", "evt:ihit:detector:sector:module:enem:totalenergy:hitenergy:x:y:z");
0172 }
0173 
0174 CastorTestAnalysis::~CastorTestAnalysis() {
0175   //destructor of CastorTestAnalysis
0176 
0177   Finish();
0178   if (verbosity > 0) {
0179     edm::LogVerbatim("ForwardSim") << std::endl << "End of CastorTestAnalysis";
0180   }
0181 
0182   edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: End of process";
0183 }
0184 
0185 //=================================================================== per EVENT
0186 void CastorTestAnalysis::update(const BeginOfJob *job) { edm::LogVerbatim("ForwardSim") << " Starting new job "; }
0187 
0188 //==================================================================== per RUN
0189 void CastorTestAnalysis::update(const BeginOfRun *run) {
0190   edm::LogVerbatim("ForwardSim") << std::endl << "CastorTestAnalysis: Starting Run";
0191   if (doNTcastorstep) {
0192     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output step root file created";
0193     TString stepfilename = stepNtFileName;
0194     castorOutputStepFile = new TFile(stepfilename, "RECREATE");
0195   }
0196 
0197   if (doNTcastorevent) {
0198     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output event root file created";
0199     TString stepfilename = eventNtFileName;
0200     castorOutputEventFile = new TFile(stepfilename, "RECREATE");
0201   }
0202 
0203   eventIndex = 0;
0204 }
0205 
0206 void CastorTestAnalysis::update(const BeginOfEvent *evt) {
0207   edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Processing Event Number: " << eventIndex;
0208   eventIndex++;
0209   stepIndex = 0;
0210 }
0211 
0212 void CastorTestAnalysis::update(const G4Step *aStep) {
0213   stepIndex++;
0214 
0215   if (doNTcastorstep) {
0216     G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0217     //  G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
0218     G4double stepL = aStep->GetStepLength();
0219     G4double stepE = aStep->GetTotalEnergyDeposit();
0220 
0221     if (verbosity >= 2)
0222       edm::LogVerbatim("ForwardSim") << "Step " << stepL << ", " << stepE;
0223 
0224     G4Track *theTrack = aStep->GetTrack();
0225     G4int theTrackID = theTrack->GetTrackID();
0226     G4double theCharge = theTrack->GetDynamicParticle()->GetCharge();
0227     //  G4String particleType = theTrack->GetDefinition()->GetParticleName();
0228     G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
0229 
0230     const G4ThreeVector &vert_mom = theTrack->GetVertexMomentumDirection();
0231     G4double vpx = vert_mom.x();
0232     G4double vpy = vert_mom.y();
0233     G4double vpz = vert_mom.z();
0234     double eta = 0.5 * log((1. + vpz) / (1. - vpz));
0235     double phi = atan2(vpy, vpx);
0236 
0237     const G4ThreeVector &hitPoint = preStepPoint->GetPosition();
0238 
0239     // Fill-in ntuple
0240     //  castorsteparray[ntcastors_evt] = (*evt)()->GetEventID();
0241     castorsteparray[ntcastors_evt] = (float)eventIndex;
0242     castorsteparray[ntcastors_trackid] = (float)theTrackID;
0243     castorsteparray[ntcastors_charge] = theCharge;
0244     castorsteparray[ntcastors_pdgcode] = pdgcode;
0245     castorsteparray[ntcastors_x] = hitPoint.x();
0246     castorsteparray[ntcastors_y] = hitPoint.y();
0247     castorsteparray[ntcastors_z] = hitPoint.z();
0248     castorsteparray[ntcastors_stepl] = stepL;
0249     castorsteparray[ntcastors_stepe] = stepE;
0250     castorsteparray[ntcastors_eta] = eta;
0251     castorsteparray[ntcastors_phi] = phi;
0252     castorsteparray[ntcastors_vpx] = vpx;
0253     castorsteparray[ntcastors_vpy] = vpy;
0254     castorsteparray[ntcastors_vpz] = vpz;
0255 
0256     /*
0257   edm::LogVerbatim("ForwardSim") << "TrackID: " << theTrackID;
0258   edm::LogVerbatim("ForwardSim") << "   StepN: "<< theTrack->GetCurrentStepNumber();
0259   edm::LogVerbatim("ForwardSim") << "      ParentID: " << aStep->GetTrack()->GetParentID();
0260   edm::LogVerbatim("ForwardSim") << "      PDG: " << pdgcode;
0261   edm::LogVerbatim("ForwardSim") << "      X,Y,Z (mm): " << theTrack->GetPosition().x() << "," << theTrack->GetPosition().y() << "," << theTrack->GetPosition().z();
0262   edm::LogVerbatim("ForwardSim") << "      KE (MeV): " << theTrack->GetKineticEnergy();
0263   edm::LogVerbatim("ForwardSim") << "      Total EDep (MeV): " << aStep->GetTotalEnergyDeposit();
0264   edm::LogVerbatim("ForwardSim") << "      StepLength (mm): " << aStep->GetStepLength();
0265   edm::LogVerbatim("ForwardSim") << "      TrackLength (mm): " << theTrack->GetTrackLength();
0266 
0267   if ( theTrack->GetNextVolume() != 0 )
0268       edm::LogVerbatim("ForwardSim") <<"      NextVolume: " << theTrack->GetNextVolume()->GetName();
0269   else 
0270       edm::LogVerbatim("ForwardSim") <<"      NextVolume: OutOfWorld";
0271   
0272   if(aStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL)
0273       edm::LogVerbatim("ForwardSim") << "      ProcessName: "<< aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
0274   else
0275       edm::LogVerbatim("ForwardSim") <<"      ProcessName: UserLimit";
0276   
0277 
0278    edm::LogVerbatim("ForwardSim") << std::endl;
0279   */
0280 
0281 #ifdef EDM_ML_DEBUG
0282     if (theTrack->GetNextVolume() != 0)
0283       edm::LogVerbatim("ForwardSim") << " NextVolume: " << theTrack->GetNextVolume()->GetName();
0284     else
0285       edm::LogVerbatim("ForwardSim") << " NextVolume: OutOfWorld";
0286 #endif
0287 
0288     //fill ntuple with step level information
0289     castorstepntuple->Fill(castorsteparray);
0290   }
0291 }
0292 
0293 //================= End of EVENT ===============
0294 void CastorTestAnalysis::update(const EndOfEvent *evt) {
0295   // Look for the Hit Collection
0296   edm::LogVerbatim("ForwardSim") << std::endl
0297                                  << "CastorTest::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
0298 
0299   // access to the G4 hit collections
0300   G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
0301   edm::LogVerbatim("ForwardSim") << "update(*evt) --> accessed all HC";
0302 
0303   int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
0304 
0305   CaloG4HitCollection *theCAFI = (CaloG4HitCollection *)allHC->GetHC(CAFIid);
0306 
0307   theCastorNumScheme = new CastorNumberingScheme();
0308   // CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
0309 
0310   /*
0311   unsigned int volumeID=0;
0312   int det, zside, sector, zmodule;
0313   std::map<int,float,std::less<int> > themap;
0314   double totalEnergy = 0;
0315   double hitEnergy = 0;
0316   double en_in_fi = 0.;
0317   double en_in_pl = 0.;
0318 */
0319   //  double en_in_bu = 0.;
0320   //  double en_in_tu = 0.;
0321 
0322   if (doNTcastorevent) {
0323     eventGlobalHit = 0;
0324     // int eventGlobalHit = 0 ;
0325 
0326     //  Check FI TBranch for Hits
0327     if (theCAFI->entries() > 0)
0328       getCastorBranchData(theCAFI);
0329 
0330     // Find Primary info:
0331     int trackID = 0;
0332 #ifdef EDM_ML_DEBUG
0333     int particleType = 0;
0334 #endif
0335     G4PrimaryParticle *thePrim = nullptr;
0336     G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0337     edm::LogVerbatim("ForwardSim") << "Event has " << nvertex << " vertex";
0338     if (nvertex == 0)
0339       edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event  ERROR: no vertex";
0340 
0341     for (int i = 0; i < nvertex; i++) {
0342       G4PrimaryVertex *avertex = (*evt)()->GetPrimaryVertex(i);
0343       if (avertex == nullptr) {
0344         edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: pointer to vertex = 0";
0345         continue;
0346       }
0347       edm::LogVerbatim("ForwardSim") << "Vertex number :" << i;
0348       int npart = avertex->GetNumberOfParticle();
0349       if (npart == 0)
0350         edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: no primary!";
0351       if (thePrim == nullptr)
0352         thePrim = avertex->GetPrimary(trackID);
0353     }
0354 
0355     double px = 0., py = 0., pz = 0., pInit = 0;
0356 #ifdef EDM_ML_DEBUG
0357     double eta = 0., phi = 0.;
0358 #endif
0359     if (thePrim != nullptr) {
0360       px = thePrim->GetPx();
0361       py = thePrim->GetPy();
0362       pz = thePrim->GetPz();
0363       pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0364       if (pInit == 0) {
0365         edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event  ERR: primary has p=0 ";
0366 #ifdef EDM_ML_DEBUG
0367       } else {
0368         float costheta = pz / pInit;
0369         float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
0370         eta = -log(tan(theta / 2));
0371 
0372         if (px != 0)
0373           phi = atan(py / px);
0374 #endif
0375       }
0376 #ifdef EDM_ML_DEBUG
0377       particleType = thePrim->GetPDGcode();
0378 #endif
0379     } else {
0380       edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: could not find primary ";
0381     }
0382 #ifdef EDM_ML_DEBUG
0383     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Particle Type " << particleType << " p/eta/phi " << pInit
0384                                    << ", " << eta << ", " << phi;
0385 #endif
0386   }
0387 
0388   int iEvt = (*evt)()->GetEventID();
0389   if (iEvt < 10)
0390     edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
0391   else if ((iEvt < 100) && (iEvt % 10 == 0))
0392     edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
0393   else if ((iEvt < 1000) && (iEvt % 100 == 0))
0394     edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
0395   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0396     edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
0397 
0398   edm::LogVerbatim("ForwardSim") << std::endl << "===>>> Done writing user histograms ";
0399 }
0400 
0401 void CastorTestAnalysis::update(const EndOfRun *run) { ; }
0402 
0403 //===================================================================
0404 void CastorTestAnalysis::getCastorBranchData(const CaloG4HitCollection *hc) {
0405   int nentries = hc->entries();
0406 
0407   if (nentries > 0) {
0408     unsigned int volumeID = 0;
0409     int det = 0, zside, sector, zmodule;
0410     std::map<int, float, std::less<int> > themap;
0411     double totalEnergy = 0;
0412     double hitEnergy = 0;
0413     double en_in_sd = 0.;
0414 
0415     for (int ihit = 0; ihit < nentries; ihit++) {
0416       CaloG4Hit *aHit = (*hc)[ihit];
0417       totalEnergy += aHit->getEnergyDeposit();
0418     }
0419 
0420     for (int ihit = 0; ihit < nentries; ihit++) {
0421       CaloG4Hit *aHit = (*hc)[ihit];
0422       volumeID = aHit->getUnitID();
0423       hitEnergy = aHit->getEnergyDeposit();
0424       en_in_sd += aHit->getEnergyDeposit();
0425       //    double enEm = aHit->getEM();
0426       //    double enHad = aHit->getHadr();
0427 
0428       themap[volumeID] += aHit->getEnergyDeposit();
0429       // int det, zside, sector, zmodule;
0430       theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
0431 
0432       // det = 2 ;  //  det=2/3 for CAFI/CAPL
0433 
0434       castoreventarray[ntcastore_evt] = (float)eventIndex;
0435       //    castoreventarray[ntcastore_ihit]      = (float)ihit;
0436       castoreventarray[ntcastore_ihit] = (float)eventGlobalHit;
0437       castoreventarray[ntcastore_detector] = (float)det;
0438       castoreventarray[ntcastore_sector] = (float)sector;
0439       castoreventarray[ntcastore_module] = (float)zmodule;
0440       castoreventarray[ntcastore_enem] = en_in_sd;
0441       castoreventarray[ntcastore_enhad] = totalEnergy;
0442       castoreventarray[ntcastore_hitenergy] = hitEnergy;
0443       castoreventarray[ntcastore_x] = aHit->getPosition().x();
0444       castoreventarray[ntcastore_y] = aHit->getPosition().y();
0445       castoreventarray[ntcastore_z] = aHit->getPosition().z();
0446       //    castoreventarray[ntcastore_x]         = aHit->getEntry().x();
0447       //    castoreventarray[ntcastore_y]         = aHit->getEntry().y();
0448       //    castoreventarray[ntcastore_z]         = aHit->getEntry().z();
0449 
0450       castoreventntuple->Fill(castoreventarray);
0451 
0452       eventGlobalHit++;
0453     }
0454   }  // nentries > 0
0455 }
0456 
0457 //===================================================================
0458 
0459 void CastorTestAnalysis::Finish() {
0460   if (doNTcastorstep) {
0461     castorOutputStepFile->cd();
0462     castorstepntuple->Write();
0463     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple step  written";
0464     castorOutputStepFile->Close();
0465     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Step file closed";
0466   }
0467 
0468   if (doNTcastorevent) {
0469     castorOutputEventFile->cd();
0470     castoreventntuple->Write("", TObject::kOverwrite);
0471     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple event written";
0472     castorOutputEventFile->Close();
0473     edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Event file closed";
0474   }
0475 }
0476 
0477 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0478 #include "FWCore/PluginManager/interface/ModuleDef.h"
0479 
0480 DEFINE_SIMWATCHER(CastorTestAnalysis);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!