40 #include "G4SDManager.hh"
41 #include "G4VProcess.hh"
42 #include "G4HCofThisEvent.hh"
43 #include "CLHEP/Units/GlobalSystemOfUnits.h"
44 #include "CLHEP/Units/GlobalPhysicalConstants.h"
59 double beamEta = (fMaxEta+fMinEta)/2.;
60 double beamPhi = (fMaxPhi+fMinPhi)/2.;
61 double beamThet= 2*atan(
exp(-beamEta));
62 if (beamPhi < 0) beamPhi += twopi;
63 iceta = (int)(beamEta/0.087) + 1;
64 icphi = (int)(fabs(beamPhi)/0.087) + 5;
67 produces<PHcalTB06Info>();
70 beamline_RM->rotateZ(-beamPhi);
71 beamline_RM->rotateY(-beamThet);
73 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06:: Initialised as observer of BeginOf"
74 <<
"Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
75 <<
" with Parameter values:\n \tbeamOffset = "
76 << beamOffset <<
"\ticeta = " <<
iceta
77 <<
"\ticphi = " <<
icphi <<
"\n\tbeamline_RM = "
87 edm::LogInfo(
"HcalTBSim") <<
"\n --------> Total number of selected entries"
88 <<
" : " <<
count <<
"\nPointers:: Histo " <<
histo;
116 int irun = (*run)()->GetRunID();
117 edm::LogInfo(
"HcalTBSim") <<
" =====> Begin of Run = " << irun;
123 evNum = (*evt) ()->GetEventID ();
125 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis: =====> Begin of event = "
133 G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
134 G4ThreeVector thePostStepPoint;
137 G4Track* aTrack = aStep->GetTrack();
138 int trackID = aTrack->GetTrackID();
139 int parentID = aTrack->GetParentID();
140 G4ThreeVector
position = aTrack->GetPosition();
141 G4ThreeVector momentum = aTrack->GetMomentum();
142 G4String partType = aTrack->GetDefinition()->GetParticleType();
143 G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
144 int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
146 bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
148 double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
149 double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
152 double stepDeltaEnergy = aStep->GetDeltaEnergy ();
153 double kinEnergy = aTrack->GetKineticEnergy ();
156 if (trackID == 1 && parentID == 0 &&
157 ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
159 if (kinEnergy == 0.) {
162 if (fabs (stepDeltaEnergy / kinEnergy) > 0.1)
pvType = 1;
171 G4String thePostPVname =
"NoName";
172 G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
174 thePostStepPoint = thePostPoint->GetPosition();
175 G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
176 if (thePostPV) thePostPVname = thePostPV->GetName ();
179 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: V1 found at: "
180 << thePostStepPoint <<
" G4VPhysicalVolume: "
183 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::fill_v1Pos: Primary Track "
189 if ((trackID != 1 && parentID == 1 &&
190 (aTrack->GetCurrentStepNumber () == 1) &&
192 (trackID == 1 && thePreStepPoint ==
pvPosition)) {
194 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::A secondary... PDG:"
195 << partPDGEncoding <<
" TrackID:" << trackID
196 <<
" ParentID:" << parentID <<
" stable: "
197 << isPDGStable <<
" Tau: " << pDGlifetime
199 << c_light*pDGlifetime*gammaFactor*1000.
200 <<
"um" <<
" GammaFactor: " << gammaFactor;
205 secEkin.push_back(aTrack->GetKineticEnergy());
208 double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
209 if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {
217 if (aTrack->GetCurrentStepNumber() == 1) {
222 std::vector<int>::iterator
pos;
223 for (pos = pos1; pos != pos2; pos++) {
227 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: A tertiary... PDG:"
228 << partPDGEncoding <<
" TrackID:" <<trackID
229 <<
" ParentID:" << parentID <<
" stable: "
230 << isPDGStable <<
" Tau: " << pDGlifetime
232 << c_light*pDGlifetime*gammaFactor*1000.
233 <<
"um GammaFactor: " << gammaFactor;
248 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Fill event "
249 << (*evt)()->GetEventID();
253 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Final analysis";
256 int iEvt = (*evt)()->GetEventID();
258 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
259 else if ((iEvt < 100) && (iEvt%10 == 0))
260 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
261 else if ((iEvt < 1000) && (iEvt%100 == 0))
262 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
263 else if ((iEvt < 10000) && (iEvt%1000 == 0))
264 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis:: Event " << iEvt;
269 std::vector<CaloHit> hhits;
272 std::map<int,float,std::less<int> > primaries;
273 double etot1=0, etot2=0;
276 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
278 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
280 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hit Collection for " << sdName
281 <<
" of ID " << idHC <<
" is obtained at " << theHC;
283 if (idHC >= 0 && theHC > 0) {
284 hhits.reserve(theHC->entries());
285 for (j = 0; j < theHC->entries(); j++) {
291 double theta = pos.theta();
293 double phi = pos.phi();
295 hhits.push_back(hit);
299 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hcal Hit i/p " << j
300 <<
" ID 0x" << std::hex <<
id << std::dec
301 <<
" time " << std::setw(6) << time <<
" theta "
302 << std::setw(8) << theta <<
" eta " << std::setw(8)
303 << eta <<
" phi " << std::setw(8) << phi <<
" e "
304 << std::setw(8) <<
e;
310 std::vector<CaloHit>::iterator itr;
311 int nHit = hhits.size();
312 std::vector<CaloHit*> hits(nHit);
313 for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
317 std::vector<CaloHit*>::iterator k1, k2;
319 for (k1 = hits.begin(); k1 != hits.end(); k1++) {
320 int det = (**k1).det();
321 int layer = (**k1).layer();
322 double ehit = (**k1).e();
323 double eta = (**k1).eta();
324 double phi = (**k1).phi();
325 double jitter = (**k1).t();
326 uint32_t unitID = (**k1).id();
328 for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
329 unitID==(**k2).id(); k2++) {
334 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
339 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hcal Hit store " << nhit
340 <<
" ID 0x" << std::hex << unitID << std::dec
341 <<
" time " << std::setw(6) << jitter <<
" eta "
342 << std::setw(8) << eta <<
" phi " << std::setw(8)
343 << phi <<
" e " << std::setw(8) << ehit;
346 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Stores " << nhit <<
" HCal hits"
347 <<
" from " << nHit <<
" input hits E(Hcal) " << etot1
351 std::vector<CaloHit> ehits;
353 idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
356 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Hit Collection for " << sdName
357 <<
" of ID " << idHC <<
" is obtained at " << theHC;
358 if (idHC >= 0 && theHC > 0) {
359 ehits.reserve(theHC->entries());
360 for (j = 0; j < theHC->entries(); j++) {
366 double theta = pos.theta();
368 double phi = pos.phi();
369 if (e < 0 || e > 100000.) e = 0;
371 ehits.push_back(hit);
375 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Ecal Hit i/p " << j
376 <<
" ID 0x" << std::hex <<
id << std::dec
377 <<
" time " << std::setw(6) << time <<
" theta "
378 << std::setw(8) << theta <<
" eta " <<std::setw(8)
379 << eta <<
" phi " << std::setw(8) << phi <<
" e "
380 << std::setw(8) <<
e;
387 std::vector<CaloHit*> hite(nHit);
388 for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
393 for (k1 = hite.begin(); k1 != hite.end(); k1++) {
394 int det = (**k1).det();
395 int layer = (**k1).layer();
396 double ehit = (**k1).e();
397 double eta = (**k1).eta();
398 double phi = (**k1).phi();
399 double jitter = (**k1).t();
400 uint32_t unitID = (**k1).id();
402 for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
403 unitID==(**k2).id(); k2++) {
408 CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
413 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Ecal Hit store " << nhit
414 <<
" ID 0x" << std::hex << unitID << std::dec
415 <<
" time " << std::setw(6) << jitter <<
" eta "
416 << std::setw(8) << eta <<
" phi " << std::setw(8)
417 << phi <<
" e " << std::setw(8) << ehit;
420 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Stores " << nhit <<
" ECal hits"
421 <<
" from " << nHit <<
" input hits E(Ecal) " << etot1
427 G4PrimaryParticle* thePrim=0;
428 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
429 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Event has " << nvertex
432 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERROR: no "
433 <<
"vertex found for event " <<
evNum;
435 for (
int i = 0 ;
i<nvertex;
i++) {
436 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(
i);
438 edm::LogInfo(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERR: pointer "
439 <<
"to vertex = 0 for event " <<
evNum;
441 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis::Vertex number :" <<
i <<
" "
442 << avertex->GetPosition();
443 int npart = avertex->GetNumberOfParticle();
447 if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
452 double px = thePrim->GetPx();
453 double py = thePrim->GetPy();
454 double pz = thePrim->GetPz();
459 <<
"primary has p=0 ";
461 double costheta = pz/
p;
464 if (px != 0 || py != 0)
phiInit = atan2(py,px);
468 edm::LogWarning(
"HcalTBSim") <<
"HcalTB06Analysis::EndOfEvent ERR: could "
469 <<
"not find primary";
487 LogDebug(
"HcalTBSim") <<
"HcalTB06Analysis:: Energy deposit at Sim Level "
528 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)
static int position[TOTALCHAMBERS][3]
virtual void produce(edm::Event &, const edm::EventSetup &)
const T & max(const T &a, const T &b)
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
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