CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HFShowerParam.cc
Go to the documentation of this file.
1 // File: HFShowerParam.cc
3 // Description: Parametrized version of HF hits
5 
10 
14 
15 #include "G4VPhysicalVolume.hh"
16 #include "G4Step.hh"
17 #include "G4Track.hh"
18 #include "G4NavigationHistory.hh"
19 #include "Randomize.hh"
20 
21 #include "CLHEP/Units/GlobalPhysicalConstants.h"
22 #include "CLHEP/Units/GlobalSystemOfUnits.h"
23 #include<iostream>
24 
25 //#define DebugLog
26 //#define mkdebug
27 
28 HFShowerParam::HFShowerParam(std::string & name, const DDCompactView & cpv,
29  edm::ParameterSet const & p) : showerLibrary(0),
30  fibre(0), gflash(0),
31  fillHisto(false) {
32  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
33  pePerGeV = m_HF.getParameter<double>("PEPerGeV");
34  trackEM = m_HF.getParameter<bool>("TrackEM");
35  bool useShowerLibrary = m_HF.getParameter<bool>("UseShowerLibrary");
36  bool useGflash = m_HF.getParameter<bool>("UseHFGflash");
37  edMin = m_HF.getParameter<double>("EminLibrary");
38  onlyLong = m_HF.getParameter<bool>("OnlyLong");
39  ref_index = m_HF.getParameter<double>("RefIndex");
40  double lambdaMean = m_HF.getParameter<double>("LambdaMean");
41  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
42  parametrizeLast = m_HF.getUntrackedParameter<bool>("ParametrizeLast",false);
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;
53 
54 #ifdef DebugLog
56  if (tfile.isAvailable()) {
57  fillHisto = true;
58  LogDebug("HFShower") << "HFShowerParam::Save histos in directory "
59  << "ProfileFromParam";
60  TFileDirectory showerDir = tfile->mkdir("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);
71  } else {
72  fillHisto = false;
73  LogDebug("HFShower") << "HFShowerParam::No file is available for "
74  << "saving histos so the flag is set to false";
75  }
76 #endif
77 
78  G4String attribute = "ReadOutName";
79  G4String value = name;
81  DDValue ddv(attribute,value,0);
83  DDFilteredView fv(cpv);
84  fv.addFilter(filter);
85  bool dodet = fv.firstChild();
86  if (dodet) {
88  //Special Geometry parameters
89  gpar = getDDDArray("gparHF",sv);
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";
94  } else {
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";
99  }
100 
101  if (useShowerLibrary) showerLibrary = new HFShowerLibrary(name, cpv, p);
102  if (useGflash) gflash = new HFGflash(p);
103  fibre = new HFFibre(name, cpv, p);
104  attLMeanInv = fibre->attLength(lambdaMean);
105  edm::LogInfo("HFShower") << "att. length used for (lambda=" << lambdaMean
106  << ") = " << 1/(attLMeanInv*cm) << " cm";
107 }
108 
110  if (fibre) delete fibre;
111  if (gflash) delete gflash;
112  if (showerLibrary) delete showerLibrary;
113 }
114 
115 void HFShowerParam::initRun(G4ParticleTable * theParticleTable) {
116  emPDG = theParticleTable->FindParticle("e-")->GetPDGEncoding();
117  epPDG = theParticleTable->FindParticle("e+")->GetPDGEncoding();
118  gammaPDG = theParticleTable->FindParticle("gamma")->GetPDGEncoding();
119 #ifdef DebugLog
120  edm::LogInfo("HFShower") << "HFShowerParam: Particle code for e- = " << emPDG
121  << " for e+ = " << epPDG << " for gamma = " << gammaPDG;
122 #endif
123  if (showerLibrary) showerLibrary->initRun(theParticleTable);
124 }
125 
126 std::vector<HFShowerParam::Hit> HFShowerParam::getHits(G4Step * aStep) {
127  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
128  G4Track * track = aStep->GetTrack();
129  G4ThreeVector hitPoint = preStepPoint->GetPosition();
130  G4int particleCode = track->GetDefinition()->GetPDGEncoding();
131  double zv = std::abs(hitPoint.z()) - gpar[4] - 0.5*gpar[1];
132  G4ThreeVector localPoint = G4ThreeVector(hitPoint.x(),hitPoint.y(),zv);
133 
134  double pin = (preStepPoint->GetTotalEnergy())/GeV;
135  double zint = hitPoint.z();
136  double zz = std::abs(zint) - gpar[4];
137 
138 #ifdef DebugLog
139  edm::LogInfo("HFShower") << "HFShowerParam: getHits "
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 "
146  << localPoint;
147 #endif
148  std::vector<HFShowerParam::Hit> hits;
150  hit.position = hitPoint;
151 
152  // look for other charged particles
153  bool other = false;
154  double pBeta = track->GetDynamicParticle()->GetTotalMomentum() / track->GetDynamicParticle()->GetTotalEnergy();
155  if (particleCode != emPDG && particleCode != epPDG && particleCode != gammaPDG ) {
156  if (track->GetDefinition()->GetPDGCharge() != 0 && pBeta > (1/ref_index) &&
157  aStep->GetTotalEnergyDeposit() > 0) other = true;
158  }
159 
160  // take only e+-/gamma/or special particles
161  if (particleCode == emPDG || particleCode == epPDG ||
162  particleCode == gammaPDG || other) {
163  // Leave out the last part
164  double edep = 0.;
165  bool kill = false;
166  if ((!trackEM) && ((zz<(gpar[1]-gpar[2])) || parametrizeLast) && (!other)){
167  edep = pin;
168  kill = true;
169  } else if (track->GetDefinition()->GetPDGCharge() != 0 &&
170  pBeta > (1/ref_index)) {
171  edep = (aStep->GetTotalEnergyDeposit())/GeV;
172  }
173  std::string path = "ShowerLibrary";
174 #ifdef DebugLog
175  edm::LogInfo("HFShower") << "HFShowerParam: getHits edep = " << edep
176  << ", Kill = " << kill << ", pin = " << pin
177  << ", edMin = " << edMin << " Other " << other;
178 #endif
179  if (edep > 0) {
180  if ((showerLibrary || gflash) && kill && pin > edMin && (!other)) {
181  if (showerLibrary) {
182  std::vector<HFShowerLibrary::Hit> hitSL = showerLibrary->getHits(aStep,kill,onlyLong);
183  for (unsigned int i=0; i<hitSL.size(); i++) {
184  bool ok = true;
185 #ifdef DebugLog
186  edm::LogInfo("HFShower")<<"HFShowerParam: getHits applyFidCut = "<<applyFidCut;
187 #endif
188  if (applyFidCut) { // @@ For showerlibrary no z-cut for Short (no z)
189  int npmt = HFFibreFiducial:: PMTNumber(hitSL[i].position);
190  if (npmt <= 0) ok = false;
191  }
192  if (ok) {
193  hit.position = hitSL[i].position;
194  hit.depth = hitSL[i].depth;
195  hit.time = hitSL[i].time;
196  hit.edep = 1;
197  hits.push_back(hit);
198 #ifdef DebugLog
199  if (fillHisto) {
200  hzv->Fill(zv);
201  em_long_sl->Fill(hit.position.z()/cm);
202  double sq = sqrt(pow(hit.position.x()/cm,2)+pow(hit.position.y()/cm,2));
203  double zp = hit.position.z()/cm;
204  if (hit.depth == 1) {
205  em_2d_1->Fill(zp, sq);
206  em_lateral_1->Fill(sq);
207  em_long_1->Fill(zp);
208  } else if (hit.depth == 2) {
209  em_2d_2->Fill(zp, sq);
210  em_lateral_2->Fill(sq);
211  em_long_2->Fill(zp);
212  }
213  }
214  edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth "
215  << hit.depth << " with edep " << hit.edep
216  << " Time " << hit.time;
217 #endif
218  }
219  }
220  } else { // GFlash clusters with known z
221  std::vector<HFGflash::Hit>hitSL=gflash->gfParameterization(aStep,kill, onlyLong);
222  for (unsigned int i=0; i<hitSL.size(); ++i) {
223  bool ok = true;
224  G4ThreeVector pe_effect(hitSL[i].position.x(), hitSL[i].position.y(),
225  hitSL[i].position.z());
226  double zv = std::abs(pe_effect.z()) - gpar[4];
227  //depth
228  int depth = 1;
229  int npmt = 0;
230  if (zv < 0. || zv > gpar[1]) {
231 #ifdef mkdebug
232  std::cout<<"-#Zcut-HFShowerParam::getHits:z="<<zv<<",m="<<gpar[1]<<std::endl;
233 #endif
234  ok = false;
235  }
236  if (ok && applyFidCut) {
237  npmt = HFFibreFiducial:: PMTNumber(pe_effect);
238 #ifdef DebugLog
239  edm::LogInfo("HFShower") << "HFShowerParam::getHits:#PMT= "
240  << npmt << ",z = " << zv;
241 #endif
242  if (npmt <= 0) {
243 #ifdef DebugLog
244  edm::LogInfo("HFShower") << "-#PMT=0 cut-HFShowerParam::"
245  << "getHits: npmt = " << npmt;
246 #endif
247  ok = false;
248  } else if (npmt > 24) { // a short fibre
249  if (zv > gpar[0]) {
250  depth = 2;
251  } else {
252 #ifdef DebugLog
253  edm::LogInfo("HFShower") << "-SHORT cut-HFShowerParam::"
254  << "getHits:zMin=" << gpar[0];
255 #endif
256  ok = false;
257  }
258  }
259 #ifdef DebugLog
260  edm::LogInfo("HFShower") << "HFShowerParam: npmt " << npmt
261  << " zv " << std::abs(pe_effect.z())
262  << ":" << gpar[4] << ":" << zv << ":"
263  << gpar[0] << " ok " << ok << " depth "
264  << depth;
265 #endif
266  } else {
267  if (G4UniformRand() > 0.5) depth = 2;
268  if (depth == 2 && zv < gpar[0]) ok = false;
269  }
270  //attenuation
271  double dist = fibre->zShift(localPoint,depth,0); // distance to PMT
272  double r1 = G4UniformRand();
273 #ifdef DebugLog
274  edm::LogInfo("HFShower") << "HFShowerParam:Distance to PMT (" <<npmt
275  << ") " << dist << ", exclusion flag "
276  << (r1 > exp(-attLMeanInv*zv));
277 #endif
278  if (r1 > exp(-attLMeanInv*dist)) ok = false;
279  if (ok) {
280  double time = fibre->tShift(localPoint,depth,0); //remaining part
281 
282  hit.position = hitSL[i].position;
283  hit.depth = depth;
284  hit.time = time + hitSL[i].time;
285  hit.edep = 1;
286  hits.push_back(hit);
287 #ifdef DebugLog
288  if (fillHisto) {
289  em_long_gflash->Fill(pe_effect.z()/cm, hitSL[i].edep);
290  hzv->Fill(zv);
291  double sq = sqrt(pow(hit.position.x()/cm,2)+pow(hit.position.y()/cm,2));
292  double zp = hit.position.z()/cm;
293  if (hit.depth == 1) {
294  em_2d_1->Fill(zp, sq);
295  em_lateral_1->Fill(s);
296  em_long_1->Fill(zp);
297  } else if (hit.depth == 2) {
298  em_2d_2->Fill(zp, sq);
299  em_lateral_2->Fill(sq);
300  em_long_2->Fill(zp);
301  }
302  }
303  edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth "
304  << hit.depth << " with edep "
305  << hit.edep << " Time " << hit.time;
306 #endif
307  }
308  }
309  }
310  } else {
311  path = "Rest";
312  edep *= pePerGeV;
313  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
314  double time = fibre->tShift(localPoint,1,0); // remaining part
315  hit.depth = 1;
316  hit.time = tSlice+time;
317  hit.edep = edep;
318  hits.push_back(hit);
319 #ifdef DebugLog
320  edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth 1 with edep "
321  << edep << " Time " << tSlice << ":" << time
322  << ":" << hit.time;
323 #endif
324  if (zz >= gpar[0]) {
325  time = fibre->tShift(localPoint,2,0);
326  hit.depth = 2;
327  hit.time = tSlice+time;
328  hits.push_back(hit);
329 #ifdef DebugLog
330  edm::LogInfo("HFShower") <<"HFShowerParam: Hit at depth 2 with edep "
331  << edep << " Time " << tSlice << ":" << time
332  << hit.time;
333 #endif
334  }
335  }
336 #ifdef DebugLog
337  for (unsigned int ii=0; ii<hits.size(); ++ii) {
338  double zv = std::abs(hits[ii].position.z());
339  if (zv > 12790) edm::LogInfo("HFShower")<< "HFShowerParam: Abnormal hit along "
340  << path << " in "
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;
346  }
347 #endif
348  if (kill) {
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);
354  }
355  }
356 #ifdef DebugLog
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];
362 #endif
363  }
364  }
365  return hits;
366 }
367 
368 std::vector<double> HFShowerParam::getDDDArray(const std::string & str,
369  const DDsvalues_type & sv)
370 {
371 #ifdef DebugLog
372  LogDebug("HFShower") << "HFShowerParam:getDDDArray called for " << str;
373 #endif
374  DDValue value(str);
375  if (DDfetch(&sv,value))
376  {
377 #ifdef DebugLog
378  LogDebug("HFShower") << value;
379 #endif
380  const std::vector<double> & fvec = value.doubles();
381  int nval = fvec.size();
382  if (nval < 2) {
383  edm::LogError("HFShower") << "HFShowerParam : # of " << str
384  << " bins " << nval << " < 2 ==> illegal";
385  throw cms::Exception("Unknown", "HFShowerParam") << "nval < 2 for array "
386  << str << "\n";
387  }
388  return fvec;
389  } else {
390  edm::LogError("HFShower") << "HFShowerParam : cannot get array " << str;
391  throw cms::Exception("Unknown", "HFShowerParam") << "cannot get array "
392  << str << "\n";
393  }
394 }
#define LogDebug(id)
std::vector< Hit > gfParameterization(G4Step *aStep, bool &ok, bool onlyLong=false)
Definition: HFGflash.cc:88
T getParameter(std::string const &) const
double attLMeanInv
Definition: HFShowerParam.h:54
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:134
bool parametrizeLast
Definition: HFShowerParam.h:55
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)
virtual ~HFShowerParam()
double attLength(double lambda)
Definition: HFFibre.cc:113
std::vector< double > gpar
Definition: HFShowerParam.h:57
#define abs(x)
Definition: mlp_lapack.h:159
HFShowerLibrary * showerLibrary
Definition: HFShowerParam.h:51
HFShowerParam(std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p)
type of data representation of DDCompactView
Definition: DDCompactView.h:77
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
TH1F * em_lateral_2
Definition: HFShowerParam.h:59
HFFibre * fibre
Definition: HFShowerParam.h:52
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
double zShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:143
list path
Definition: scaleCards.py:51
T sqrt(T t)
Definition: SSEVec.h:46
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
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::...
Definition: DDsvalues.h:19
bool isAvailable() const
Definition: Service.h:47
HFGflash * gflash
Definition: HFShowerParam.h:53
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &)
G4ThreeVector position
Definition: HFShowerParam.h:38
static const G4LogicalVolume * GetVolume(const std::string &name)
double ref_index
Definition: HFShowerParam.h:54
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
DDsvalues_type mergedSpecifics() const
TH1F * em_long_sl
Definition: HFShowerParam.h:60
void initRun(G4ParticleTable *theParticleTable)
TH1F * em_long_gflash
Definition: HFShowerParam.h:60
TH1F * em_long_1
Definition: HFShowerParam.h:59
void initRun(G4ParticleTable *)
TH1F * em_long_1_tuned
Definition: HFShowerParam.h:60
bool firstChild()
set the current node to the first child ...
TH1F * em_lateral_1
Definition: HFShowerParam.h:59
T * make() const
make new ROOT object
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
tuple cout
Definition: gather_cfg.py:121
TH1F * em_long_2
Definition: HFShowerParam.h:59
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37