15 #include "G4VPhysicalVolume.hh"
18 #include "G4NavigationHistory.hh"
19 #include "Randomize.hh"
21 #include "CLHEP/Units/GlobalPhysicalConstants.h"
22 #include "CLHEP/Units/GlobalSystemOfUnits.h"
35 bool useShowerLibrary = m_HF.
getParameter<
bool>(
"UseShowerLibrary");
40 double lambdaMean = m_HF.
getParameter<
double>(
"LambdaMean");
43 edm::LogInfo(
"HFShower") <<
"HFShowerParam::Use of shower library is set to "
44 << useShowerLibrary <<
" Use of Gflash is set to "
45 << useGflash <<
" P.E. per GeV " << pePerGeV
46 <<
", ref. index of fibre " << ref_index
47 <<
", Track EM Flag " << trackEM <<
", edMin "
48 << edMin <<
" GeV, use of Short fibre info in"
49 <<
" shower library set to " << !(
onlyLong)
50 <<
", use of parametrization for last part set to "
51 << parametrizeLast <<
", Mean lambda " << lambdaMean
52 <<
", Application of Fiducial Cut " << applyFidCut;
58 LogDebug(
"HFShower") <<
"HFShowerParam::Save histos in directory "
59 <<
"ProfileFromParam";
61 hzv = showerDir.
make<TH1F>(
"hzv",
"Longitudinal Profile;Radiation Length;Number of Spots",330,0.0,1650.0);
62 em_2d_1 = showerDir.
make<TH2F>(
"em_2d_1",
"Lateral Profile vs. Shower Depth;cm;Events",800,800.0,1600.0,100,50.0,150.0);
63 em_long_1 = showerDir.
make<TH1F>(
"em_long_1",
"Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
64 em_long_1_tuned = showerDir.
make<TH1F>(
"em_long_1_tuned",
"Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
65 em_lateral_1 = showerDir.
make<TH1F>(
"em_lateral_1",
"Lateral Profile;cm;Events",100,50.0,150.0);
66 em_2d_2 = showerDir.
make<TH2F>(
"em_2d_2",
"Lateral Profile vs. Shower Depth;cm;Events",800,800.0,1600.0,100,50.0,150.0);
67 em_long_2 = showerDir.
make<TH1F>(
"em_long_2",
"Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
68 em_lateral_2 = showerDir.
make<TH1F>(
"em_lateral_2",
"Lateral Profile;cm;Events",100,50.0,150.0);
69 em_long_gflash = showerDir.
make<TH1F>(
"em_long_gflash",
"Longitudinal Profile From GFlash;cm;Number of Spots",800,800.0,1600.0);
70 em_long_sl = showerDir.
make<TH1F>(
"em_long_sl",
"Longitudinal Profile From Shower Library;cm;Number of Spots",800,800.0,1600.0);
73 LogDebug(
"HFShower") <<
"HFShowerParam::No file is available for "
74 <<
"saving histos so the flag is set to false";
78 G4String attribute =
"ReadOutName";
90 edm::LogInfo(
"HFShower") <<
"HFShowerParam: " <<gpar.size() <<
" gpar (cm)";
91 for (
unsigned int ig=0; ig<gpar.size(); ig++)
92 edm::LogInfo(
"HFShower") <<
"HFShowerParam: gpar[" << ig <<
"] = "
93 << gpar[ig]/cm <<
" cm";
95 edm::LogError(
"HFShower") <<
"HFShowerParam: cannot get filtered "
96 <<
" view for " << attribute <<
" matching " <<
name;
97 throw cms::Exception(
"Unknown",
"HFShowerParam") <<
"cannot match " << attribute
98 <<
" to " << name <<
"\n";
105 edm::LogInfo(
"HFShower") <<
"att. length used for (lambda=" << lambdaMean
116 emPDG = theParticleTable->FindParticle(
"e-")->GetPDGEncoding();
117 epPDG = theParticleTable->FindParticle(
"e+")->GetPDGEncoding();
118 gammaPDG = theParticleTable->FindParticle(
"gamma")->GetPDGEncoding();
127 G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
128 G4Track * track = aStep->GetTrack();
129 G4ThreeVector hitPoint = preStepPoint->GetPosition();
130 G4int particleCode = track->GetDefinition()->GetPDGEncoding();
132 G4ThreeVector localPoint = G4ThreeVector(hitPoint.x(),hitPoint.y(),zv);
134 double pin = (preStepPoint->GetTotalEnergy())/GeV;
135 double zint = hitPoint.z();
140 << track->GetDefinition()->GetParticleName()
141 <<
" of energy " << pin <<
" GeV"
142 <<
" Pos x,y,z = " << hitPoint.x() <<
","
143 << hitPoint.y() <<
"," << zint <<
" (" << zz <<
","
144 << localPoint.z() <<
", "
145 << (localPoint.z()+0.5*
gpar[1]) <<
") Local "
148 std::vector<HFShowerParam::Hit> hits;
154 double pBeta = track->GetDynamicParticle()->GetTotalMomentum() / track->GetDynamicParticle()->GetTotalEnergy();
156 if (track->GetDefinition()->GetPDGCharge() != 0 && pBeta > (1/
ref_index) &&
157 aStep->GetTotalEnergyDeposit() > 0) other =
true;
161 if (particleCode ==
emPDG || particleCode ==
epPDG ||
162 particleCode ==
gammaPDG || other) {
169 }
else if (track->GetDefinition()->GetPDGCharge() != 0 &&
171 edep = (aStep->GetTotalEnergyDeposit())/GeV;
173 std::string
path =
"ShowerLibrary";
175 edm::LogInfo(
"HFShower") <<
"HFShowerParam: getHits edep = " << edep
176 <<
", Kill = " << kill <<
", pin = " << pin
177 <<
", edMin = " <<
edMin <<
" Other " << other;
183 for (
unsigned int i=0;
i<hitSL.size();
i++) {
190 if (npmt <= 0) ok =
false;
194 hit.
depth = hitSL[
i].depth;
195 hit.
time = hitSL[
i].time;
204 if (hit.
depth == 1) {
208 }
else if (hit.
depth == 2) {
214 edm::LogInfo(
"HFShower") <<
"HFShowerParam: Hit at depth "
215 << hit.
depth <<
" with edep " << hit.
edep
216 <<
" Time " << hit.
time;
222 for (
unsigned int i=0;
i<hitSL.size(); ++
i) {
224 G4ThreeVector pe_effect(hitSL[
i].
position.x(), hitSL[
i].position.y(),
225 hitSL[
i].position.z());
230 if (zv < 0. || zv >
gpar[1]) {
232 std::cout<<
"-#Zcut-HFShowerParam::getHits:z="<<zv<<
",m="<<gpar[1]<<std::endl;
239 edm::LogInfo(
"HFShower") <<
"HFShowerParam::getHits:#PMT= "
240 << npmt <<
",z = " << zv;
244 edm::LogInfo(
"HFShower") <<
"-#PMT=0 cut-HFShowerParam::"
245 <<
"getHits: npmt = " << npmt;
248 }
else if (npmt > 24) {
253 edm::LogInfo(
"HFShower") <<
"-SHORT cut-HFShowerParam::"
254 <<
"getHits:zMin=" << gpar[0];
260 edm::LogInfo(
"HFShower") <<
"HFShowerParam: npmt " << npmt
261 <<
" zv " <<
std::abs(pe_effect.z())
262 <<
":" << gpar[4] <<
":" << zv <<
":"
263 << gpar[0] <<
" ok " << ok <<
" depth "
267 if (G4UniformRand() > 0.5) depth = 2;
268 if (depth == 2 && zv < gpar[0]) ok =
false;
272 double r1 = G4UniformRand();
274 edm::LogInfo(
"HFShower") <<
"HFShowerParam:Distance to PMT (" <<npmt
275 <<
") " << dist <<
", exclusion flag "
284 hit.
time = time + hitSL[
i].time;
293 if (hit.
depth == 1) {
297 }
else if (hit.
depth == 2) {
303 edm::LogInfo(
"HFShower") <<
"HFShowerParam: Hit at depth "
304 << hit.
depth <<
" with edep "
305 << hit.
edep <<
" Time " << hit.
time;
313 double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
320 edm::LogInfo(
"HFShower") <<
"HFShowerParam: Hit at depth 1 with edep "
321 << edep <<
" Time " << tSlice <<
":" << time
330 edm::LogInfo(
"HFShower") <<
"HFShowerParam: Hit at depth 2 with edep "
331 << edep <<
" Time " << tSlice <<
":" << time
337 for (
unsigned int ii=0; ii<hits.size(); ++ii) {
339 if (zv > 12790)
edm::LogInfo(
"HFShower")<<
"HFShowerParam: Abnormal hit along "
341 << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
342 <<
" at " << hits[ii].position <<
" zz "
343 << zv <<
" Edep " << edep <<
" due to "
344 <<track->GetDefinition()->GetParticleName()
345 <<
" time " << hit.
time;
349 track->SetTrackStatus(fStopAndKill);
350 G4TrackVector tv = *(aStep->GetSecondary());
351 for (
unsigned int kk=0;
kk<tv.size(); ++
kk) {
352 if (tv[
kk]->
GetVolume() == preStepPoint->GetPhysicalVolume())
353 tv[
kk]->SetTrackStatus(fStopAndKill);
357 edm::LogInfo(
"HFShower") <<
"HFShowerParam: getHits kill (" << kill
358 <<
") track " << track->GetTrackID()
359 <<
" at " << hitPoint
360 <<
" and deposit " << edep <<
" " << hits.size()
361 <<
" times" <<
" ZZ " << zz <<
" " <<
gpar[0];
372 LogDebug(
"HFShower") <<
"HFShowerParam:getDDDArray called for " << str;
380 const std::vector<double> & fvec = value.
doubles();
381 int nval = fvec.size();
384 <<
" bins " << nval <<
" < 2 ==> illegal";
385 throw cms::Exception(
"Unknown",
"HFShowerParam") <<
"nval < 2 for array "
390 edm::LogError(
"HFShower") <<
"HFShowerParam : cannot get array " << str;
391 throw cms::Exception(
"Unknown",
"HFShowerParam") <<
"cannot get array "
std::vector< Hit > gfParameterization(G4Step *aStep, bool &ok, bool onlyLong=false)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
std::vector< Hit > getHits(G4Step *aStep, bool &ok, bool onlyLong=false)
void addFilter(const DDFilter &, log_op op=AND)
std::vector< Hit > getHits(G4Step *aStep)
static int PMTNumber(G4ThreeVector pe_effect)
double attLength(double lambda)
std::vector< double > gpar
HFShowerLibrary * showerLibrary
HFShowerParam(std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p)
type of data representation of DDCompactView
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
static int position[TOTALCHAMBERS][3]
double zShift(G4ThreeVector point, int depth, int fromEndAbs=0)
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &)
static const G4LogicalVolume * GetVolume(const std::string &name)
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
DDsvalues_type mergedSpecifics() const
void initRun(G4ParticleTable *theParticleTable)
void initRun(G4ParticleTable *)
bool firstChild()
set the current node to the first child ...
T * make() const
make new ROOT object
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Power< A, B >::type pow(const A &a, const B &b)
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.