All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
SBNDMCFlash Class Reference
Inheritance diagram for SBNDMCFlash:

Public Member Functions

 SBNDMCFlash (fhicl::ParameterSet const &p)
 
 SBNDMCFlash (SBNDMCFlash const &)=delete
 
 SBNDMCFlash (SBNDMCFlash &&)=delete
 
SBNDMCFlash & operator= (SBNDMCFlash const &)=delete
 
SBNDMCFlash & operator= (SBNDMCFlash &&)=delete
 
void produce (art::Event &e) override
 

Private Member Functions

void GetFlashLocation (std::vector< double >, double &, double &, double &, double &)
 

Private Attributes

std::string _mctruth_label
 
std::string _trigger_label
 
std::string _simphot_label
 
std::string _simphot_insta
 
std::vector< std::string > _pd_to_use
 
int _tpc
 
float _qe_direct
 
float _qe_refl
 
float _window_length
 
float _pre_window
 
bool _debug
 
opdet::sbndPDMapAlg _pds_map
 
int _run
 
int _subrun
 
int _event
 
int _pe_total = 0
 
int _pe_scintillation = 0
 
int _pe_cherenkov = 0
 
std::vector< double > _cherenkov_time_v
 
std::vector< int > _cherenkov_pmt_v
 
std::vector< double > _scintillation_time_v
 
std::vector< int > _scintillation_pmt_v
 
TTree * _tree1
 

Detailed Description

Definition at line 46 of file SBNDMCFlash_module.cc.

Constructor & Destructor Documentation

SBNDMCFlash::SBNDMCFlash ( fhicl::ParameterSet const &  p)
explicit

Definition at line 92 of file SBNDMCFlash_module.cc.

93  : EDProducer{p} // ,
94  // More initializers here.
95 {
96  _mctruth_label = p.get<std::string>("MCTruthProduct", "generator");
97  _trigger_label = p.get<std::string>("TriggerProduct", "triggersim");
98  _simphot_label = p.get<std::string>("SimPhotProduct", "largeant");
99  _simphot_insta = p.get<std::string>("SimPhotProductInstance", "");
100  _pd_to_use = p.get<std::vector<std::string>>("PD", _pd_to_use);
101  _tpc = p.get<int>("TPC", -1);
102  _qe_direct = p.get<float>("QEDirect", 0.03);
103  _qe_refl = p.get<float>("QERefl", 0.03);
104  _window_length = p.get<float>("WindowLength", 8); // us
105  _pre_window = p.get<float>("PreWindow", 0.1); // us
106  _debug = p.get<bool>("DebugMode", false);
107 
108  art::ServiceHandle<art::TFileService> fs;
109  _tree1 = fs->make<TTree>("tree","");
110  _tree1->Branch("run", &_run, "run/I");
111  _tree1->Branch("subrun", &_subrun, "subrun/I");
112  _tree1->Branch("event", &_event, "event/I");
113  _tree1->Branch("pe_total", &_pe_total, "pe_total/I");
114  _tree1->Branch("pe_scintillation", &_pe_scintillation, "pe_scintillation/I");
115  _tree1->Branch("pe_cherenkov", &_pe_cherenkov, "pe_cherenkov/I");
116  _tree1->Branch("cherenkov_time_v", "std::vector<double>", &_cherenkov_time_v);
117  _tree1->Branch("cherenkov_pmt_v", "std::vector<int>", &_cherenkov_pmt_v);
118  _tree1->Branch("scintillation_time_v", "std::vector<double>", &_scintillation_time_v);
119  _tree1->Branch("scintillation_pmt_v", "std::vector<int>", &_scintillation_pmt_v);
120 
121  produces< std::vector<recob::OpFlash> >();
122 }
SBNDMCFlash::_scintillation_time_v
std::vector< double > _scintillation_time_v
Definition: SBNDMCFlash_module.cc:83
SBNDMCFlash::_simphot_label
std::string _simphot_label
Definition: SBNDMCFlash_module.cc:65
p
pdgs p
Definition: selectors.fcl:22
SBNDMCFlash::_scintillation_pmt_v
std::vector< int > _scintillation_pmt_v
Definition: SBNDMCFlash_module.cc:84
SBNDMCFlash::_mctruth_label
std::string _mctruth_label
Definition: SBNDMCFlash_module.cc:63
SBNDMCFlash::_pd_to_use
std::vector< std::string > _pd_to_use
Definition: SBNDMCFlash_module.cc:67
SBNDMCFlash::_cherenkov_pmt_v
std::vector< int > _cherenkov_pmt_v
Definition: SBNDMCFlash_module.cc:82
SBNDMCFlash::_cherenkov_time_v
std::vector< double > _cherenkov_time_v
Definition: SBNDMCFlash_module.cc:81
SBNDMCFlash::_simphot_insta
std::string _simphot_insta
Definition: SBNDMCFlash_module.cc:66
SBNDMCFlash::_trigger_label
std::string _trigger_label
Definition: SBNDMCFlash_module.cc:64
SBNDMCFlash::SBNDMCFlash ( SBNDMCFlash const &  )
delete
SBNDMCFlash::SBNDMCFlash ( SBNDMCFlash &&  )
delete

