Inheritance diagram for HcalTB02Analysis:

Public Member Functions
	HcalTB02Analysis (const edm::ParameterSet &p)

	HcalTB02Analysis (const HcalTB02Analysis &)=delete

const HcalTB02Analysis &	operator= (const HcalTB02Analysis &)=delete

void	produce (edm::Event &, const edm::EventSetup &) override

	~HcalTB02Analysis () override

Public Member Functions inherited from SimProducer
const SimProducer &	operator= (const SimProducer &)=delete

void	registerProducts (edm::ProducesCollector producesCollector)

	SimProducer ()

	SimProducer (const SimProducer &)=delete

Public Member Functions inherited from SimWatcher
virtual void	beginRun (edm::EventSetup const &)

bool	isMT () const

const SimWatcher &	operator= (const SimWatcher &)=delete

virtual void	registerConsumes (edm::ConsumesCollector)

	SimWatcher ()

	SimWatcher (const SimWatcher &)=delete

virtual	~SimWatcher ()

Public Member Functions inherited from Observer< const BeginOfEvent * >
	Observer ()

void	slotForUpdate (const BeginOfEvent *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfEvent * >
	Observer ()

void	slotForUpdate (const EndOfEvent *iT)

virtual	~Observer ()

Private Member Functions
void	clear ()

void	fillEvent (HcalTB02HistoClass &)

void	finish ()

void	update (const BeginOfEvent *evt) override
	This routine will be called when the appropriate signal arrives. More...

void	update (const EndOfEvent *evt) override
	This routine will be called when the appropriate signal arrives. More...

Private Attributes
float	E5x5Matrix

float	E5x5MatrixN

float	E7x7Matrix

float	E7x7MatrixN

std::map< int, float >	energyInCrystals

std::map< int, float >	energyInScints

double	eta

std::string	fileNameTuple

bool	hcalOnly

std::unique_ptr< HcalTB02Histo >	histo

double	incidentEnergy

int	maxTime

std::vector< std::string >	names

int	particleType

double	phi

double	pInit

std::map< int, float >	primaries

float	SEnergy

float	SEnergyN

float	xE3x3Matrix

float	xE3x3MatrixN

float	xE5x5Matrix

float	xE5x5MatrixN

double	xIncidentEnergy

float	xSEnergy

float	xSEnergyN

Additional Inherited Members
Protected Member Functions inherited from SimProducer
template<class T >
void	produces ()

template<class T >
void	produces (const std::string &instanceName)

Protected Member Functions inherited from SimWatcher
void	setMT (bool val)

Detailed Description

Definition at line 61 of file HcalTB02Analysis.cc.

Constructor & Destructor Documentation

HcalTB02Analysis::HcalTB02Analysis ( const edm::ParameterSet & p )

Definition at line 106 of file HcalTB02Analysis.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hcalOnly, histo, and names.

                                                            {
   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB02Analysis");
   hcalOnly = m_Anal.getUntrackedParameter<bool>("HcalClusterOnly", true);
   names = m_Anal.getParameter<std::vector<std::string> >("Names");
 
   produces<HcalTB02HistoClass>();
 
   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis:: Initialised as observer of "
                                 << "BeginOfJob/BeginOfEvent/EndOfEvent with "
                                 << "Parameter values:\n \thcalOnly = " << hcalOnly;
 
   histo = std::make_unique<HcalTB02Histo>(m_Anal);
 }

HcalTB02Analysis::HcalTB02Analysis ( const HcalTB02Analysis & )

delete

HcalTB02Analysis::~HcalTB02Analysis ( )

override

Definition at line 120 of file HcalTB02Analysis.cc.

References finish().

120 { finish(); }

HcalTB02Analysis::finish

void finish()

Definition: HcalTB02Analysis.cc:443

Member Function Documentation

void HcalTB02Analysis::clear ( void )

private

Definition at line 427 of file HcalTB02Analysis.cc.

References E5x5Matrix, E5x5MatrixN, E7x7Matrix, E7x7MatrixN, energyInCrystals, energyInScints, eta, incidentEnergy, maxTime, particleType, phi, pInit, primaries, SEnergy, SEnergyN, xE3x3Matrix, xE3x3MatrixN, xE5x5Matrix, xE5x5MatrixN, xIncidentEnergy, xSEnergy, and xSEnergyN.

Referenced by update().

                              {
   primaries.clear();
   particleType = 0;
   pInit = eta = phi = incidentEnergy = 0;
 
   SEnergy = E7x7Matrix = E5x5Matrix = SEnergyN = 0;
   E7x7MatrixN = E5x5MatrixN = 0;
   energyInScints.clear();
   maxTime = 0;
 
   xIncidentEnergy = 0;
   energyInCrystals.clear();
   xSEnergy = xSEnergyN = xE5x5Matrix = xE3x3Matrix = 0;
   xE5x5MatrixN = xE3x3MatrixN = 0;
 }

void HcalTB02Analysis::fillEvent ( HcalTB02HistoClass & product )

private

Definition at line 397 of file HcalTB02Analysis.cc.

Referenced by produce().

                                                             {
   //Beam information
   product.set_Nprim(float(primaries.size()));
   product.set_partType(particleType);
   product.set_Einit(pInit / GeV);
   product.set_eta(eta);
   product.set_phi(phi);
   product.set_Eentry(incidentEnergy);
 
   //Calorimeter energy
   product.set_ETot(SEnergy / GeV);
   product.set_E7x7(E7x7Matrix / GeV);
   product.set_E5x5(E5x5Matrix / GeV);
   product.set_ETotN(SEnergyN / GeV);
   product.set_E7x7N(E7x7MatrixN / GeV);
   product.set_E5x5N(E5x5MatrixN / GeV);
   product.set_NUnit(float(energyInScints.size()));
   product.set_Ntimesli(float(maxTime));
 
   //crystal information
   product.set_xEentry(xIncidentEnergy);
   product.set_xNUnit(float(energyInCrystals.size()));
   product.set_xETot(xSEnergy / GeV);
   product.set_xETotN(xSEnergyN / GeV);
   product.set_xE5x5(xE5x5Matrix / GeV);
   product.set_xE3x3(xE3x3Matrix / GeV);
   product.set_xE5x5N(xE5x5MatrixN / GeV);
   product.set_xE3x3N(xE3x3MatrixN / GeV);
 }

void HcalTB02Analysis::finish ( )

private

Definition at line 443 of file HcalTB02Analysis.cc.

Referenced by progressbar.ProgressBar::__next__(), and ~HcalTB02Analysis().

                               {
   /*
   //Profile 
   std::ofstream   oFile;
   oFile.open("profile.dat");
   float st[19] = {0.8,
                   0.4,  0.4,  0.4,  0.4,  0.4,   
                   0.4,  0.4,  0.4,  0.4,  0.4,
                   0.4,  0.4,  0.4,  0.4,  0.4, 
                   0.8,  1.0,  1.0};
                  
   //cm of material (brass) in front of scintillator layer i:
  
   float w[19] = {7.45,                         //iron !
          6.00, 6.00, 6.00, 6.00, 6.00, //brass
          6.00, 6.00, 6.00, 6.00, 6.60, //brass
          6.60, 6.60, 6.60, 6.60, 6.60, //brass
          8.90, 20.65, 19.5};            //brass,iron !
 
   for (int ilayer = 0; ilayer<19; ilayer++) {
 
     // Histogram mean and sigma calculated from the ROOT histos
     edm::LogVerbatim("HcalTBSim") << "Layer number: " << ilayer << " Mean = " 
                   << histo->getMean(ilayer) << " sigma = "   
                   << histo->getRMS(ilayer) << " LThick= "   
                   << w[ilayer] << " SThick= "   << st[ilayer];
       
     oFile << ilayer << " "  << histo->getMean(ilayer) << " " 
       << histo->getRMS(ilayer)  << " " << w[ilayer] << " " << st[ilayer] 
       << std::endl;
 
   } 
   oFile.close(); 
   */
 }

const HcalTB02Analysis& HcalTB02Analysis::operator= ( const HcalTB02Analysis & )

delete

void HcalTB02Analysis::produce	(	edm::Event &	e,
		const edm::EventSetup &
	)

overridevirtual

Implements SimProducer.

Definition at line 126 of file HcalTB02Analysis.cc.

References fillEvent(), eostools::move(), and edm::Event::put().

                                                                 {
   std::unique_ptr<HcalTB02HistoClass> product(new HcalTB02HistoClass);
   fillEvent(*product);
   e.put(std::move(product));
 }

void HcalTB02Analysis::update ( const BeginOfEvent * )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 132 of file HcalTB02Analysis.cc.

References clear().

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                      {
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis: =====> Begin of event = " << (*evt)()->GetEventID();
 #endif
   clear();
 }

void HcalTB02Analysis::update ( const EndOfEvent * )

overrideprivatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 139 of file HcalTB02Analysis.cc.

References E5x5Matrix, E5x5MatrixN, E7x7Matrix, E7x7MatrixN, energyInCrystals, energyInScints, eta, CaloG4Hit::getEM(), CaloG4Hit::getHadr(), CaloG4Hit::getIncidentEnergy(), CaloG4Hit::getTimeSliceID(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), GeV, hcalOnly, histo, mps_fire::i, incidentEnergy, phase1PixelTopology::layer, log, LogDebug, SiStripPI::max, maxTime, min(), names, npart, particleType, phi, pInit, funct::pow(), primaries, sd, SEnergy, SEnergyN, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), theta(), xE3x3Matrix, xE3x3MatrixN, xE5x5Matrix, xE5x5MatrixN, xIncidentEnergy, and xSEnergy.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                    {
   CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
   CLHEP::RandGaussQ randGauss(*engine);
 
   // Look for the Hit Collection
   LogDebug("HcalTBSim") << "HcalTB02Analysis::Fill event " << (*evt)()->GetEventID();
 
   // access to the G4 hit collections
   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
   int ihit = 0;
 
   // HCAL
   std::string sd = names[0];
   int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
   CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
   LogDebug("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << HCHCid << " is obtained at "
                         << theHCHC;
 
   int nentries = 0;
   nentries = theHCHC->entries();
   if (nentries == 0)
     return;
 
   int xentries = 0;
   int XTALSid = 0;
   CaloG4HitCollection* theXTHC = nullptr;
 
   if (!hcalOnly) {
     // XTALS
     sd = names[1];
     XTALSid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
     theXTHC = (CaloG4HitCollection*)allHC->GetHC(XTALSid);
     LogDebug("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << XTALSid
                           << " is obtained at " << theXTHC;
     xentries = theXTHC->entries();
     if (xentries == 0)
       return;
   }
 
   LogDebug("HcalTBSim") << "HcalTB02Analysis :: There are " << nentries << " HCal hits, and" << xentries
                         << " xtal hits";
 
   float ETot = 0., xETot = 0.;
   int scintID = 0, xtalID = 0;
 
   // HCAL
   std::unique_ptr<HcalTB02HcalNumberingScheme> org(new HcalTB02HcalNumberingScheme());
 
   if (HCHCid >= 0 && theHCHC != nullptr) {
     for (ihit = 0; ihit < nentries; ihit++) {
       CaloG4Hit* aHit = (*theHCHC)[ihit];
       scintID = aHit->getUnitID();
       int layer = org->getlayerID(scintID);
       float enEm = aHit->getEM();
       float enhad = aHit->getHadr();
 
       if (layer == 0) {
         enEm = enEm / 2.;
         enhad = enhad / 2.;
       }
 
       energyInScints[scintID] += enEm + enhad;
       primaries[aHit->getTrackID()] += enEm + enhad;
       float time = aHit->getTimeSliceID();
       if (time > maxTime)
         maxTime = (int)time;
       histo->fillAllTime(time);
     }
 
     //
     // Profile
     //
 
     float TowerEne[8][18], TowerEneCF[8][18], LayerEne[19], EnRing[100];
     for (int iphi = 0; iphi < 8; iphi++) {
       for (int jeta = 0; jeta < 18; jeta++) {
         TowerEne[iphi][jeta] = 0.;
         TowerEneCF[iphi][jeta] = 0.;
       }
     }
 
     for (int ilayer = 0; ilayer < 19; ilayer++)
       LayerEne[ilayer] = 0.;
     for (int iring = 0; iring < 100; iring++)
       EnRing[iring] = 0.;
 
     for (std::map<int, float>::iterator is = energyInScints.begin(); is != energyInScints.end(); is++) {
       ETot = (*is).second;
 
       int layer = org->getlayerID((*is).first);
 
       if ((layer != 17) && (layer != 18)) {
         float eta = org->getetaID((*is).first);
         float phi = org->getphiID((*is).first);
 
         SEnergy += ETot;
         TowerEne[(int)phi][(int)eta] += ETot;
 
         TowerEneCF[(int)phi][(int)eta] += ETot * (1. + 0.1 * randGauss.fire());
         double dR = 0.08727 * std::sqrt((eta - 8.) * (eta - 8.) + (phi - 3.) * (phi - 3.));
         EnRing[(int)(dR / 0.01)] += ETot;
       }
 
       LayerEne[layer] += ETot;
     }
 
     for (int ilayer = 0; ilayer < 19; ilayer++) {
       histo->fillProfile(ilayer, LayerEne[ilayer] / GeV);
     }
 
     for (int iring = 0; iring < 100; iring++)
       histo->fillTransProf(0.01 * iring + 0.005, EnRing[iring] / GeV);
 
     for (int iphi = 0; iphi < 8; iphi++) {
       for (int jeta = 0; jeta < 18; jeta++) {
         SEnergyN += TowerEneCF[iphi][jeta] + 3. * randGauss.fire();  // QGSP
 
         double Rand = 3. * randGauss.fire();  // QGSP
 
         if ((iphi >= 0) && (iphi < 7)) {
           if ((jeta >= 5) && (jeta < 12)) {
             E7x7Matrix += TowerEne[iphi][jeta];
             E7x7MatrixN += TowerEneCF[iphi][jeta] + Rand;
 
             if ((iphi >= 1) && (iphi < 6)) {
               if ((jeta >= 6) && (jeta < 11)) {
                 E5x5Matrix += TowerEne[iphi][jeta];
                 E5x5MatrixN += TowerEneCF[iphi][jeta] + Rand;
               }
             }
           }
         }
       }
     }
 
     //
     // Find Primary info:
     //
     int trackID = 0;
     G4PrimaryParticle* thePrim = nullptr;
     G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
     LogDebug("HcalTBSim") << "HcalTB02Analysis :: Event has " << nvertex << " vertex";
     if (nvertex == 0)
       edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event  "
                                    << "ERROR: no vertex";
 
     for (int i = 0; i < nvertex; i++) {
       G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
       if (avertex == nullptr) {
         edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
                                      << "ERROR: pointer to vertex = 0";
       } else {
         int npart = avertex->GetNumberOfParticle();
         LogDebug("HcalTBSim") << "HcalTB02Analysis::Vertex number :" << i << " with " << npart << " particles";
         if (thePrim == nullptr)
           thePrim = avertex->GetPrimary(trackID);
       }
     }
 
     double px = 0., py = 0., pz = 0.;
 
     if (thePrim != nullptr) {
       px = thePrim->GetPx();
       py = thePrim->GetPy();
       pz = thePrim->GetPz();
       pInit = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
       if (pInit == 0) {
         edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
                                      << " ERROR: primary has p=0 ";
       } else {
         float costheta = pz / pInit;
         float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
         eta = -log(tan(theta / 2));
         if (px != 0)
           phi = atan(py / px);
       }
       particleType = thePrim->GetPDGcode();
     } else {
       LogDebug("HcalTBSim") << "HcalTB02Analysis:: End Of Event ERROR: could"
                             << " not find primary ";
     }
 
     CaloG4Hit* firstHit = (*theHCHC)[0];
     incidentEnergy = (firstHit->getIncidentEnergy() / GeV);
 
   }  // number of Hits > 0
 
   if (!hcalOnly) {
     // XTALS
 
     if (XTALSid >= 0 && theXTHC != nullptr) {
       for (int xihit = 0; xihit < xentries; xihit++) {
         CaloG4Hit* xaHit = (*theXTHC)[xihit];
 
         float xenEm = xaHit->getEM();
         float xenhad = xaHit->getHadr();
         xtalID = xaHit->getUnitID();
 
         energyInCrystals[xtalID] += xenEm + xenhad;
       }
 
       float xCrysEne[7][7];
       for (int irow = 0; irow < 7; irow++) {
         for (int jcol = 0; jcol < 7; jcol++) {
           xCrysEne[irow][jcol] = 0.;
         }
       }
 
       for (std::map<int, float>::iterator is = energyInCrystals.begin(); is != energyInCrystals.end(); is++) {
         int xtalID = (*is).first;
         xETot = (*is).second;
 
         int irow = (int)(xtalID / 100.);
         int jcol = (int)(xtalID - 100. * irow);
 
         xSEnergy += xETot;
         xCrysEne[irow][jcol] = xETot;
 
         float dR = std::sqrt(0.01619 * 0.01619 * (jcol - 3) * (jcol - 3) + 0.01606 * 0.01606 * (irow - 3) * (irow - 3));
         histo->fillTransProf(dR, xETot * 1.05);
 
         if ((irow > 0) && (irow < 6)) {
           if ((jcol > 0) && (jcol < 6)) {
             xE5x5Matrix += xCrysEne[irow][jcol];
             xE5x5MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
 
             if ((irow > 1) && (irow < 5)) {
               if ((jcol > 1) && (jcol < 5)) {
                 xE3x3Matrix += xCrysEne[irow][jcol];
                 xE3x3MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
               }
             }
           }
         }
       }
 
       if (!hcalOnly) {
         //  assert(theXTHC);
         if (theXTHC != nullptr) {
           CaloG4Hit* xfirstHit = (*theXTHC)[0];
           xIncidentEnergy = xfirstHit->getIncidentEnergy() / GeV;
         }
       }
 
     }  // number of Hits > 0
   }
 
   int iEvt = (*evt)()->GetEventID();
   if (iEvt < 10)
     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
   else if ((iEvt < 100) && (iEvt % 10 == 0))
     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
   else if ((iEvt < 1000) && (iEvt % 100 == 0))
     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
 }