39 #include "G4SDManager.hh"
40 #include "G4VProcess.hh"
41 #include "G4HCofThisEvent.hh"
42 #include "CLHEP/Units/GlobalSystemOfUnits.h"
43 #include "CLHEP/Units/GlobalPhysicalConstants.h"
58 double beamEta = (fMaxEta+fMinEta)/2.;
59 double beamPhi = (fMaxPhi+fMinPhi)/2.;
60 double beamThet= 2*atan(
exp(-beamEta));
61 if (beamPhi < 0) beamPhi += twopi;
62 iceta = (int)(beamEta/0.087) + 1;
63 icphi = (int)(fabs(beamPhi)/0.087) + 5;
66 produces<PHcalTB06Info>();
69 beamline_RM->rotateZ(-beamPhi);
70 beamline_RM->rotateY(-beamThet);
72 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06:: Initialised as observer of BeginOf"
73 <<
"Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
74 <<
" with Parameter values:\n \tbeamOffset = "
75 << beamOffset <<
"\ticeta = " <<
iceta
76 <<
"\ticphi = " <<
icphi <<
"\n\tbeamline_RM = "
86 edm::LogInfo(
"HcalTBSim") <<
"\n --------> Total number of selected entries"
87 <<
" : " <<
count <<
"\nPointers:: Histo " <<
histo;
115 int irun = (*run)()->GetRunID();
116 edm::LogInfo(
"HcalTBSim") <<
" =====> Begin of Run = " << irun;
122 evNum = (*evt) ()->GetEventID ();
124 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis: =====> Begin of event = "
132 G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
133 G4ThreeVector thePostStepPoint;
136 G4Track* aTrack = aStep->GetTrack();
137 int trackID = aTrack->GetTrackID();
138 int parentID = aTrack->GetParentID();
139 G4ThreeVector
position = aTrack->GetPosition();
140 G4ThreeVector momentum = aTrack->GetMomentum();
141 G4String partType = aTrack->GetDefinition()->GetParticleType();
142 G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
143 int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
145 bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
147 double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
148 double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
151 double stepDeltaEnergy = aStep->GetDeltaEnergy ();
152 double kinEnergy = aTrack->GetKineticEnergy ();
155 if (trackID == 1 && parentID == 0 &&
156 ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
158 if (kinEnergy == 0.) {
161 if (fabs (stepDeltaEnergy / kinEnergy) > 0.1)
pvType = 1;
170 G4String thePostPVname =
"NoName";
171 G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
173 thePostStepPoint = thePostPoint->GetPosition();
174 G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
175 if (thePostPV) thePostPVname = thePostPV->GetName ();
178 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: V1 found at: "
179 << thePostStepPoint <<
" G4VPhysicalVolume: "
182 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::fill_v1Pos: Primary Track "
188 if ((trackID != 1 && parentID == 1 &&
189 (aTrack->GetCurrentStepNumber () == 1) &&
191 (trackID == 1 && thePreStepPoint ==
pvPosition)) {
193 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::A secondary... PDG:"
194 << partPDGEncoding <<
" TrackID:" << trackID
195 <<
" ParentID:" << parentID <<
" stable: "
196 << isPDGStable <<
" Tau: " << pDGlifetime
198 << c_light*pDGlifetime*gammaFactor*1000.
199 <<
"um" <<
" GammaFactor: " << gammaFactor;
204 secEkin.push_back(aTrack->GetKineticEnergy());
207 double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
208 if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {
216 if (aTrack->GetCurrentStepNumber() == 1) {
221 std::vector<int>::iterator pos;
222 for (pos = pos1; pos != pos2; pos++) {
226 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: A tertiary... PDG:"
227 << partPDGEncoding <<
" TrackID:" <<trackID
228 <<
" ParentID:" << parentID <<
" stable: "
229 << isPDGStable <<
" Tau: " << pDGlifetime
231 << c_light*pDGlifetime*gammaFactor*1000.
232 <<
"um GammaFactor: " << gammaFactor;
247 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Fill event "
248 << (*evt)()->GetEventID();
252 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Final analysis";
255 int iEvt = (*evt)()->GetEventID();
257 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
258 else if ((iEvt < 100) && (iEvt%10 == 0))
259 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
260 else if ((iEvt < 1000) && (iEvt%100 == 0))
261 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
262 else if ((iEvt < 10000) && (iEvt%1000 == 0))
263 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
268 std::vector<CaloHit> hhits;
271 std::map<int,float,std::less<int> > primaries;
272 double etot1=0, etot2=0;
275 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
277 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
279 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hit Collection for " << sdName
280 <<
" of ID " << idHC <<
" is obtained at " << theHC;
282 if (idHC >= 0 && theHC > 0) {
283 hhits.reserve(theHC->entries());
284 for (j = 0; j < theHC->entries(); j++) {
290 double theta = pos.theta();
292 double phi = pos.phi();
294 hhits.push_back(hit);
298 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hcal Hit i/p " << j
299 <<
" ID 0x" << std::hex <<
id << std::dec
300 <<
" time " << std::setw(6) << time <<
" theta "
301 << std::setw(8) << theta <<
" eta " << std::setw(8)
302 << eta <<
" phi " << std::setw(8) << phi <<
" e "
303 << std::setw(8) <<
e;
309 std::vector<CaloHit>::iterator itr;
310 int nHit = hhits.size();
311 std::vector<CaloHit*> hits(nHit);
312 for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
316 std::vector<CaloHit*>::iterator k1, k2;
318 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
319 int det = (**k1).det();
320 int layer = (**k1).layer();
321 double ehit = (**k1).e();
322 double eta = (**k1).eta();
323 double phi = (**k1).phi();
324 double jitter = (**k1).t();
325 uint32_t unitID = (**k1).id();
327 for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
328 unitID==(**k2).id(); k2++) {
333 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
338 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hcal Hit store " << nhit
339 <<
" ID 0x" << std::hex << unitID << std::dec
340 <<
" time " << std::setw(6) << jitter <<
" eta "
341 << std::setw(8) << eta <<
" phi " << std::setw(8)
342 << phi <<
" e " << std::setw(8) << ehit;
345 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Stores " << nhit <<
" HCal hits"
346 <<
" from " << nHit <<
" input hits E(Hcal) " << etot1
350 std::vector<CaloHit> ehits;
352 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
355 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hit Collection for " << sdName
356 <<
" of ID " << idHC <<
" is obtained at " << theHC;
357 if (idHC >= 0 && theHC > 0) {
358 ehits.reserve(theHC->entries());
359 for (j = 0; j < theHC->entries(); j++) {
365 double theta = pos.theta();
367 double phi = pos.phi();
368 if (e < 0 || e > 100000.) e = 0;
370 ehits.push_back(hit);
374 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Ecal Hit i/p " << j
375 <<
" ID 0x" << std::hex <<
id << std::dec
376 <<
" time " << std::setw(6) << time <<
" theta "
377 << std::setw(8) << theta <<
" eta " <<std::setw(8)
378 << eta <<
" phi " << std::setw(8) << phi <<
" e "
379 << std::setw(8) <<
e;
386 std::vector<CaloHit*> hite(nHit);
387 for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
392 for (k1 = hite.begin(); k1 != hite.end(); k1++) {
393 int det = (**k1).det();
394 int layer = (**k1).layer();
395 double ehit = (**k1).e();
396 double eta = (**k1).eta();
397 double phi = (**k1).phi();
398 double jitter = (**k1).t();
399 uint32_t unitID = (**k1).id();
401 for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
402 unitID==(**k2).id(); k2++) {
407 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
412 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Ecal Hit store " << nhit
413 <<
" ID 0x" << std::hex << unitID << std::dec
414 <<
" time " << std::setw(6) << jitter <<
" eta "
415 << std::setw(8) << eta <<
" phi " << std::setw(8)
416 << phi <<
" e " << std::setw(8) << ehit;
419 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Stores " << nhit <<
" ECal hits"
420 <<
" from " << nHit <<
" input hits E(Ecal) " << etot1
426 G4PrimaryParticle* thePrim=0;
427 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
428 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Event has " << nvertex
431 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERROR: no "
432 <<
"vertex found for event " <<
evNum;
434 for (
int i = 0 ;
i<nvertex;
i++) {
435 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(
i);
437 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERR: pointer "
438 <<
"to vertex = 0 for event " <<
evNum;
440 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Vertex number :" <<
i <<
" "
441 << avertex->GetPosition();
442 int npart = avertex->GetNumberOfParticle();
446 if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
451 double px = thePrim->GetPx();
452 double py = thePrim->GetPy();
453 double pz = thePrim->GetPz();
458 <<
"primary has p=0 ";
460 double costheta = pz/
p;
463 if (px != 0 || py != 0)
phiInit = atan2(py,px);
467 edm::LogWarning(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERR: could "
468 <<
"not find primary";
486 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Energy deposit at Sim Level "
527 pvUVW = G4ThreeVector();
T getParameter(std::string const &) const
#define DEFINE_SIMWATCHER(type)
void update(const BeginOfRun *run)
This routine will be called when the appropriate signal arrives.
std::vector< std::string > names
std::vector< CaloHit > ecalHitCache
void setEdep(double simtot, double sime, double simh)
void fillEdep(double etots, double eecals, double ehcals)
Geom::Theta< T > theta() const
void setVtxPrim(int evNum, int type, double x, double y, double z, double u, double v, double w, double px, double py, double pz)
virtual void produce(edm::Event &, const edm::EventSetup &)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
void fillEvent(PHcalTB06Info &)
Tan< T >::type tan(const T &t)
void setVtxSec(int id, int pdg, double px, double py, double pz, double ek)
void saveHit(unsigned int det, double eta, double phi, double e, double t)
std::vector< CaloHit > hcalHitCache
std::vector< G4ThreeVector > secMomentum
std::vector< int > shortLivedSecondaries
XYZPointD XYZPoint
point in space with cartesian internal representation
void setPrimary(int primary, int id, double energy, double eta, double phi)
void fillBuffer(const EndOfEvent *evt)
HcalTB06Analysis(const edm::ParameterSet &p)
std::vector< double > secEkin
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
static int position[264][3]
double getTimeSlice() const
math::XYZPoint getEntry() const
std::vector< int > secPartID
uint32_t getUnitID() const
void fillPrimary(double energy, double eta, double phi)
virtual ~HcalTB06Analysis()
std::vector< int > secTrackID
G4RotationMatrix * beamline_RM
Power< A, B >::type pow(const A &a, const B &b)
double getEnergyDeposit() const