Member Function Documentation

void SBNDMCFlash::GetFlashLocation ( std::vector< double >  pePerOpChannel,
double &  Ycenter,
double &  Zcenter,
double &  Ywidth,
double &  Zwidth 
)
private

Definition at line 351 of file SBNDMCFlash_module.cc.

356 {
357 
358  // Reset variables
359  Ycenter = Zcenter = 0.;
360  Ywidth = Zwidth = -999.;
361  double totalPE = 0.;
362  double sumy = 0., sumz = 0., sumy2 = 0., sumz2 = 0.;
363 
364  for (unsigned int opch = 0; opch < pePerOpChannel.size(); opch++) {
365 
366  // Get physical detector location for this opChannel
367  double PMTxyz[3];
368  ::art::ServiceHandle<geo::Geometry> geo;
369  geo->OpDetGeoFromOpChannel(opch).GetCenter(PMTxyz);
370 
371  // Add up the position, weighting with PEs
372  sumy += pePerOpChannel[opch]*PMTxyz[1];
373  sumy2 += pePerOpChannel[opch]*PMTxyz[1]*PMTxyz[1];
374  sumz += pePerOpChannel[opch]*PMTxyz[2];
375  sumz2 += pePerOpChannel[opch]*PMTxyz[2]*PMTxyz[2];
376 
377  totalPE += pePerOpChannel[opch];
378  }
379 
380  Ycenter = sumy/totalPE;
381  Zcenter = sumz/totalPE;
382 
383  // This is just sqrt(<x^2> - <x>^2)
384  if ( (sumy2*totalPE - sumy*sumy) > 0. )
385  Ywidth = std::sqrt(sumy2*totalPE - sumy*sumy)/totalPE;
386 
387  if ( (sumz2*totalPE - sumz*sumz) > 0. )
388  Zwidth = std::sqrt(sumz2*totalPE - sumz*sumz)/totalPE;
389 }
SBNDMCFlash& SBNDMCFlash::operator= ( SBNDMCFlash const &  )
delete
SBNDMCFlash& SBNDMCFlash::operator= ( SBNDMCFlash &&  )
delete
void SBNDMCFlash::produce ( art::Event &  e)
override

Definition at line 124 of file SBNDMCFlash_module.cc.

