#include <SimG4CMS/HcalTestBeam/interface/HcalTB06Analysis.h>

Inheritance diagram for HcalTB06Analysis:

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

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

virtual	~HcalTB06Analysis ()

Public Member Functions inherited from SimProducer
void	registerProducts (edm::ProducerBase &iProd)

	SimProducer ()

Public Member Functions inherited from SimWatcher
	SimWatcher ()

virtual	~SimWatcher ()

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

void	slotForUpdate (const BeginOfRun *iT)

virtual	~Observer ()

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 ()

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

void	slotForUpdate (const G4Step *iT)

virtual	~Observer ()

Private Member Functions
void	clear ()

void	fillBuffer (const EndOfEvent *evt)

void	fillEvent (PHcalTB06Info &)

void	finalAnalysis ()

	HcalTB06Analysis (const HcalTB06Analysis &)

void	init ()

const HcalTB06Analysis &	operator= (const HcalTB06Analysis &)

void	update (const BeginOfRun *run)
	This routine will be called when the appropriate signal arrives. More...

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

void	update (const G4Step *step)
	This routine will be called when the appropriate signal arrives. More...

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

Private Attributes
G4RotationMatrix *	beamline_RM

double	beamOffset

int	count

std::vector< CaloHit >	ecalHitCache

double	eecals

double	ehcals

double	etaInit

double	etots

int	evNum

std::vector< CaloHit >	hcalHitCache

std::vector< CaloHit >	hcalHitLayer

HcalTB06Histo *	histo

int	iceta

int	icphi

std::vector< std::string >	names

int	nPrimary

int	particleType

double	phiInit

double	pInit

bool	pvFound

G4ThreeVector	pvMomentum

G4ThreeVector	pvPosition

int	pvType

G4ThreeVector	pvUVW

std::vector< double >	secEkin

std::vector< G4ThreeVector >	secMomentum

std::vector< int >	secPartID

std::vector< int >	secTrackID

std::vector< int >	shortLivedSecondaries

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

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

Detailed Description

Description: Analysis of 2004 Hcal Test beam simulation

Usage: A Simwatcher class and can be activated from Oscarproducer module

Definition at line 44 of file HcalTB06Analysis.h.

Constructor & Destructor Documentation

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

Definition at line 49 of file HcalTB06Analysis.cc.

References beamline_RM, beamOffset, create_public_lumi_plots::exp, edm::ParameterSet::getParameter(), histo, iceta, icphi, init(), and names.

                                                           : histo(0) {
 
   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB06Analysis");
   names          = m_Anal.getParameter<std::vector<std::string> >("Names");
   beamOffset     =-m_Anal.getParameter<double>("BeamPosition")*cm;
   double fMinEta = m_Anal.getParameter<double>("MinEta");
   double fMaxEta = m_Anal.getParameter<double>("MaxEta");
   double fMinPhi = m_Anal.getParameter<double>("MinPhi");
   double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
   double beamEta = (fMaxEta+fMinEta)/2.;
   double beamPhi = (fMaxPhi+fMinPhi)/2.;
   double beamThet= 2*atan(exp(-beamEta));
   if (beamPhi < 0) beamPhi += twopi;
   iceta          = (int)(beamEta/0.087) + 1;
   icphi          = (int)(fabs(beamPhi)/0.087) + 5;
   if (icphi > 72) icphi -= 73;
 
   produces<PHcalTB06Info>();
 
   beamline_RM = new G4RotationMatrix;
   beamline_RM->rotateZ(-beamPhi);
   beamline_RM->rotateY(-beamThet);
  
   edm::LogInfo("HcalTBSim") << "HcalTB06:: Initialised as observer of BeginOf"
                 << "Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
                 << " with Parameter values:\n \tbeamOffset = " 
                 << beamOffset << "\ticeta = " << iceta 
                 << "\ticphi = " << icphi << "\n\tbeamline_RM = "
                 << *beamline_RM;
 
   init();
 
   histo  = new HcalTB06Histo(m_Anal);
 } 

HcalTB06Analysis::~HcalTB06Analysis ( )

virtual

Definition at line 84 of file HcalTB06Analysis.cc.

References count, and histo.

                                     {
 
   edm::LogInfo("HcalTBSim") << "\n -------->  Total number of selected entries"
                 << " : " << count << "\nPointers:: Histo " <<histo;
   if (histo)   {
     delete histo;
     histo  = 0;
   }
 }

HcalTB06Analysis::HcalTB06Analysis ( const HcalTB06Analysis & )

private

Member Function Documentation

void HcalTB06Analysis::clear ( void )

private

Definition at line 521 of file HcalTB06Analysis.cc.

