FLSHitListAction.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file FLSHitListAction.cxx
3 /// \brief Use Geant4's user "hooks" to maintain a list of FLSHits
4 ///
5 /// \version $Id: FLSHitListAction.cxx,v 1.24 2012-09-20 21:45:42 greenc Exp $
6 /// \author brebel@fnal.gov
7 ////////////////////////////////////////////////////////////////////////
10 #include "Simulation/FLSHit.h"
11 #include "Simulation/FLSHitList.h"
12 
13 #include <vector>
14 #include <map>
15 
16 // ROOT includes
17 #include "TGeoMaterial.h"
18 #include <TGeoManager.h>
19 #include "TH2D.h"
20 
21 // G4 includes
22 #include "Geant4/G4Event.hh"
23 #include "Geant4/G4Track.hh"
24 #include "Geant4/G4Step.hh"
25 #include "Geant4/G4DynamicParticle.hh"
26 #include "Geant4/G4StepPoint.hh"
27 #include "Geant4/G4EnergyLossForExtrapolator.hh"
28 #include "Geant4/globals.hh"
29 
30 // ART includes
32 #include "cetlib_except/exception.h"
35 
36 #include "nug4/G4Base/UserActionFactory.h"
37 USERACTIONREG3(g4n,FLSHitListAction,g4n::FLSHitListAction)
38 
39 #include "GeometryObjects/CellGeo.h"
41 
42 #include <iostream>
43 #include <iomanip>
44 
45 namespace g4n
46 {
47 
48  //-------------------------------------------------------------
49  // Constructor.
51  fFLSHit(nullptr),
52  fFLSHitList(nullptr),
53  fEnergyCut(0),
54  fBirksConstant(0.0),
55  fChouConstant(0.0),
56  fIsRockVetoUsed(false),
57  fWavelengthLo(0),
58  fWavelengthHi(0),
59  fEnergyExtrapolator(nullptr),
60  fPerformFLSHitListSuppression(false),
61  fCerenkovPhotonsVBeta(0)
62  {
63  // Create the particle list that we'll (re-)use during the course
64  // of the Geant4 simulation.
65  fFLSHit = new sim::FLSHit();
67 
68  }
69 
70  //-------------------------------------------------------------
71  // Constructor.
72  FLSHitListAction::FLSHitListAction(const double min_energy, const bool is_rockveto_used) :
73  fFLSHit(nullptr),
74  fFLSHitList(nullptr),
75  fEnergyCut(min_energy*CLHEP::GeV),
76  fBirksConstant(0.0),
77  fChouConstant(0.0),
78  fIsRockVetoUsed(is_rockveto_used),
79  fEnergyExtrapolator(nullptr),
82  {
83  // Create the particle list that we'll (re-)use during the course
84  // of the Geant4 simulation.
85  fFLSHit = new sim::FLSHit();
87 
88  // Define the EnergyLossExtrapolator class if rock veto is used
89  if(is_rockveto_used)fEnergyExtrapolator = new G4EnergyLossForExtrapolator(0);
90 
91  }
92 
93  //-------------------------------------------------------------
94  // Destructor.
96  {
97  // Delete anything that we created with "new'.
98  delete fFLSHitList;
99  delete fFLSHit;
101  }
102 
103  //-------------------------------------------------------------
105  {
106  fEnergyCut = pset.get< double >("G4EnergyThreshold")*CLHEP::GeV;
107  fBirksConstant = pset.get< G4double >("BirksConstant")*CLHEP::cm/CLHEP::MeV;
108  fChouConstant = pset.get< G4double >("ChouConstant")*CLHEP::cm2/CLHEP::MeV/CLHEP::MeV;
109  fIsRockVetoUsed = pset.get< bool >("G4CheckRockVeto");
110  fPerformFLSHitListSuppression = pset.get< bool >("PerformFLSHitListSuppression");
111  fCauchyPars = pset.get< std::vector<double> >("CauchyPars");
112  fWavelengthLo = pset.get< G4double >("WavelengthLo")*CLHEP::nm;
113  fWavelengthHi = pset.get< G4double >("WavelengthHi")*CLHEP::nm;
114 
115  //units of c_i from fCauchyPars is (nm)^(2i)
116  G4double nm2 = CLHEP::nm*CLHEP::nm;
117  G4double parUnit = 1.0;
118  for (auto&& rPar : fCauchyPars ) {
119  rPar *= parUnit;
120  parUnit *= nm2;
121  }
122 
124 
125  }
126 
127  //-------------------------------------------------------------
129  {
130  mf::LogInfo("FLSHitListAction")
131  << "FLSHitListAction::PrintConfig \n"
132  << " EnergyCut " << fEnergyCut << "\n"
133  << " BirksConstant " << fBirksConstant*CLHEP::MeV/CLHEP::cm << "\n"
134  << " ChouConstant " << fChouConstant*CLHEP::MeV*CLHEP::MeV/CLHEP::cm2 << "\n"
135  << " IsRockVetoUsed " << (fIsRockVetoUsed?"true":"false") << "\n"
136  << " WavelengthLo " << fWavelengthLo/CLHEP::nm << "\n"
137  << " WavelengthHi " << fWavelengthHi/CLHEP::nm << "\n";
138  //figure out how to print CauchyPars
139  }
140 
141  //-------------------------------------------------------------
143  {
144  // Clear any previous particle information.
145  fFLSHitList->fHits.clear();
146  }
147 
148  //-------------------------------------------------------------
149  // Create our initial sim::FLSHit object and add it to the sim::FLSHitList.
151  {
152  // Clear FLS hits and make some other initializations
154  fIsRockVetoChecked = false;
155  fFLSHit->Clear();
157 
158  LOG_DEBUG("FLSHitListAction") << "pretracking step with track id: "
160  }
161 
162  //-------------------------------------------------------------
163  void FLSHitListAction::PostTrackingAction( const G4Track* /*track*/)
164  {
165  LOG_DEBUG("FLSHitListAction") << "done tracking for this g4 track, now "
166  << fFLSHitList->fHits.size() << " hits";
167  }
168 
169  double FLSHitListAction::GetRIndex(G4double wavelength)
170  {
171  double answer(0.0);
172  double xval(1.0);
173  G4double wavelength2 = wavelength*wavelength;
174  for (auto rPar : fCauchyPars) {
175  answer += rPar*xval;
176  xval /= wavelength2;
177  }
178  return answer;
179  }
180 
182  {
183  //check to see if fCerenkovPhotonsVBeta was ever initialized
184  bool firstTime = true;
185  if (fCerenkovPhotonsVBeta) {
186  delete fCerenkovPhotonsVBeta;
187  firstTime = false;
188  }
189  TH2D* hCerenkov2D = new TH2D("Cerenkov2D","", 2000, 0.0, 1.0, 100, fWavelengthLo, fWavelengthHi);
190  for (int iBeta = 1; iBeta <= hCerenkov2D->GetXaxis()->GetNbins(); ++iBeta) {
191  for (int iLambda = 1; iLambda <= hCerenkov2D->GetYaxis()->GetNbins(); ++iLambda) {
192  double beta = hCerenkov2D->GetXaxis()->GetBinCenter(iBeta);
193  G4double lambda = hCerenkov2D->GetYaxis()->GetBinCenter(iLambda);
194  G4double delLambda = hCerenkov2D->GetYaxis()->GetBinWidth(iLambda);
195  double rindex = GetRIndex(lambda);
196  double nPho = 0.0459/(lambda/CLHEP::cm*lambda/CLHEP::cm)*(1.0 - 1.0/pow(beta*rindex,2))*delLambda/CLHEP::cm;
197  if (nPho < 0.0) nPho = 0.0;
198  hCerenkov2D->SetBinContent(iBeta, iLambda, nPho);
199  }
200  }
201 
202  fCerenkovPhotonsVBeta = hCerenkov2D->ProjectionX("hCerenkov1D");
203 
204  if (firstTime) {
206  TH1D* hCerenkov = tfs->make<TH1D>("hCerenkov", ";#beta;Cerenkov Photons/cm", 2000, 0.0, 1.0);
207  for (int iBeta = 1; iBeta <= hCerenkov2D->GetXaxis()->GetNbins(); ++iBeta) {
208  hCerenkov->SetBinContent(iBeta, fCerenkovPhotonsVBeta->GetBinContent(iBeta));
209  }
210  }
211 
212  delete hCerenkov2D;
213  return;
214  }
215 
216  double FLSHitListAction::GetCerenkovPhotons(double beta, G4double pathlength, G4double charge)
217  {
218  double charge2 = std::pow( charge/CLHEP::eplus, 2 );
219  if (charge2 > 0) {
220  return fCerenkovPhotonsVBeta->Interpolate(beta)*charge2*(pathlength/CLHEP::cm);
221  }
222  else return 0.0;
223  }
224 
225  //-------------------------------------------------------------
226  // With every step, add to the particle's trajectory.
228  {
229  //mf::LogInfo("FLSHitListAction") << "FLSHitListAction::SteppingAction";
230 
231  /// Get the pointer to the track
232  G4Track *track = step->GetTrack();
233 
234  // Checking the RockVeto. This is to stop tracking low energy charged
235  // particles in the rock, since they are not going to produce any hits
236  // in the detector. Check it only for the first step.
237  // The reason, this function appears in the Stepping and not in the PreTrackingAction,
238  // where it is supposed to be is kind of stupid.
239  // The input step and track pointers in the Actions must be declared as const.
240  // Hence we cannot modify these objects. However step->GetTrack() returns non-const
241  // value (thank God) and thus we can modify it.
243 
244  // Check the rock veto
245  if(IsRockVetoed(track)){
246  // Kill the track
247  track->SetTrackStatus(fKillTrackAndSecondaries);
248  return;
249  }
250 
251  fIsRockVetoChecked = true;//! If it's already checked, don't need to check it on the next step
252  }// end of rock veto
253 
254 
255  // make sure the FLSHit object has been initialized
256  if(fFLSHit == 0){
257  mf::LogWarning("FLSHitListAction") << "WARNING - FLSHit object not initialized in FLSHitListAction";
258  return;
259  }
260 
261  const CLHEP::Hep3Vector &pos0 = step->GetPreStepPoint()->GetPosition(); // Start of the step
262  const CLHEP::Hep3Vector &pos = track->GetPosition(); // End of the step
263  const CLHEP::Hep3Vector &mom = track->GetMomentum();
264 
265  LOG_DEBUG("FLSHitListAction") << "flshits momentum = "
266  << mom.x() << " " << mom.y() << " " << mom.z()
267  << " " << mom.mag();
268 
269  double tpos0[3] = {pos0.x()/CLHEP::cm, pos0.y()/CLHEP::cm, pos0.z()/CLHEP::cm}; ///< Start of the step
270  double tpos1[3] = {pos.x()/CLHEP::cm , pos.y()/CLHEP::cm , pos.z()/CLHEP::cm}; ///< End of the step
271  double tpos2[3] = {0.}; ///< Some temporary vector that we'll make use of
272 
273  // If it's a null step, don't use it. Otherwise it may induce an additional
274  // FLSHit, which is wrong
275  if(tpos0[0]==tpos1[0] &&
276  tpos0[1]==tpos1[1] &&
277  tpos0[2]==tpos1[2])return;
278 
279  // If DetectorBigBox is used, we look if the particle is inside the box and
280  // going to be outside of the box. If so, stop tracking it.
281  if(fGeo->isDetectorBigBoxUsed()){
282  // tpos1 is the point after the step
283  bool isParticleGoingToBeOutsideDetectorBigBox = fGeo->isInsideDetectorBigBox(tpos1[0], tpos1[1], tpos1[2]);
284 
285  if(fIsParticleInsideDetectorBigBox && !isParticleGoingToBeOutsideDetectorBigBox)
286  track->SetTrackStatus(fStopAndKill);
287  else
288  fIsParticleInsideDetectorBigBox = isParticleGoingToBeOutsideDetectorBigBox;
289  }// end of check if the detectorbigbox used
290 
291 
292  //check that we are in the correct material to record a hit - ie scintillator
293  std::string material = track->GetMaterial()->GetName();
294  if(material.compare("Scintillator") != 0 ) {
295  return;
296  }
297 
298  //check if the current cell pointer is empty, that means
299  //this is the first step for this track
300  double dcos[3] = {mom.x()/mom.mag(), mom.y()/mom.mag(), mom.z()/mom.mag()};
301  geo::CellUniqueId cellid = fGeo->CellId(tpos1[0], tpos1[1], tpos1[2],
302  dcos[0], dcos[1], dcos[2],
303  step->GetStepLength()/CLHEP::cm);
304 
305  if(cellid == 0) cellid = fGeo->CellId(tpos0[0], tpos0[1], tpos0[2],
306  dcos[0], dcos[1], dcos[2],
307  step->GetStepLength()/CLHEP::cm);
308 
309  LOG_DEBUG("FSHitListAction") << "attempting to get cell id pos is "
310  << tpos0[0] << ", " << tpos0[1] << ", " << tpos0[2]
311  << " result is " << cellid;
312 
313  /// Get density for Birks' Law calculation
314  TGeoManager *geomanager = fGeo->ROOTGeoManager();
315  if ( ! geomanager ) {
316  throw cet::exception("NoTGeoManager")
317  << "FLSHitListAction: no TGeoManager given by fGeo. "
318  << "I need this to get the scintillator density. Failing.\n"
319  << __FILE__ << ":" << __LINE__ << "\n";
320  }
321  TGeoMaterial *mat = geomanager->GetMaterial("Scintillator");
322  if ( ! mat ) {
323  throw cet::exception("NoTGeoMaterial")
324  << "FLSHitListAction: no TGeoMaterial names 'Scintillator' found. "
325  << "I need this to get the scintillator density. Failing.\n"
326  << __FILE__ << ":" << __LINE__ << "\n";
327  }
328 
329  /// Setting Birks Constants
330  /// \todo We should not be setting Birks Constant at each Geant4 step.
331  /// 1. rho could be calculated in the very beginning (material is allways Scintillator) once
332  /// 2. SetBirksConstant may probably also be done once somehow,
333  /// but for sure we could do it in PreTrackingAction
334  {
335  const double rho = mat->GetDensity();
336  track->GetMaterial()->GetIonisation()->SetBirksConstant( (fBirksConstant/rho));
337  }// end of setting Birks Constants
338 
339  /// Getting Energy depositions
340  const double edep = step->GetTotalEnergyDeposit()/CLHEP::GeV;
341  fNovaG4EmSaturation.SetChouConstant(fChouConstant);
342  const double edepBirks = fNovaG4EmSaturation.VisibleEnergyDeposition(step)/CLHEP::GeV;
343 
344  /// Cerenkov photons
345  double cerenkov = GetCerenkovPhotons(0.5*(step->GetPreStepPoint()->GetBeta()+step->GetPostStepPoint()->GetBeta()), track->GetStepLength(), track->GetDynamicParticle()->GetCharge());
346 
347  /// Protect against geometry skew between G4/ROOT geometries
348  /// if cellid returns 0 it means that the ROOT geometry couldnt
349  /// find the same cell as the G4 Geometry - check the prestep point
350  /// volume name against the current volume and if they are the
351  /// same name then just reuse the current cell id
352  if(cellid == 0
353  && step->GetPreStepPoint()->GetPhysicalVolume()->GetName() == track->GetVolume()->GetName()
354  && fFLSHit->GetCellUniqueId() > 0
355  ) cellid = fFLSHit->GetCellUniqueId();
356 
357  if(cellid == 0){
358  /// oh dear, can't find which cell we are in with the ROOT geometry and
359  /// this is a brand new track with the fFLSHit having been reset.
360  /// add it to the flshit in the hopes that the next step will be in the same
361  /// cell
362 
363  mf::LogWarning("CantFindCell") << "Can't find cell in ROOT Geometry, first energy deposition"
364  << " dE = " << step->GetTotalEnergyDeposit()/CLHEP::GeV
365  << " (GeV) for this track will be added to first located hit.";
366 
367  fFLSHit->AddEdep( edep );
368  fFLSHit->AddEdepBirks( edepBirks );
369  fFLSHit->AddNCerenkov( cerenkov );
370  return;
371  }// end of check if the cell is not found
372 
373 
374  /// plane and cell number holders used below
375  int ip = 0, ic = 0;
376  /// Get local positions for the start of the step (tpos0) and store them in tpos2
377  /// Also, get ip and ic
378  fGeo->IdToCell(cellid, &ip, &ic)->WorldToLocal(tpos0, tpos2);
379 
380  //check to see if we are in a new cell with this step
381  if( fFLSHit->GetCellUniqueId() != cellid ){
382 
383  // the first step in this volume - clear the hit
384  fFLSHit->Clear();
385 
386  // Get the position and time the particle first enters
387  // the volume, as well as the pdg code. that information is
388  // obtained from the G4Track object that corresponds to this step
389  // account for the fact that we use cm, ns, GeV rather than the G4 defaults
390  fFLSHit->SetPDG ( track->GetDefinition()->GetPDGEncoding() );
392 
393  //figure out the geometry information for the nova geometry
394  fFLSHit->SetCellUniqueId(cellid);
395 
396  fFLSHit->SetPlane(ip);
397  fFLSHit->SetCell(ic);
398  }//end if first step in the volume
399 
400  fFLSHit->AddEdep( edep );
401  fFLSHit->AddEdepBirks( edepBirks );
402  fFLSHit->AddNCerenkov( cerenkov );
403  /// Add position
404  fFLSHit->AddPos(tpos2[0], tpos2[1], tpos2[2], (double)step->GetPreStepPoint()->GetGlobalTime()/CLHEP::ns, step->GetPreStepPoint()->GetKineticEnergy() / CLHEP::GeV);
405 
406  /// maybe this is the last step in the volume - either the particle exits,
407  /// disappears, or simply stops
408  if(track->GetVolume()->GetName() != step->GetPostStepPoint()->GetPhysicalVolume()->GetName()
409  || track->GetKineticEnergy() == 0.){
410 
411  //going into a new volume with the next step or no kinetic energy left for this particle
412  if(fFLSHit->GetEdep() > 0.){
413  /// Get local positions for the end of the step (tpos1) and store them in tpos2
414  fGeo->IdToCell(fFLSHit->GetCellUniqueId(), &ip, &ic)->WorldToLocal(tpos1, tpos2);
415 
416  fFLSHit->AddPos(tpos2[0], tpos2[1], tpos2[2], (double)step->GetPostStepPoint()->GetGlobalTime()/CLHEP::ns, step->GetPostStepPoint()->GetKineticEnergy() / CLHEP::GeV);
417 
418  // Fill the FLSHitList
419  fFLSHitList->fHits.push_back( *fFLSHit );
420 
421  /// Zero out ID now since it's already recorded
422  /// It needs to be zeroed out, because next step
423  /// may arrive to the same cell
424  /// and in that case we want a new FLSHit
426  }// end of check if the energy was deposited into the cell
427 
428  }//end if last hit in the volume
429 
430 
431  }// end of FLSHitListAction::SteppingAction
432 
433 
434  //------------------------------------------------------------------------------
435  // There's one last thing to do: All the particles have their
436  // parent IDs set (in PostTrackingAction), but we haven't set the
437  // daughters yet. That's done in this method.
439  {
440 
441  /// Perform the suppression if requested
443  }
444 
445 
446  //------------------------------------------------------------------------------
447  // Suppress fFLSHitList
449  /// collect the hits coming from the same track ids in the same cells
450  std::map< geo::CellUniqueId, std::vector<sim::FLSHit> > cellToHits;
451  std::map< geo::CellUniqueId, std::vector<sim::FLSHit> >::iterator cthitr;
452 
453  const int nflshits = fFLSHitList->fHits.size();
454  /// Loop over all FLSHits and construct a map of vectors of hits with CellID as a key
455  for(int h = 0; h < nflshits; ++h){
456  const sim::FLSHit& hit = fFLSHitList->fHits[h];
457  const geo::CellUniqueId cid = hit.GetCellUniqueId();
458  cthitr = cellToHits.find(cid);
459 
460  if( cthitr != cellToHits.end() )
461  (*cthitr).second.push_back(hit);
462  else{
463  std::vector<sim::FLSHit> vec;
464  vec.push_back(hit);
465  cellToHits[cid] = vec;
466  }
467  }
468 
469  // clear the fFLSHitList as all the hits are now in the map
470  fFLSHitList->fHits.clear();
471 
472  /// now loop over the map and collect hits from the
473  /// same track id in the same cell.
474  /// This loop is over all CellIDs
475  for(cthitr = cellToHits.begin(); cthitr != cellToHits.end(); ++cthitr){
476  /// This will be a map of vector FLSHits with the same (CellID and TrackID)
477  std::map< int, std::vector<sim::FLSHit> > trackToHits;
478  std::map< int, std::vector<sim::FLSHit> >::iterator tthitr;
479 
480  /// Loop over hits with the same Cell ID
481  /// and fill the map
482  for(size_t f = 0; f < (*cthitr).second.size(); ++f){
483  tthitr = trackToHits.find((*cthitr).second.at(f).GetTrackID());
484  if( tthitr != trackToHits.end() )
485  (*tthitr).second.push_back((*cthitr).second.at(f));
486  else{
487  std::vector<sim::FLSHit> vec;
488  vec.push_back((*cthitr).second.at(f));
489  trackToHits[(*cthitr).second.at(f).GetTrackID()] = vec;
490  }
491  }// end of loop over hits with the same Cell ID and filling the map
492 
493  // now that all the hits from this cell are sorted by track id,
494  // combine them into one hit for each id and add them to the fFLSHitList
495  // weight the various values by the energy deposited
496  for(tthitr = trackToHits.begin(); tthitr != trackToHits.end(); ++tthitr){
497  const int trackId = (*tthitr).first;
498  const std::vector<sim::FLSHit>& vec = (*tthitr).second;
499 
500  /// Construct the empty FLSHit
501  sim::FLSHit out_hit;
502 
503  /// Set the CellID, plane etc from the first FLSHit (these entries are identical within this loop)
504  out_hit.SetCellUniqueId(vec[0].GetCellUniqueId());
505  out_hit.SetPlane (vec[0].GetPlaneID());
506  out_hit.SetCell (vec[0].GetCellID());
507  out_hit.SetPDG (vec[0].GetPDG());
508  out_hit.SetTrackId (trackId);
509 
510  float edep_total = 0.0;
511  float edepbirks_total = 0.0;
512  float ncerenkov_total = 0.0;
513 
514 
515  float minT(vec.front().GetEntryT()), maxT(vec.back().GetExitT());
516  float entry_x(vec.front().GetEntryX()), exit_x(vec.back().GetExitX());
517  float entry_y(vec.front().GetEntryY()), exit_y(vec.back().GetExitY());
518  float entry_z(vec.front().GetEntryZ()), exit_z(vec.back().GetExitZ());
519  float entry_t(vec.front().GetEntryT()), exit_t(vec.back().GetExitT());
520  float entry_e(vec.front().GetEntryEnergy()), exit_e(vec.back().GetExitEnergy());
521 
522  /// Loop over FLSHits with the same CellID and TrackID
523  for(size_t f = 0; f < vec.size(); ++f){
524  const sim::FLSHit& h = vec[f];
525  edep_total += h.GetEdep();
526  edepbirks_total += h.GetEdepBirks();
527  ncerenkov_total += h.GetNCerenkov();
528 
529  if( h.GetEntryT() < minT ){
530  //std::cout << "out of order min" << std::endl;
531  minT = h.GetEntryT();
532  entry_x = h.GetEntryX();
533  entry_y = h.GetEntryY();
534  entry_z = h.GetEntryZ();
535  entry_e = h.GetEntryEnergy();
536  }
537  if( h.GetExitT() > maxT ){
538  //std::cout << "out of order max" << std::endl;
539  maxT = h.GetExitT();
540  exit_x = h.GetExitX();
541  exit_y = h.GetExitY();
542  exit_z = h.GetExitZ();
543  exit_e = h.GetExitEnergy();
544  }
545  }// end of the loop over FLSHits with the same Cell ID and TrackID
546 
547  // New hit with earliest entry point and latest exit point
548  // Yes, there are problems with this, but it's best for
549  // muon, proton, pion, kaon... _tracks_.
550  out_hit.AddPos(entry_x, entry_y, entry_z, entry_t, entry_e);
551  out_hit.AddPos(exit_x , exit_y , exit_z , exit_t , exit_e);
552 
553  /// Set the total edep, summed over all FLSHits with the same CellID and TrackID
554  out_hit.AddEdep(edep_total);
555  out_hit.AddEdepBirks(edepbirks_total);
556  out_hit.AddNCerenkov(ncerenkov_total);
557  // add the hit to fFLSHitList
558  fFLSHitList->fHits.push_back(out_hit);
559 
560  } // end loop over track id to hit map
561  }// end loop over cell id to hit map
562 
563  return true;
564  }
565 
566  //------------------------------------------------------------------------------
567  // Checking if particle is in the rock and that it's not going to reach the detector.
568  // If that's the case, veto it (kill the track and it's secondaries)
569  bool FLSHitListAction::IsRockVetoed(const G4Track* track) {
570 
571  //G4double energy = track->GetKineticEnergy();
572  // No veto'ing for uncharged particles
573 
574  if (! (track->GetDefinition()->GetPDGCharge()) ) return false;
575 
576  G4String material = track->GetMaterial()->GetName();
577  if(material != "Granite") return false;
578 
579  const CLHEP::Hep3Vector &pos = track->GetStep()->GetPreStepPoint()->GetPosition();
580 
581  double tpos[3] = {pos.x()/CLHEP::cm, pos.y()/CLHEP::cm, pos.z()/CLHEP::cm};
582  double tdir[3] = {track->GetMomentumDirection().getX(),
583  track->GetMomentumDirection().getY(),
584  track->GetMomentumDirection().getZ()};
585 
586 
587  gGeoManager->InitTrack(tpos, tdir);
588 
589  gGeoManager->FindNextBoundary();
590 
591  // Calculate the energy in GeV (I think)
592  double energyAfterStep = fEnergyExtrapolator->EnergyAfterStep(track->GetKineticEnergy()/CLHEP::MeV,
593  gGeoManager->GetStep() / CLHEP::cm*CLHEP::mm,
594  track->GetMaterial(),
595  track->GetDefinition());
596 
597  return (energyAfterStep < fEnergyCut);
598 
599  }
600 
601  //------------------------------------------------------------------------------
602  // For now, it's obsolete
604 
606 
607  double x1, x2;
608  double y1, y2;
609  double z1, z2;
610  fGeo->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
611 
612  std::cout<<"World X "<<x1<<" "<<x2<<"\n";
613  std::cout<<"World Y "<<y1<<" "<<y2<<"\n";
614  std::cout<<"World Z "<<z1<<" "<<z2<<"\n";
615 
616  fGeo->DetectorEnclosureBox(&x1, &x2, &y1, &y2, &z1, &z2);
617 
618  std::cout<<"Detector X "<<x1<<" "<<x2<<"\n";
619  std::cout<<"Detector Y "<<y1<<" "<<y2<<"\n";
620  std::cout<<"Detector Z "<<z1<<" "<<z2<<"\n";
621 
622  return false;
623 
624  const CLHEP::Hep3Vector &pos = track->GetStep()->GetPreStepPoint()->GetPosition();
625 
626  double tpos[3] = {pos.x()/CLHEP::cm, pos.y()/CLHEP::cm, pos.z()/CLHEP::cm};
627  double tdir[3] = {track->GetMomentumDirection().getX(),
628  track->GetMomentumDirection().getY(),
629  track->GetMomentumDirection().getZ()};
630 
631  gGeoManager->InitTrack(tpos, tdir);
632  gGeoManager->FindNextBoundary();
633 
634  G4double steplength = gGeoManager->GetStep() / CLHEP::cm*CLHEP::mm;
635 
636  //Calculate energy in GeV (I think)
637  double energyAfterStep = fEnergyExtrapolator->EnergyAfterStep(track->GetKineticEnergy()/CLHEP::MeV,
638  steplength,
639  track->GetMaterial(),
640  track->GetDefinition());
641 
642  return (energyAfterStep < fEnergyCut);
643  }
644 
645 }//end namespace
geo::CellUniqueId GetCellUniqueId() const
Unique Cell ID.
Definition: FLSHit.h:41
void AddNCerenkov(const float ncerenkov)
Definition: FLSHit.h:107
TruthSlim – remove generated objects that don&#39;t contribute.
G4double fWavelengthLo
Low wavelength for integrating over Cerenkov light productions.
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
float GetEntryX() const
Entry point of the particle (position, time and energy)
Definition: FLSHit.h:48
void SetTrackId(const int trackid)
Definition: FLSHit.h:101
double GetRIndex(double wavelength)
bool isDetectorBigBoxUsed() const
Do we use the Detector Big Box.
Definition: GeometryBase.h:190
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
double GetCerenkovPhotons(double beta, G4double pathLength, G4double charge)
NovaG4EmSaturation fNovaG4EmSaturation
Determines the visible energy (after Birks suppression) Modified to include Birks-Chou.
void WorldToLocal(const double *local, double *world) const
Definition: CellGeo.cxx:100
Float_t y1[n_points_granero]
Definition: compare.C:5
void PreTrackingAction(const G4Track *)
bool fPerformFLSHitListSuppression
Do we want to perform FLSHitList suppression?
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
float GetEntryEnergy() const
Definition: FLSHit.h:52
Float_t x1[n_points_granero]
Definition: compare.C:5
float GetExitT() const
Definition: FLSHit.h:57
constexpr T pow(T x)
Definition: pow.h:75
Definition: event.h:19
void BeginOfEventAction(const G4Event *)
TGeoManager * ROOTGeoManager() const
float GetNCerenkov() const
Get total N Cerenkov photons.
Definition: FLSHit.h:35
static constexpr double ns
double z() const
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
TString ip
Definition: loadincs.C:5
TH1D * fCerenkovPhotonsVBeta
Precomputed Cerenkov Photons/cm vs. beta.
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
Use Geant4&#39;s user "hooks" to maintain a list of particles generated by Geant4.
Double_t beta
void SteppingAction(const G4Step *)
void PostTrackingAction(const G4Track *)
float GetExitX() const
Exit point of the particle (position, time and energy)
Definition: FLSHit.h:54
void SetCell(const int ic)
Definition: FLSHit.h:103
void SetCellUniqueId(const geo::CellUniqueId &id)
Definition: FLSHit.h:104
G4double fBirksConstant
constant used in Birks suppression
void EndOfEventAction(const G4Event *)
string material
Definition: eplot.py:19
static constexpr double cm
Definition: SystemOfUnits.h:99
list of energy deposits from Geant4
void SetPDG(const int pdg)
Set methods.
Definition: FLSHit.h:100
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
static const int GetCurrentTrackID()
T get(std::string const &key) const
Definition: ParameterSet.h:231
art::ServiceHandle< geo::Geometry > fGeo
static constexpr double MeV
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
static constexpr Double_t GeV
Definition: Munits.h:291
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
G4double fWavelengthHi
High wavelength for integrating over Cerenkov light productions.
Eigen::VectorXd vec
Double_t edep
Definition: macro.C:13
sim::FLSHitList * fFLSHitList
The accumulated information for hits in the event.
void PrintConfig(std::string const &opt)
void Config(fhicl::ParameterSet const &pset)
float GetEdepBirks() const
Get total Energy with Birks suppression deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:33
static constexpr double cm2
OStream cout
Definition: OStream.cxx:6
float GetEntryZ() const
Definition: FLSHit.h:50
float GetExitEnergy() const
Definition: FLSHit.h:58
float GetEntryT() const
Definition: FLSHit.h:51
double mag() const
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
void AddPos(const float x, const float y, const float z, const double t, const float energy)
Definition: FLSHit.h:108
static constexpr double nm
Definition: SystemOfUnits.h:92
T * make(ARGS...args) const
void AddEdep(const float edep)
Definition: FLSHit.h:105
void SetPlane(const int ip)
Definition: FLSHit.h:102
A vector of FLSHit from single neutrino interaction.
Definition: FLSHitList.h:13
Definition: structs.h:12
bool ParticleProjection(G4Track *)
static constexpr double mm
Definition: SystemOfUnits.h:95
std::vector< double > fCauchyPars
Coefficients of the Cauchy parameterization of the refractive index.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double x() const
void Clear()
Clear the FLS hit.
Definition: FLSHit.cxx:39
bool IsRockVetoed(const G4Track *)
bool fIsRockVetoUsed
Do we ise the RockVeto method?
void DetectorEnclosureBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
double y() const
float GetExitY() const
Definition: FLSHit.h:55
sim::FLSHit * fFLSHit
The information for a single hit.
float GetExitZ() const
Definition: FLSHit.h:56
const CellGeo * IdToCell(const CellUniqueId &id, int *iplane, int *icell) const
void AddEdepBirks(const float edepbirks)
Definition: FLSHit.h:106
bool fIsRockVetoChecked
Is rock veto already checked?
static constexpr double GeV
G4double fChouConstant
2nd order constant used in Birks suppression
bool fIsParticleInsideDetectorBigBox
Is the particle inside the Big Box?
static constexpr double eplus
Eigen::MatrixXd mat
G4EnergyLossForExtrapolator * fEnergyExtrapolator
Used in the RockVeto method.
bool supressFLSHitList()
Suppress fFLSHitList.
float GetEntryY() const
Definition: FLSHit.h:49
enum BeamMode string