125 {
126 
127  ::art::ServiceHandle<geo::Geometry> geo;
128  auto const *lar_prop = lar::providerFrom<detinfo::LArPropertiesService>();
129  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
130  auto const det_prop = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clock_data);
131 
132  _qe_direct = _qe_direct / lar_prop->ScintPreScale();
133  _qe_refl = _qe_refl / lar_prop->ScintPreScale();
134 
135  _run = e.id().run();
136  _subrun = e.id().subRun();
137  _event = e.id().event();
138 
139  // produce OpFlash data-product to be filled within module
140  std::unique_ptr< std::vector<recob::OpFlash> > opflashes(new std::vector<recob::OpFlash>);
141 
142  if (e.isRealData()) {
143  e.put(std::move(opflashes));
144  return;
145  }
146 
147  // art::Handle<std::vector<raw::Trigger> > evt_trigger_h;
148  // e.getByLabel(_trigger_label,evt_trigger_h);
149 
150  // if( !evt_trigger_h.isValid() || evt_trigger_h->empty() ) {
151  // std::cout << "[NeutrinoMCFlash] Trigger product is not valid or empty." << std::endl;
152  // e.put(std::move(opflashes));
153  // return;
154  // }
155 
156  art::Handle<std::vector<simb::MCTruth> > evt_mctruth_h;
157  e.getByLabel(_mctruth_label,evt_mctruth_h);
158 
159  if( !evt_mctruth_h.isValid() || evt_mctruth_h->empty() ) {
160  std::cerr << "[NeutrinoMCFlash] MCTruth product is not valid or empty." << std::endl;
161  e.put(std::move(opflashes));
162  return;
163  }
164 
165  // art::Handle<std::vector<sim::SimPhotons> > evt_simphot_h;
166  // e.getByLabel(_simphot_label,evt_simphot_h);
167 
168  // if( !evt_simphot_h.isValid() || evt_simphot_h->empty() ) {
169  // std::cerr << "[NeutrinoMCFlash] SimPhotons product is not valid or empty." << std::endl;
170  // e.put(std::move(opflashes));
171  // return;
172  // }
173 
174 
175  // if(evt_simphot_h->size() != geo->NOpDets()) {
176  // std::cout << "[NeutrinoMCFlash] Unexpected # of channels in simphotons!" << std::endl;
177  // e.put(std::move(opflashes));
178  // return;
179  // }
180 
181  // art::Handle<std::vector<sim::SimPhotonsLite> > evt_simphot_h;
182  // e.getByLabel(_simphot_label, _simphot_insta, evt_simphot_h);
183 
184  // if( !evt_simphot_h.isValid() || evt_simphot_h->empty() ) {
185  // std::cerr << "[NeutrinoMCFlash] SimPhotons product is not valid or empty." << std::endl;
186  // e.put(std::move(opflashes));
187  // return;
188  // }
189 
190  //std::vector<art::Handle<std::vector<sim::SimPhotonsLite> >> evt_simphot_hs;
191  //e.getManyByType(evt_simphot_hs);
192  auto evt_simphot_hs = e.getMany<std::vector<sim::SimPhotonsLite> >();
193  if( evt_simphot_hs.size() == 0 ) {
194  std::cerr << "[NeutrinoMCFlash] SimPhotonsLite product is not valid or empty." << std::endl;
195  e.put(std::move(opflashes));
196  return;
197  }
198 
199 
200  // opdet=>opchannel mapping
201  std::vector<size_t> opdet2opch(geo->NOpDets(),0);
202  for(size_t opch=0; opch<opdet2opch.size(); ++opch){
203  opdet2opch[geo->OpDetFromOpChannel(opch)] = opch;
204  }
205 
206  // Get the OpChannel of the PD to use
207  std::vector<int> opch_to_use = lightana::PDNamesToList(_pd_to_use);
208 
209  // auto const & evt_trigger = (*evt_trigger_h)[0];
210  // auto const trig_time = evt_trigger.TriggerTime();
211  auto const trig_time = clock_data.TriggerOffsetTPC();
212 
213  if (_debug) std::cout << "trig_time: " << trig_time << std::endl;
214  if (_debug) std::cout << "G4ToElecTime(0): " << clock_data.G4ToElecTime(0) << std::endl;
215  if (_debug) std::cout << "G4ToElecTime(1000): " << clock_data.G4ToElecTime(1000) << std::endl;
216  if (_debug) std::cout << "Number of OpDets: " << geo->NOpDets() << std::endl;
217 
218  double nuTime = -1.e9;
219  if (_debug) std::cout << "We have " << evt_mctruth_h->size() << " mctruth events." << std::endl;
220  for (size_t n = 0; n < evt_mctruth_h->size(); n++) {
221 
222  simb::MCTruth const& evt_mctruth = (*evt_mctruth_h)[n];
223  if (_debug) std::cout << "Origin: " << evt_mctruth.Origin() << std::endl;
224  if (evt_mctruth.Origin() != 1 ) continue;
225  if (_debug) std::cout << "We have " << evt_mctruth.NParticles() << " particles." << std::endl;
226  for (int p = 0; p < evt_mctruth.NParticles(); p++) {
227 
228  simb::MCParticle const& par = evt_mctruth.GetParticle(p);
229  //if (par.PdgCode() != 14) continue;
230  if (_debug){
231  std::cout << "Particle pdg: " << par.PdgCode() << std::endl;
232  std::cout << "new Particle time: " << par.T() << std::endl;
233  std::cout << "new converted: " << clock_data.G4ToElecTime(par.T()) - trig_time << std::endl;
234  std::cout << std::endl;
235  }
236  if (std::abs(par.PdgCode()) == 14 || std::abs(par.PdgCode()) == 12) {
237  nuTime = par.T();
238  }
239  }
240  }
241 
242  if (nuTime == -1.e9) {
243  std::cout << "[NeutrinoMCFlash] No neutrino found." << std::endl;
244  e.put(std::move(opflashes));
245  return;
246  }
247 
248  std::cout << "[NeutrinoMCFlash] Neutrino G4 interaction time: " << nuTime << std::endl;
249 
250  std::vector<std::vector<double> > pmt_v(1,std::vector<double>(geo->NOpDets(),0));
251  _pe_total = _pe_cherenkov = _pe_scintillation = 0;
252  _cherenkov_time_v.clear();
253  _cherenkov_pmt_v.clear();
254  _scintillation_time_v.clear();
255  _scintillation_pmt_v.clear();
256 
257  // float sampling = time_service->OpticalClock().Frequency() / 1000.0; // GHz
258  float start_window = clock_data.TriggerOffsetTPC(); // us
259  if (_debug) std::cout << "start_window " << start_window << std::endl;
260 
261 
262  // auto const& Photons = *(evt_simphot_h);
263  // for (sim::SimPhotonsLite const& photons: Photons) {
264  // unsigned int const nPhotons = std::accumulate(
265  // photons.DetectedPhotons.begin(), photons.DetectedPhotons.end(),
266  // 0U, [](auto sum, auto const& entry){ return sum + entry.second; }
267  // );
268  // std::cout << "Channel=" << photons.OpChannel << " has " << nPhotons << " photons (format: [tick] photons):" << std::endl;
269  // for (auto const& pair: photons.DetectedPhotons) {
270  // std::cout << "\t [" << pair.first << "] " << pair.second << std::endl;
271  // }
272  // }
273  std::cout << "We have " << evt_simphot_hs.size() << " SimPhoton collections" << std::endl;
274 
275  float nuTime_elec = clock_data.G4ToElecTime(nuTime) - trig_time;
276 
277  for (const art::Handle<std::vector<sim::SimPhotonsLite>> &evt_simphot_h: evt_simphot_hs) {
278 
279  bool reflected = (evt_simphot_h.provenance()->productInstanceName() == "Reflected");
280 
281  // float qe = _qe_direct;
282  // if (reflected) qe = _qe_refl;
283 
284  for(sim::SimPhotonsLite const& photons: *(evt_simphot_h)) {
285 
286  int opch = photons.OpChannel;
287 
288 
289  for(auto const& pair: photons.DetectedPhotons) {
290 
291  float photon_time = pair.first - start_window;
292 
293  float photon_time_elec = clock_data.G4ToElecTime(photon_time) - trig_time;
294 
295  if (_debug && !reflected) {
296  std::cout << "Photon time: " << photon_time_elec
297  << " - Neutrino time: " << nuTime_elec << std::endl;
298  }
299 
300 
301  if (photon_time_elec > nuTime_elec + _window_length) continue;
302  if (photon_time_elec > nuTime_elec - _pre_window){
303 
304  auto iter = std::find(opch_to_use.begin(), opch_to_use.end(), opch);
305  if (iter == opch_to_use.end()) continue;
306 
307  auto const& pt = geo->OpDetGeoFromOpChannel(opch).GetCenter();
308 
309  if(pt.X() < 0 && _tpc == 1) continue;
310  if(pt.X() > 0 && _tpc == 0) continue;
311 
312  pmt_v[0][opch] += pair.second;
313  _pe_total++;
314 
315  // for (int i = 0; i < pair.second; i++) {
316  // float r = _random.Uniform(1.);
317  // if(r < qe) {
318  // pmt_v[0][opch] += 1;
319  // _pe_total ++;
320  // }
321  // }
322  }
323  }
324  }
325  }
326 
327  double Ycenter, Zcenter, Ywidth, Zwidth;
328  GetFlashLocation(pmt_v[0], Ycenter, Zcenter, Ywidth, Zwidth);
329 
330  recob::OpFlash flash(nuTime_elec, // time w.r.t. trigger
331  0, // time width
332  nuTime_elec, // flash time in elec clock
333  0., // frame (?)
334  pmt_v[0], // pe per pmt
335  0, 0, 1, // this are just default values
336  Ycenter, Ywidth, Zcenter, Zwidth); // flash location
337 
338  std::cout << "[NeutrinoMCFlash] MC Flash Time: " << flash.Time() << std::endl;
339  std::cout << "[NeutrinoMCFlash] MC Flash PE: " << flash.TotalPE() << std::endl;
340  // for (size_t i = 0; i < pmt_v[0].size(); i++) {
341  // std::cout << "ch " << i << " => " << pmt_v[0][i] << std::endl;
342  // }
343 
344  opflashes->emplace_back(std::move(flash));
345 
346  e.put(std::move(opflashes));
347 
348  _tree1->Fill();
349 }
SBNDMCFlash::_scintillation_time_v
std::vector< double > _scintillation_time_v
Definition: SBNDMCFlash_module.cc:83
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
p
pdgs p
Definition: selectors.fcl:22
SBNDMCFlash::_scintillation_pmt_v
std::vector< int > _scintillation_pmt_v
Definition: SBNDMCFlash_module.cc:84
lightana::PDNamesToList
std::vector< int > PDNamesToList(std::vector< std::string > pd_names)
Definition: sbndcode/sbndcode/OpDetReco/OpFlash/FlashFinder/FlashFinderFMWKInterface.cxx:91
SBNDMCFlash::_mctruth_label
std::string _mctruth_label
Definition: SBNDMCFlash_module.cc:63
SBNDMCFlash::_pd_to_use
std::vector< std::string > _pd_to_use
Definition: SBNDMCFlash_module.cc:67
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
SBNDMCFlash::_cherenkov_pmt_v
std::vector< int > _cherenkov_pmt_v
Definition: SBNDMCFlash_module.cc:82
SBNDMCFlash::_cherenkov_time_v
std::vector< double > _cherenkov_time_v
Definition: SBNDMCFlash_module.cc:81
sim::SimPhotonsLite
Compact representation of photons on a channel.
Definition: SimPhotons.h:103
SBNDMCFlash::GetFlashLocation
void GetFlashLocation(std::vector< double >, double &, double &, double &, double &)
Definition: SBNDMCFlash_module.cc:351
e
do i e
Definition: updateVersionList.sh:68
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10