References ecalHitCache, etaInit, hcalHitCache, nPrimary, particleType, phiInit, pInit, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

Referenced by init(), and update().

                              {
 
   pvFound = false;
   pvType  =-2;
   pvPosition = G4ThreeVector();
   pvMomentum = G4ThreeVector();
   pvUVW      = G4ThreeVector();
   secTrackID.clear();
   secPartID.clear();
   secMomentum.clear();
   secEkin.clear();
   shortLivedSecondaries.clear();
 
   ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end()); 
   hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end()); 
   nPrimary = particleType = 0;
   pInit = etaInit = phiInit = 0;
 }

void HcalTB06Analysis::fillBuffer ( const EndOfEvent * evt )

private

Definition at line 266 of file HcalTB06Analysis.cc.

References TauDecayModes::dec, alignCSCRings::e, ecalHitCache, eta(), etaInit, evNum, CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getEntry(), CaloG4Hit::getTimeSlice(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), GeV, hcalHitCache, i, j, create_public_lumi_plots::log, LogDebug, bookConverter::max, min(), names, npart, nPrimary, AlCaHLTBitMon_ParallelJobs::p, particleType, phi, phiInit, pInit, funct::pow(), python.multivaluedict::sort(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), theta(), and cond::rpcobgas::time.

Referenced by update().

                                                         {
 
   std::vector<CaloHit> hhits;
   int                  idHC, j;
   CaloG4HitCollection* theHC;
   std::map<int,float,std::less<int> > primaries;
   double               etot1=0, etot2=0;
 
   // Look for the Hit Collection of HCal
   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
   std::string sdName = names[0];
   idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
   theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
             << " of ID " << idHC << " is obtained at " << theHC;
 
   if (idHC >= 0 && theHC > 0) {
     hhits.reserve(theHC->entries());
     for (j = 0; j < theHC->entries(); j++) {
       CaloG4Hit* aHit = (*theHC)[j]; 
       double e        = aHit->getEnergyDeposit()/GeV;
       double time     = aHit->getTimeSlice();
       math::XYZPoint pos  = aHit->getEntry();
       unsigned int id = aHit->getUnitID();
       double theta    = pos.theta();
       double eta      = -log(tan(theta * 0.5));
       double phi      = pos.phi();
       CaloHit hit(2,1,e,eta,phi,time,id);
       hhits.push_back(hit);
       primaries[aHit->getTrackID()]+= e;
       etot1 += e;
 #ifdef ddebug
       LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit i/p " << j 
                 << "  ID 0x" << std::hex << id << std::dec 
                 << " time " << std::setw(6) << time << " theta "
                 << std::setw(8) << theta << " eta " << std::setw(8)
                 << eta << " phi " << std::setw(8) << phi << " e " 
                 << std::setw(8) << e;
 #endif
     }
   }
 
   // Add hits in the same channel within same time slice
   std::vector<CaloHit>::iterator itr;
   int nHit = hhits.size();
   std::vector<CaloHit*> hits(nHit);
   for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
     hits[j] = &hhits[j];
   }
   sort(hits.begin(),hits.end(),CaloHitIdMore());
   std::vector<CaloHit*>::iterator k1, k2;
   int nhit = 0;
   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
     int      det    = (**k1).det();
     int      layer  = (**k1).layer();
     double   ehit   = (**k1).e();
     double   eta    = (**k1).eta();
     double   phi    = (**k1).phi();
     double   jitter = (**k1).t();
     uint32_t unitID = (**k1).id();
     int      jump  = 0;
     for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
            unitID==(**k2).id(); k2++) {
       ehit += (**k2).e();
       jump++;
     }
     nhit++;
     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
     hcalHitCache.push_back(hit);
     etot2 += ehit;
     k1    += jump;
 #ifdef ddebug
     LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit store " << nhit 
               << "  ID 0x" << std::hex  << unitID  << std::dec 
               << " time " << std::setw(6) << jitter << " eta "
               << std::setw(8) << eta << " phi " << std::setw(8) 
               << phi  << " e " << std::setw(8) << ehit;
 #endif
   }
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " HCal hits"
             << " from " << nHit << " input hits E(Hcal) " << etot1 
             << " " << etot2;
   
   // Look for the Hit Collection of ECal
   std::vector<CaloHit> ehits;
   sdName= names[1];
   idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
   theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
   etot1 = etot2 = 0;
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
             << " of ID " << idHC << " is obtained at " << theHC;
   if (idHC >= 0 && theHC > 0) {
     ehits.reserve(theHC->entries());
     for (j = 0; j < theHC->entries(); j++) {
       CaloG4Hit* aHit = (*theHC)[j]; 
       double e        = aHit->getEnergyDeposit()/GeV;
       double time     = aHit->getTimeSlice();
       math::XYZPoint pos  = aHit->getEntry();
       unsigned int id = aHit->getUnitID();
       double theta    = pos.theta();
       double eta      = -log(tan(theta * 0.5));
       double phi      = pos.phi();
       if (e < 0 || e > 100000.) e = 0;
       CaloHit hit(1,0,e,eta,phi,time,id);
       ehits.push_back(hit);
       primaries[aHit->getTrackID()]+= e;
       etot1 += e;
 #ifdef ddebug
       LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit i/p " << j 
                 << "  ID 0x" << std::hex << id  << std::dec 
                 << " time " << std::setw(6) << time << " theta " 
                 << std::setw(8) << theta  << " eta " <<std::setw(8)
                 << eta  << " phi " << std::setw(8) << phi << " e "
                 << std::setw(8) << e;
 #endif
     }
   }
 
   // Add hits in the same channel within same time slice
   nHit = ehits.size();
   std::vector<CaloHit*> hite(nHit);
   for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
     hite[j] = &ehits[j];
   }
   sort(hite.begin(),hite.end(),CaloHitIdMore());
   nhit = 0;
   for (k1 = hite.begin(); k1 != hite.end(); k1++) {
     int      det    = (**k1).det();
     int      layer  = (**k1).layer();
     double   ehit   = (**k1).e();
     double   eta    = (**k1).eta();
     double   phi    = (**k1).phi();
     double   jitter = (**k1).t();
     uint32_t unitID = (**k1).id();
     int      jump  = 0;
     for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
            unitID==(**k2).id(); k2++) {
       ehit += (**k2).e();
       jump++;
     }
     nhit++;
     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
     ecalHitCache.push_back(hit);
     etot2 += ehit;
     k1    += jump;
 #ifdef ddebug
     LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit store " << nhit
               << "  ID 0x" << std::hex << unitID  << std::dec 
               << " time " << std::setw(6) << jitter << " eta "
               << std::setw(8) << eta << " phi " << std::setw(8)
               << phi << " e " << std::setw(8) << ehit;
 #endif
   }
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " ECal hits"
             << " from " << nHit << " input hits E(Ecal) " << etot1 
             << " " << etot2;
 
   // Find Primary info:
   nPrimary    = (int)(primaries.size());
   int trackID = 0;
   G4PrimaryParticle* thePrim=0;
   int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Event has " << nvertex 
             << " verteices";
   if (nvertex<=0)
     edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERROR: no "
                   << "vertex found for event " << evNum;
 
   for (int i = 0 ; i<nvertex; i++) {
     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
     if (avertex == 0) {
       edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: pointer "
                 << "to vertex = 0 for event " << evNum;
     } else {
       LogDebug("HcalTBSim") << "HcalTB06Analysis::Vertex number :" << i << " "
                 << avertex->GetPosition();
       int npart = avertex->GetNumberOfParticle();
       if (npart == 0)
     edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::End Of Event ERR: "
                      << "no primary!";
       if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
     }
   }
     
   if (thePrim != 0) {
     double px = thePrim->GetPx();
     double py = thePrim->GetPy();
     double pz = thePrim->GetPz();
     double p  = std::sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
     pInit     = p/GeV;
     if (p==0) 
       edm::LogWarning("HcalTBSim") << "HcalTB06Analysis:: EndOfEvent ERR: "
                    << "primary has p=0 ";
     else {
       double costheta = pz/p;
       double theta = acos(std::min(std::max(costheta,-1.),1.));
       etaInit = -log(tan(theta/2));
       if (px != 0 || py != 0) phiInit = atan2(py,px);  
     }
     particleType = thePrim->GetPDGcode();
   } else 
     edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: could "
                  << "not find primary";
 
 }

