CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CaloTrkProcessing.cc
Go to the documentation of this file.
1 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
2 #include "SimG4Core/Notification/interface/TrackWithHistory.h"
3 #include "SimG4Core/Notification/interface/TrackInformation.h"
4 
5 
6 #include "SimG4CMS/Calo/interface/CaloTrkProcessing.h"
7 
8 #include "SimG4Core/Application/interface/SimTrackManager.h"
9 #include "DetectorDescription/Core/interface/DDFilter.h"
10 #include "DetectorDescription/Core/interface/DDFilteredView.h"
11 #include "DetectorDescription/Core/interface/DDSolid.h"
12 #include "DetectorDescription/Core/interface/DDValue.h"
13 #include "FWCore/Utilities/interface/Exception.h"
14 
15 #include "G4EventManager.hh"
16 
17 #include "G4LogicalVolumeStore.hh"
18 #include "G4LogicalVolume.hh"
19 #include "G4Step.hh"
20 #include "G4Track.hh"
21 
22 #include "G4SystemOfUnits.hh"
23 
24 //#define DebugLog
25 
26 CaloTrkProcessing::CaloTrkProcessing(G4String name,
27  const DDCompactView & cpv,
28  const SensitiveDetectorCatalog & clg,
29  edm::ParameterSet const & p,
30  const SimTrackManager* manager) :
31  SensitiveCaloDetector(name, cpv, clg, p), lastTrackID(-1),
32  m_trackManager(manager) {
33 
34  //Initialise the parameter set
35  edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("CaloTrkProcessing");
36  testBeam = m_p.getParameter<bool>("TestBeam");
37  eMin = m_p.getParameter<double>("EminTrack")*MeV;
38  putHistory = m_p.getParameter<bool>("PutHistory");
39 
40  edm::LogInfo("CaloSim") << "CaloTrkProcessing: Initailised with TestBeam = "
41  << testBeam << " Emin = " << eMin << " MeV and"
42  << " History flag " << putHistory;
43 
44  //Get the names
45  G4String attribute = "ReadOutName";
46  DDSpecificsFilter filter;
47  DDValue ddv(attribute,name,0);
48  filter.setCriteria(ddv,DDCompOp::equals);
49  DDFilteredView fv(cpv);
50  fv.addFilter(filter);
51  fv.firstChild();
52  DDsvalues_type sv(fv.mergedSpecifics());
53 
54  G4String value = "Calorimeter";
55  std::vector<std::string> caloNames = getNames (value, sv);
56 #ifdef DebugLog
57  LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
58  for (unsigned int i=0; i<caloNames.size(); i++)
59  LogDebug("CaloSim") << " (" << i << ") " << caloNames[i];
60 #endif
61 
62  value = "Levels";
63  std::vector<double> levels = getNumbers (value, sv);
64 #ifdef DebugLog
65  LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
66  for (unsigned int i=0; i<levels.size(); i++)
67  LogDebug("CaloSim") << " (" << i << ") " << levels[i];
68 #endif
69 
70  value = "Neighbours";
71  std::vector<double> neighbours = getNumbers (value, sv);
72 #ifdef DebugLog
73  LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
74  for (unsigned int i=0; i<neighbours.size(); i++)
75  LogDebug("CaloSim") << " (" << i << ") " << neighbours[i];
76 #endif
77 
78  value = "Inside";
79  std::vector<std::string> insideNames = getNames (value, sv);
80 #ifdef DebugLog
81  LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
82  for (unsigned int i=0; i<insideNames.size(); i++)
83  LogDebug("CaloSim") << " (" << i << ") " << insideNames[i];
84 #endif
85 
86  value = "InsideLevel";
87  std::vector<double> insideLevel = getNumbers (value, sv);
88 #ifdef DebugLog
89  LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
90  for (unsigned int i=0; i<insideLevel.size(); i++)
91  LogDebug("CaloSim") << " (" << i << ") " << insideLevel[i];
92 #endif
93 
94  if (caloNames.size() < neighbours.size()) {
95  edm::LogError("CaloSim") << "CaloTrkProcessing: # of Calorimeter bins "
96  << caloNames.size() << " does not match with "
97  << neighbours.size() << " ==> illegal ";
98  throw cms::Exception("Unknown", "CaloTrkProcessing")
99  << "Calorimeter array size does not match with size of neighbours\n";
100  }
101 
102  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
103  std::vector<G4LogicalVolume *>::const_iterator lvcite;
104  int istart = 0;
105  for (unsigned int i=0; i<caloNames.size(); i++) {
106  G4LogicalVolume* lv = 0;
107  G4String name = caloNames[i];
108  int number = static_cast<int>(neighbours[i]);
109  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
110  if ((*lvcite)->GetName() == name) {
111  lv = (*lvcite);
112  break;
113  }
114  }
115  if (lv != 0) {
116  CaloTrkProcessing::Detector detector;
117  detector.name = name;
118  detector.lv = lv;
119  detector.level = static_cast<int>(levels[i]);
120  if (istart+number > (int)(insideNames.size())) {
121  edm::LogError("CaloSim") << "CaloTrkProcessing: # of InsideNames bins "
122  << insideNames.size() <<" too few compaerd to "
123  << istart+number << " requested ==> illegal ";
124  throw cms::Exception("Unknown", "CaloTrkProcessing")
125  << "InsideNames array size does not match with list of neighbours\n";
126  }
127  std::vector<std::string> inside;
128  std::vector<G4LogicalVolume*> insideLV;
129  std::vector<int> insideLevels;
130  for (int k = 0; k < number; k++) {
131  lv = 0;
132  name = insideNames[istart+k];
133  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
134  if ((*lvcite)->GetName() == name) {
135  lv = (*lvcite);
136  break;
137  }
138  inside.push_back(name);
139  insideLV.push_back(lv);
140  insideLevels.push_back(static_cast<int>(insideLevel[istart+k]));
141  }
142  detector.fromDets = inside;
143  detector.fromDetL = insideLV;
144  detector.fromLevels = insideLevels;
145  detectors.push_back(detector);
146  }
147  istart += number;
148  }
149 
150  edm::LogInfo("CaloSim") << "CaloTrkProcessing: with " << detectors.size()
151  << " calorimetric volumes";
152  for (unsigned int i=0; i<detectors.size(); i++) {
153  edm::LogInfo("CaloSim") << "CaloTrkProcessing: Calorimeter volume " << i
154  << " " << detectors[i].name << " LV "
155  << detectors[i].lv << " at level "
156  << detectors[i].level << " with "
157  << detectors[i].fromDets.size() << " neighbours";
158  for (unsigned int k=0; k<detectors[i].fromDets.size(); k++)
159  edm::LogInfo("CaloSim") << " Element " << k << " "
160  << detectors[i].fromDets[k] << " LV "
161  << detectors[i].fromDetL[k] << " at level "
162  << detectors[i].fromLevels[k];
163  }
164 }
165 
166 CaloTrkProcessing::~CaloTrkProcessing() {
167  edm::LogInfo("CaloSim") << "CaloTrkProcessing: Deleted";
168 }
169 
170 void CaloTrkProcessing::update(const BeginOfEvent * evt) {
171  lastTrackID = -1;
172 }
173 
174 void CaloTrkProcessing::update(const G4Step * aStep) {
175 
176  // define if you are at the surface of CALO
177 
178  G4Track* theTrack = aStep->GetTrack();
179  int id = theTrack->GetTrackID();
180 
181  TrackInformation* trkInfo = dynamic_cast<TrackInformation*>
182  (theTrack->GetUserInformation());
183 
184  if (trkInfo == 0) {
185  edm::LogError("CaloSim") << "CaloTrkProcessing: No trk info !!!! abort ";
186  throw cms::Exception("Unknown", "CaloTrkProcessing")
187  << "cannot get trkInfo for Track " << id << "\n";
188  }
189 
190  if (testBeam) {
191  if (trkInfo->getIDonCaloSurface() == 0) {
192 #ifdef DebugLog
193  LogDebug("CaloSim") << "CaloTrkProcessing set IDonCaloSurface to " << id
194  << " at step Number "
195  << theTrack->GetCurrentStepNumber();
196 #endif
197  trkInfo->setIDonCaloSurface(id,0,0,
198  theTrack->GetDefinition()->GetPDGEncoding(),
199  theTrack->GetMomentum().mag());
200  lastTrackID = id;
201  if (theTrack->GetKineticEnergy()/MeV > eMin)
202  trkInfo->putInHistory();
203  }
204  } else {
205  if (putHistory) {
206  trkInfo->putInHistory();
207  // trkInfo->setAncestor();
208  }
209 #ifdef DebugLog
210  LogDebug("CaloSim") << "CaloTrkProcessing Entered for " << id
211  << " at stepNumber "<< theTrack->GetCurrentStepNumber()
212  << " IDonCaloSur.. " << trkInfo->getIDonCaloSurface()
213  << " CaloCheck " << trkInfo->caloIDChecked();
214 #endif
215  if (trkInfo->getIDonCaloSurface() != 0) {
216  if (trkInfo->caloIDChecked() == false) {
217  G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
218  const G4VTouchable* post_touch = postStepPoint->GetTouchable();
219 
220  if (isItInside(post_touch, trkInfo->getIDCaloVolume(),
221  trkInfo->getIDLastVolume()) > 0) {
222  trkInfo->setIDonCaloSurface(0, -1, -1, 0, 0);
223  } else {
224  trkInfo->setCaloIDChecked(true);
225  }
226  }
227  } else {
228  G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
229  const G4VTouchable* post_touch = postStepPoint->GetTouchable();
230  int ical = isItCalo(post_touch);
231  if (ical >= 0) {
232  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
233  const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
234  int inside = isItInside(pre_touch, ical, -1);
235  if (inside >= 0 || (theTrack->GetCurrentStepNumber()==1)) {
236  trkInfo->setIDonCaloSurface(id,ical,inside,
237  theTrack->GetDefinition()->GetPDGEncoding(),
238  theTrack->GetMomentum().mag());
239  trkInfo->setCaloIDChecked(true);
240  lastTrackID = id;
241  if (theTrack->GetKineticEnergy()/MeV > eMin) trkInfo->putInHistory();
242 #ifdef DebugLog
243  LogDebug("CaloSim") <<"CaloTrkProcessing: set ID on Calo " << ical
244  << " surface (Inside " << inside << ") to "
245  << id << " of a Track with Kinetic Energy "
246  << theTrack->GetKineticEnergy()/MeV << " MeV";
247 #endif
248  }
249  }
250  }
251  }
252 }
253 
254 std::vector<std::string> CaloTrkProcessing::getNames(const G4String str,
255  const DDsvalues_type &sv){
256 
257 #ifdef DebugLog
258  LogDebug("CaloSim") << "CaloTrkProcessing::getNames called for " << str;
259 #endif
260  DDValue value(str);
261  if (DDfetch(&sv,value)) {
262 #ifdef DebugLog
263  LogDebug("CaloSim") << value;
264 #endif
265  const std::vector<std::string> & fvec = value.strings();
266  int nval = fvec.size();
267  if (nval < 1) {
268 
269  edm::LogError("CaloSim") << "CaloTrkProcessing: # of " << str
270  << " bins " << nval << " < 1 ==> illegal ";
271  throw cms::Exception("Unknown", "CaloTrkProcessing")
272  << "nval < 2 for array " << str << "\n";
273  }
274 
275  return fvec;
276  } else {
277  edm::LogError("CaloSim") << "CaloTrkProcessing: cannot get array " << str ;
278  throw cms::Exception("Unknown", "CaloTrkProcessing")
279  << "cannot get array " << str << "\n";
280  }
281 }
282 
283 std::vector<double> CaloTrkProcessing::getNumbers(const G4String str,
284  const DDsvalues_type &sv) {
285 
286 #ifdef DebugLog
287  LogDebug("CaloSim") << "CaloTrkProcessing::getNumbers called for " << str;
288 #endif
289  DDValue value(str);
290  if (DDfetch(&sv,value)) {
291 #ifdef DebugLog
292  LogDebug("CaloSim") << value;
293 #endif
294  const std::vector<double> & fvec = value.doubles();
295  int nval = fvec.size();
296  if (nval < 1) {
297  edm::LogError("CaloSim") << "CaloTrkProcessing: # of " << str
298  << " bins " << nval << " < 1 ==> illegal ";
299  throw cms::Exception("Unknown", "CaloTrkProcessing")
300  << "nval < 2 for array " << str << "\n";
301  }
302 
303  return fvec;
304  } else {
305  edm::LogError("CaloSim") << "CaloTrkProcessing: cannot get array " << str ;
306  throw cms::Exception("Unknown", "CaloTrkProcessing")
307  << "cannot get array " << str << "\n";
308  }
309 }
310 
311 int CaloTrkProcessing::isItCalo(const G4VTouchable* touch) {
312 
313  int lastLevel = -1;
314  G4LogicalVolume* lv=0;
315  for (unsigned int it=0; it < detectors.size(); it++) {
316  if (lastLevel != detectors[it].level) {
317  lastLevel = detectors[it].level;
318  lv = detLV(touch, lastLevel);
319 #ifdef DebugLog
320  std::string name1 = "Unknown";
321  if (lv != 0) name1 = lv->GetName();
322  LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
323  << " at Level " << lastLevel;
324  int levels = detLevels(touch);
325  if (levels > 0) {
326  G4String name2[20]; int copyno2[20];
327  detectorLevel(touch, levels, copyno2, name2);
328  for (int i2=0; i2<levels; i2++)
329  LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
330  << copyno2[i2];
331  }
332 #endif
333  }
334  bool ok = (lv == detectors[it].lv);
335  if (ok) return it;
336  }
337  return -1;
338 }
339 
340 int CaloTrkProcessing::isItInside(const G4VTouchable* touch, int idcal,
341  int idin) {
342  int lastLevel = -1;
343  G4LogicalVolume* lv=0;
344  int id1, id2;
345  if (idcal < 0) {id1 = 0; id2 = static_cast<int>(detectors.size());}
346  else {id1 = idcal; id2 = id1+1;}
347  for (int it1 = id1; it1 < id2; it1++) {
348  if (idin < 0) {
349  for (unsigned int it2 = 0; it2 < detectors[it1].fromDets.size(); it2++) {
350  if (lastLevel != detectors[it1].fromLevels[it2]) {
351  lastLevel = detectors[it1].fromLevels[it2];
352  lv = detLV(touch,lastLevel);
353 #ifdef DebugLog
354  std::string name1 = "Unknown";
355  if (lv != 0) name1 = lv->GetName();
356  LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
357  << " at Level " << lastLevel;
358  int levels = detLevels(touch);
359  if (levels > 0) {
360  G4String name2[20]; int copyno2[20];
361  detectorLevel(touch, levels, copyno2, name2);
362  for (int i2=0; i2<levels; i2++)
363  LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
364  << copyno2[i2];
365  }
366 #endif
367  }
368  bool ok = (lv == detectors[it1].fromDetL[it2]);
369  if (ok) return it2;
370  }
371  } else {
372  lastLevel = detectors[it1].fromLevels[idin];
373  lv = detLV(touch,lastLevel);
374 #ifdef DebugLog
375  std::string name1 = "Unknown";
376  if (lv != 0) name1 = lv->GetName();
377  LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
378  << " at Level " << lastLevel;
379  int levels = detLevels(touch);
380  if (levels > 0) {
381  G4String name2[20]; int copyno2[20];
382  detectorLevel(touch, levels, copyno2, name2);
383  for (int i2=0; i2<levels; i2++)
384  LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
385  << copyno2[i2];
386  }
387 #endif
388  bool ok = (lv == detectors[it1].fromDetL[idin]);
389  if (ok) return idin;
390  }
391  }
392  return -1;
393 }
394 
395 int CaloTrkProcessing::detLevels(const G4VTouchable* touch) const {
396 
397  //Return number of levels
398  if (touch)
399  return ((touch->GetHistoryDepth())+1);
400  else
401  return 0;
402 }
403 
404 G4LogicalVolume* CaloTrkProcessing::detLV(const G4VTouchable* touch,
405  int currentlevel) const {
406 
407  G4LogicalVolume* lv=0;
408  if (touch) {
409  int level = ((touch->GetHistoryDepth())+1);
410  if (level > 0 && level >= currentlevel) {
411  int ii = level - currentlevel;
412  lv = touch->GetVolume(ii)->GetLogicalVolume();
413  return lv;
414  }
415  }
416  return lv;
417 }
418 
419 void CaloTrkProcessing::detectorLevel(const G4VTouchable* touch, int& level,
420  int* copyno, G4String* name) const {
421 
422  static const std::string unknown("Unknown");
423  //Get name and copy numbers
424  if (level > 0) {
425  for (int ii = 0; ii < level; ii++) {
426  int i = level - ii - 1;
427  G4VPhysicalVolume* pv = touch->GetVolume(i);
428  if (pv != 0)
429  name[ii] = pv->GetName();
430  else
431  name[ii] = unknown;
432  copyno[ii] = touch->GetReplicaNumber(i);
433  }
434  }
435 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
CaloTrkProcessing::detLV
G4LogicalVolume * detLV(const G4VTouchable *, int) const
Definition: CaloTrkProcessing.cc:404
DDValue::doubles
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:139
CaloTrkProcessing::CaloTrkProcessing
CaloTrkProcessing(G4String aSDname, const DDCompactView &cpv, const SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *)
Definition: CaloTrkProcessing.cc:26
TrackInformation::getIDonCaloSurface
int getIDonCaloSurface() const
Definition: TrackInformation.h:36
CaloTrkProcessing::detLevels
int detLevels(const G4VTouchable *) const
Definition: CaloTrkProcessing.cc:395
DDFilteredView::addFilter
void addFilter(const DDFilter &, DDLogOp op=DDLogOp::AND)
Definition: DDFilteredView.cc:19
CaloTrkProcessing::isItCalo
int isItCalo(const G4VTouchable *)
Definition: CaloTrkProcessing.cc:311
alcazmumu_cfi.filter
tuple filter
Definition: alcazmumu_cfi.py:34
cuy.ii
int ii
Definition: cuy.py:588
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:77
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:80
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
CaloTrkProcessing::detectorLevel
void detectorLevel(const G4VTouchable *, int &, int *, G4String *) const
Definition: CaloTrkProcessing.cc:419
MeV
const double MeV
TrackInformation::getIDCaloVolume
int getIDCaloVolume() const
Definition: TrackInformation.h:40
CaloTrkProcessing::isItInside
int isItInside(const G4VTouchable *, int, int)
Definition: CaloTrkProcessing.cc:340
CaloTrkProcessing::Detector::fromLevels
std::vector< int > fromLevels
Definition: CaloTrkProcessing.h:58
CaloTrkProcessing::detectors
std::vector< Detector > detectors
Definition: CaloTrkProcessing.h:70
DDsvalues_type
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
CaloTrkProcessing::getNumbers
std::vector< double > getNumbers(G4String, const DDsvalues_type &)
Definition: CaloTrkProcessing.cc:283
CaloTrkProcessing::Detector::fromDetL
std::vector< G4LogicalVolume * > fromDetL
Definition: CaloTrkProcessing.h:57
SensitiveDetector::getNames
virtual std::vector< std::string > getNames()
Definition: SensitiveDetector.h:41
TrackInformation::getIDLastVolume
int getIDLastVolume() const
Definition: TrackInformation.h:41
CaloTrkProcessing::Detector::fromDets
std::vector< std::string > fromDets
Definition: CaloTrkProcessing.h:56
DDValue::strings
const std::vector< std::string > & strings() const
a reference to the std::string-valued values stored in the given instance of DDValue ...
Definition: DDValue.h:62
CaloTrkProcessing::update
void update(const BeginOfEvent *evt)
This routine will be called when the appropriate signal arrives.
Definition: CaloTrkProcessing.cc:170
TrackInformation::caloIDChecked
bool caloIDChecked() const
Definition: TrackInformation.h:42
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:56
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
TrackInformation::setCaloIDChecked
void setCaloIDChecked(bool f)
Definition: TrackInformation.h:43
TrackInformation::setIDonCaloSurface
void setIDonCaloSurface(int id, int ical, int last, int pdgID, double p)
Definition: TrackInformation.h:37
DDFilteredView::firstChild
bool firstChild()
set the current node to the first child ...
Definition: DDFilteredView.cc:110
testEve_cfg.level
tuple level
Definition: testEve_cfg.py:34
relval_steps.k
list k
Definition: relval_steps.py:1709
DDSpecificsFilter::setCriteria
void setCriteria(const DDValue &nameVal, DDCompOp, DDLogOp l=DDLogOp::AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:245
DDSpecificsFilter
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:32
Generated for CMSSW Reference Manual by doxygen 1.8.5