Member Data Documentation

std::vector<int> SBNDMCFlash::_cherenkov_pmt_v
private

Definition at line 82 of file SBNDMCFlash_module.cc.

std::vector<double> SBNDMCFlash::_cherenkov_time_v
private

Definition at line 81 of file SBNDMCFlash_module.cc.

bool SBNDMCFlash::_debug
private

Definition at line 71 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_event
private

Definition at line 77 of file SBNDMCFlash_module.cc.

std::string SBNDMCFlash::_mctruth_label
private

Definition at line 63 of file SBNDMCFlash_module.cc.

std::vector<std::string> SBNDMCFlash::_pd_to_use
private

Definition at line 67 of file SBNDMCFlash_module.cc.

opdet::sbndPDMapAlg SBNDMCFlash::_pds_map
private

Definition at line 73 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_pe_cherenkov = 0
private

Definition at line 80 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_pe_scintillation = 0
private

Definition at line 79 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_pe_total = 0
private

Definition at line 78 of file SBNDMCFlash_module.cc.

float SBNDMCFlash::_pre_window
private

Definition at line 70 of file SBNDMCFlash_module.cc.

float SBNDMCFlash::_qe_direct
private

Definition at line 69 of file SBNDMCFlash_module.cc.

float SBNDMCFlash::_qe_refl
private

Definition at line 69 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_run
private

Definition at line 77 of file SBNDMCFlash_module.cc.

std::vector<int> SBNDMCFlash::_scintillation_pmt_v
private

Definition at line 84 of file SBNDMCFlash_module.cc.

std::vector<double> SBNDMCFlash::_scintillation_time_v
private

Definition at line 83 of file SBNDMCFlash_module.cc.

std::string SBNDMCFlash::_simphot_insta
private

Definition at line 66 of file SBNDMCFlash_module.cc.

std::string SBNDMCFlash::_simphot_label
private

Definition at line 65 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_subrun
private

Definition at line 77 of file SBNDMCFlash_module.cc.

int SBNDMCFlash::_tpc
private

Definition at line 68 of file SBNDMCFlash_module.cc.

TTree* SBNDMCFlash::_tree1
private

Definition at line 85 of file SBNDMCFlash_module.cc.

std::string SBNDMCFlash::_trigger_label
private

Definition at line 64 of file SBNDMCFlash_module.cc.

float SBNDMCFlash::_window_length
private

Definition at line 70 of file SBNDMCFlash_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5