void HcalTB06Analysis::fillEvent ( PHcalTB06Info & product )

private

Definition at line 493 of file HcalTB06Analysis.cc.

References ecalHitCache, eecals, ehcals, etaInit, etots, evNum, hcalHitCache, i, nPrimary, particleType, phiInit, pInit, pvMomentum, pvPosition, pvType, pvUVW, PHcalTB06Info::saveHit(), secEkin, secMomentum, secPartID, secTrackID, PHcalTB06Info::setEdep(), PHcalTB06Info::setPrimary(), PHcalTB06Info::setVtxPrim(), PHcalTB06Info::setVtxSec(), x, detailsBasic3DVector::y, and detailsBasic3DVector::z.

Referenced by produce().

                                                         {
 
   //Beam Information
   product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
 
   // Total Energy
   product.setEdep(etots, eecals, ehcals);
 
   //Energy deposits in the crystals and towers
   for (unsigned int i=0; i<hcalHitCache.size(); i++) 
     product.saveHit(hcalHitCache[i].id(), hcalHitCache[i].eta(),
             hcalHitCache[i].phi(), hcalHitCache[i].e(),
             hcalHitCache[i].t());
   for (unsigned int i=0; i<ecalHitCache.size(); i++) 
     product.saveHit(ecalHitCache[i].id(), ecalHitCache[i].eta(),
             ecalHitCache[i].phi(), ecalHitCache[i].e(),
             ecalHitCache[i].t());
 
   //Vertex associated quantities
   product.setVtxPrim(evNum, pvType, pvPosition.x(), pvPosition.y(), 
              pvPosition.z(), pvUVW.x(), pvUVW.y(), pvUVW.z(),
              pvMomentum.x(), pvMomentum.y(), pvMomentum.z());
   for (unsigned int i = 0; i < secTrackID.size(); i++) {
     product.setVtxSec(secTrackID[i], secPartID[i], secMomentum[i].x(),
               secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
   }
 }

void HcalTB06Analysis::finalAnalysis ( )

private

Definition at line 472 of file HcalTB06Analysis.cc.

References ecalHitCache, eecals, ehcals, etaInit, etots, HcalTB06Histo::fillEdep(), HcalTB06Histo::fillPrimary(), hcalHitCache, histo, i, LogDebug, phiInit, and pInit.

Referenced by update().

                                      {
 
   //Beam Information
   histo->fillPrimary(pInit, etaInit, phiInit);
 
   // Total Energy
   eecals = ehcals = 0.;
   for (unsigned int i=0; i<hcalHitCache.size(); i++) {
     ehcals += hcalHitCache[i].e();
   }
   for (unsigned int i=0; i<ecalHitCache.size(); i++) {
     eecals += ecalHitCache[i].e();
   }
   etots = eecals + ehcals;
   LogDebug("HcalTBSim") << "HcalTB06Analysis:: Energy deposit at Sim Level "
             << "(Total) " << etots << " (ECal) " << eecals 
             << " (HCal) " << ehcals;
   histo->fillEdep(etots, eecals, ehcals);
 }

void HcalTB06Analysis::init ( void )

private

Definition at line 105 of file HcalTB06Analysis.cc.

References clear(), count, and evNum.

Referenced by HcalTB06Analysis().

                             {
 
   // counter 
   count = 0;
   evNum = 0;
   clear();
 }

const HcalTB06Analysis& HcalTB06Analysis::operator= ( const HcalTB06Analysis & )

private

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

virtual

Implements SimProducer.

Definition at line 98 of file HcalTB06Analysis.cc.

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

                                                                 {
 
   std::auto_ptr<PHcalTB06Info> product(new PHcalTB06Info);
   fillEvent(*product);
   e.put(product);
 }

void HcalTB06Analysis::update ( const BeginOfRun * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 113 of file HcalTB06Analysis.cc.

                                                     {
 
   int irun = (*run)()->GetRunID();
   edm::LogInfo("HcalTBSim") <<" =====> Begin of Run = " << irun;
  
 }

void HcalTB06Analysis::update ( const BeginOfEvent * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 120 of file HcalTB06Analysis.cc.

References clear(), and evNum.

                                                       {
  
   evNum = (*evt) ()->GetEventID ();
   clear();
   edm::LogInfo("HcalTBSim") << "HcalTB06Analysis: =====> Begin of event = "
                 << evNum;
 }

void HcalTB06Analysis::update ( const G4Step * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 128 of file HcalTB06Analysis.cc.

References LogDebug, NULL, sysUtil::pid, position, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

                                                   {
 
   if (aStep != NULL) {
     //Get Step properties
     G4ThreeVector thePreStepPoint  = aStep->GetPreStepPoint()->GetPosition();
     G4ThreeVector thePostStepPoint;
 
     // Get Tracks properties
     G4Track*      aTrack   = aStep->GetTrack();
     int           trackID  = aTrack->GetTrackID();
     int           parentID = aTrack->GetParentID();
     G4ThreeVector position = aTrack->GetPosition();
     G4ThreeVector momentum = aTrack->GetMomentum();
     G4String      partType = aTrack->GetDefinition()->GetParticleType();
     G4String      partSubType = aTrack->GetDefinition()->GetParticleSubType();
     int    partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
 #ifdef ddebug
     bool   isPDGStable = aTrack->GetDefinition()->GetPDGStable();
 #endif
     double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
     double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
 
     if (!pvFound) { //search for v1
       double stepDeltaEnergy = aStep->GetDeltaEnergy ();
       double kinEnergy = aTrack->GetKineticEnergy ();
       
       // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
       if (trackID == 1 && parentID == 0 && 
       ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
     pvType = -1;
     if (kinEnergy == 0.) {
       pvType = 0;
     } else {
       if (fabs (stepDeltaEnergy / kinEnergy) > 0.1) pvType = 1;
     }
     pvFound    = true;
     pvPosition = position;
     pvMomentum = momentum;
     // Rotated coord.system:
     pvUVW      = (*beamline_RM)*(pvPosition);
 
     //Volume name requires some checks:
     G4String thePostPVname = "NoName";
     G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
     if (thePostPoint) {
       thePostStepPoint = thePostPoint->GetPosition();
       G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
       if (thePostPV) thePostPVname = thePostPV->GetName ();
     }
 #ifdef ddebug
     LogDebug("HcalTBSim") << "HcalTB06Analysis:: V1 found at: " 
                   << thePostStepPoint << " G4VPhysicalVolume: " 
                   << thePostPVname;
 #endif      
     LogDebug("HcalTBSim") << "HcalTB06Analysis::fill_v1Pos: Primary Track "
                   << "momentum: " << pvMomentum << " psoition " 
                   << pvPosition << " u/v/w " << pvUVW;
       }
     } else { 
       // watch for secondaries originating @v1, including the surviving primary
       if ((trackID != 1 && parentID == 1 &&
        (aTrack->GetCurrentStepNumber () == 1) && 
        (thePreStepPoint == pvPosition)) || 
       (trackID == 1 && thePreStepPoint == pvPosition)) {
 #ifdef ddebug
     LogDebug("HcalTBSim") << "HcalTB06Analysis::A secondary...  PDG:" 
                   << partPDGEncoding << " TrackID:" << trackID
                   << " ParentID:" << parentID << " stable: "  
                   << isPDGStable << " Tau: " << pDGlifetime 
                   << " cTauGamma=" 
                   << c_light*pDGlifetime*gammaFactor*1000.
                   << "um" << " GammaFactor: " << gammaFactor;
 #endif      
     secTrackID.push_back(trackID);
     secPartID.push_back(partPDGEncoding);
     secMomentum.push_back(momentum);
     secEkin.push_back(aTrack->GetKineticEnergy());
 
     // Check for short-lived secondaries: cTauGamma<100um
     double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
     if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {//short-lived secondary
       shortLivedSecondaries.push_back(trackID);
       } else {//normal secondary - enter into the V1-calorimetric tree
     //          histos->fill_v1cSec (aTrack);
       }
       }
       // Also watch for tertiary particles coming from 
       // short-lived secondaries from V1
       if (aTrack->GetCurrentStepNumber() == 1) {
     if (shortLivedSecondaries.size() > 0) {
       int pid = parentID;
       std::vector<int>::iterator pos1= shortLivedSecondaries.begin();
       std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
       std::vector<int>::iterator pos;
       for (pos = pos1; pos != pos2; pos++) {
         if (*pos == pid) {//ParentID is on the list of short-lived 
           // secondary 
 #ifdef ddebug
           LogDebug("HcalTBSim") << "HcalTB06Analysis:: A tertiary...  PDG:"
                     << partPDGEncoding << " TrackID:" <<trackID
                     << " ParentID:" << parentID << " stable: "
                     << isPDGStable << " Tau: " << pDGlifetime
                     << " cTauGamma=" 
                     << c_light*pDGlifetime*gammaFactor*1000. 
                     << "um GammaFactor: " << gammaFactor;
 #endif
         }
       }
     }
       }
     }
   }
 }

void HcalTB06Analysis::update ( const EndOfEvent * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 242 of file HcalTB06Analysis.cc.

References count, fillBuffer(), finalAnalysis(), and LogDebug.

                                                     {
 
   count++;
 
   //fill the buffer
   LogDebug("HcalTBSim") << "HcalTB06Analysis::Fill event " 
             << (*evt)()->GetEventID();
   fillBuffer (evt);
   
   //Final Analysis
   LogDebug("HcalTBSim") << "HcalTB06Analysis::Final analysis";  
   finalAnalysis();
 
   int iEvt = (*evt)()->GetEventID();
   if (iEvt < 10) 
     edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
   else if ((iEvt < 100) && (iEvt%10 == 0)) 
     edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
   else if ((iEvt < 1000) && (iEvt%100 == 0)) 
     edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
   else if ((iEvt < 10000) && (iEvt%1000 == 0)) 
     edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
 }