Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JDomino.cc File Reference

Example program to verify Monte Carlo data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#include <sstream>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JROOT/JManager.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JWeight.hh"
#include "JROOT/JRootToolkit.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JLang/JException.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JLangToolkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to verify Monte Carlo data.

Author
mdejong

Definition in file JDomino.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 183 of file JDomino.cc.

184 {
185  using namespace std;
186  using namespace JPP;
187 
188  //-------------------------------------------------------------------------
189  // command-line arguments
190  //-------------------------------------------------------------------------
191 
192  JMultipleFileScanner<Evt> inputFile;
193  JLimit_t& numberOfEvents = inputFile.getLimit();
194  string outputFile;
195  string detectorFile;
196  set<int> hit_types;
197  int debug;
198 
199  try {
200 
201  JParser<> zap("Example program to verify Monte Carlo data.");
202 
203  zap['f'] = make_field(inputFile);
204  zap['n'] = make_field(numberOfEvents) = JLimit::max();
205  zap['o'] = make_field(outputFile) = "domino.root";
206  zap['a'] = make_field(detectorFile) = "";
207  zap['T'] = make_field(hit_types) = JPARSER::initialised();
208  zap['d'] = make_field(debug) = 2;
209 
210  zap(argc, argv);
211  }
212  catch(const exception &error) {
213  FATAL(error.what() << endl);
214  }
215 
216  cout.tie(&cerr);
217 
218 
219  //-------------------------------------------------------------------------
220  // load detector file
221  //-------------------------------------------------------------------------
222 
224 
225  if (detectorFile != "") {
226 
227  try {
228  load(detectorFile, detector);
229  }
230  catch(const JException& error) {
231  FATAL(error);
232  }
233  }
234 
235  const JPMTRouter router(detector);
236 
237  const JHead header = getHeader(inputFile);
238 
239 
240  //-------------------------------------------------------------------------
241  // try to determine coordinate origin
242  //-------------------------------------------------------------------------
243 
244  JPosition3D center = get<JPosition3D>(header);
245 
246  NOTICE("Apply detector offset " << center << endl);
247 
248  detector -= center;
249 
250 
251  //-------------------------------------------------------------------------
252  // calculate instrumented volume
253  //-------------------------------------------------------------------------
254 
255  JCylinder3D inst_vol(detector.begin(), detector.end());
256 
257  NOTICE("JDomino: Instrumented volume dimensions (x, y, r, zmin, zmax): " << inst_vol << endl);
258 
259 
260  //-------------------------------------------------------------------------
261  // histogram binning
262  //-------------------------------------------------------------------------
263 
264  typedef map<int, JWeight> map_type;
265 
266 
267  double xmin = 2.0;
268  double xmax = 8.0;
269 
270  // energy range is taken from -in order of priority- neutrino or muon data
271 
272  if (header.is_valid(&JHead::cut_nu)) {
273  xmin = log10(header.cut_nu.Emin);
274  xmax = log10(header.cut_nu.Emax);
275  } else if (header.is_valid(&JHead::cut_in)) {
276  xmin = log10(header.cut_in.Emin);
277  xmax = log10(header.cut_in.Emax);
278  }
279  const int nx = (int) ((xmax - xmin) / 0.1);
280 
281  const Double_t T[] = { -50.0, -20.0, -10.0, -5.0, -2.0, 0.0, +2.0, +5.0, +10.0, +15.0, +20.0, +30.0, +40.0, +50.0,
282  +75.0, +100.0, +150.0, +200.0, +250.0, +300.0, +400.0, +500.0, 1000.0 }; // [ns]
283 
284  //-------------------------------------------------------------------------
285  // histograms for events
286  //-------------------------------------------------------------------------
287 
288  TH1D hits("hits", NULL, 100, 0.0, 8.0); // number of mc_hits per event
289  TH1D trks("trks", NULL, 10001, -5000.5, 5000.5); // number of track_in type's per event
290  TH1D job ("job" , NULL, 10001, -5000.5, 5000.5); // number of photo-electrons by track_in type
291 
292  TProfile hits_per_E_in ("hits_per_E_in", "average number of hits per E_{#nu} inside instrumented volume", nx, xmin, xmax);
293  TProfile hits_per_E_out("hits_per_E_out", "average number of hits per E_{#nu} outside instrumented volume", nx, xmin, xmax);
294 
295  typedef JManager<int, TH1D> JManager_t;
296 
297  JManager_t npe_per_type(new TH1D("npe[%]", NULL, 5000, 0.5, 5000.5), '%', ios::fmtflags(ios::showpos));
298 
299  TH1D npe_per_pmt ("pmt", NULL, 100001, -0.5, 100000.5); // number of photo-electrons per PMT
300 
301 
302  //-------------------------------------------------------------------------
303  // histograms for neutrino
304  //-------------------------------------------------------------------------
305 
306  TH2D nuExD("nuExD", NULL, nx, xmin, xmax, 100, 0.0, 1000.0); // Enu versus distance between vertex and PMT
307  TH2D* nuExD_tmp = (TH2D*) nuExD.Clone("nuExD.tmp"); // normalisation
308 
309  TH2D nuExc("nuExc", NULL, nx, xmin, xmax, 100, -1.0, +1.0); // Enu versus emission angle
310  TH2D* nuExc_tmp = (TH2D*) nuExc.Clone("nuExc.tmp"); // normalisation
311 
312  TH2D nuDxc("nuDxc", NULL, 50, 0.0, 1000.0, 100, -1.0, +1.0); // distance between vertex and PMT versus emission angle
313  TH2D* nuDxc_tmp = (TH2D*) nuDxc.Clone("nuDxc.tmp"); // normalisation
314 
315  TH2D nuDxU("nuDxU", NULL, 50, 0.0, 1000.0, 100, -1.0, +1.0); // distance between vertex and PMT versus incidence angle
316  TH2D* nuDxU_tmp = (TH2D*) nuDxU.Clone("nuDxU.tmp"); // normalisation
317 
318  TH2D nuDxT("nuDxT", NULL, 50, 0.0, 1000.0, getSize(T) - 1, T); // distance between vertex and PMT versus time residual
319  TH2D* nuDxT_tmp = (TH2D*) nuDxT.Clone("nuDxT.tmp"); // normalisation
320 
321  nuExD.Sumw2();
322  nuExc.Sumw2();
323  nuDxc.Sumw2();
324  nuDxU.Sumw2();
325  nuDxT.Sumw2();
326 
327 
328  //-------------------------------------------------------------------------
329  // histograms for muon
330  //-------------------------------------------------------------------------
331 
332  TH2D muExR("muExR", NULL, nx, xmin, xmax, 30, 0.0, 300.0); // Emu versus distance between muon and PMT
333  TH2D* muExR_tmp = (TH2D*) muExR.Clone("muExR.tmp"); // normalisation
334 
335  TH2D muRxU("muRxU", NULL, 15, 0.0, 300.0, 100, -1.0, +1.0); // distance between muon and PMT versus incidence angle
336  TH2D* muRxU_tmp = (TH2D*) muRxU.Clone("muRxU.tmp"); // normalisation
337 
338  TH2D muRxT("muRxT", NULL, 15, 0.0, 300.0, getSize(T) - 1, T); // distance between muon and PMT versus time residual
339  TH2D* muRxT_tmp = (TH2D*) muRxT.Clone("muRxT.tmp"); // normalisation
340 
341  muExR.Sumw2();
342  muRxU.Sumw2();
343  muRxT.Sumw2();
344 
345 
346  //-------------------------------------------------------------------------
347  // loop over events
348  //-------------------------------------------------------------------------
349 
350  while (inputFile.hasNext()) {
351 
352  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
353 
354  const Evt* event = inputFile.next();
355 
356  double NPE = 0.0;
357 
358  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
359  NPE += getNPE(*hit);
360  }
361 
362  hits.Fill(log10((Double_t) event->mc_hits.size()));
363 
364  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
365  trks.Fill(track->type);
366  }
367 
368  JTrack3E neutrino;
369 
370  if (has_neutrino(*event)) {
371  neutrino = getTrack(get_neutrino(*event));
372  }
373 
374  const JVertex3D vertex = neutrino.getVertex();
375 
376  map_type npe_pmt; // temporary buffer to accumulate total npe's by individual PMTs
377  map_type npe_type; // temporary buffer to accumulate total npe's by hit origin track
378 
379 
380  //-------------------------------------------------------------------------
381  // loop over hits
382  //-------------------------------------------------------------------------
383 
384  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
385 
386  const int type = hit->type;
387  const double t1 = getTime(*hit);
388  const double npe = getNPE (*hit);
389 
390  npe_pmt[hit->pmt_id].put(npe);
391 
392  npe_type[0] .put(npe);
393  npe_type[type].put(npe);
394 
395  if (hit_types.empty() || hit_types.count(type) != 0) {
396 
397  vector<Trk>::const_iterator track = find_if(event->mc_trks.begin(),
398  event->mc_trks.end(),
399  make_predicate(&Trk::id, hit->origin));
400 
401  if (track == event->mc_trks.end()) {
402  ERROR("Hit " << *hit << " has no origin." << endl);
403  continue;
404  }
405 
406  if (count_if(event->mc_trks.begin(),
407  event->mc_trks.end(),
408  make_predicate(&Trk::id, hit->origin)) > 1) {
409  ERROR("Hit " << *hit << " has ambiguous origin." << endl);
410  continue;
411  }
412 
413  job.Fill((double) track->type, npe);
414 
415  if (router.hasPMT(hit->pmt_id)) {
416 
417  const JPMT& pmt = router.getPMT(hit->pmt_id);
418 
419  if (is_muon(*track)) {
420 
421  const JTrack3E muon = getTrack(*track);
422 
423  const double E = muon.getE();
424  const double x = log10(E);
425  const double R = muon.getDistance(pmt);
426  const double t0 = muon.getT (pmt);
427  const double ct1 = muon.getDot (pmt);
428 
429  muExR.Fill(x, R, getMuonWeight(E, npe));
430  muRxU.Fill(R, ct1, getMuonWeight(E, R, npe));
431  muRxT.Fill(R, t1 - t0, getMuonWeight(E, R, npe));
432 
433  } else if (has_neutrino(*event)) {
434 
435  const double E = neutrino.getE();
436  const double x = log10(E);
437  const double D = vertex.getDistance(pmt);
438  const double t0 = vertex.getT(pmt);
439  const double ct0 = neutrino.getDot(pmt.getPosition());
440  const double ct1 = vertex.getDot(pmt);
441 
442  nuExD.Fill(x, D, getNeutrinoWeight(E, npe));
443  nuExc.Fill(x, ct0, getNeutrinoWeight(E, npe));
444  nuDxc.Fill(D, ct0, getNeutrinoWeight(E, D, npe));
445  nuDxU.Fill(D, ct1, getNeutrinoWeight(E, D, npe));
446  nuDxT.Fill(D, t1 - t0, getNeutrinoWeight(E, D, npe));
447  }
448  }
449  } // end if PMT of hit found in PMT router
450  } // end loop over hits
451 
452 
453  //-------------------------------------------------------------------------
454  // continue with logic for event, fill histograms common to all input files
455  //-------------------------------------------------------------------------
456 
457  for (map_type::const_iterator i = npe_pmt.begin(); i != npe_pmt.end(); ++i) {
458  npe_per_pmt.Fill(i->second.getTotal()); // npe / PMT
459  }
460 
461  for (map_type::const_iterator i = npe_type.begin(); i != npe_type.end(); ++i) {
462  npe_per_type[i->first]->Fill(i->second.getTotal()); // npe / type
463  }
464 
465 
466  //-------------------------------------------------------------------------
467  // normalisation of muon histograms
468  //-------------------------------------------------------------------------
469 
470  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
471 
472  if (is_muon(*track)) {
473 
474  const JTrack3E muon = getTrack(*track);
475 
476  const double E = muon.getE();
477  const double x = log10(E);
478 
479  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
480 
481  const double R = muon.getDistance(*module);
482 
483  muExR_tmp->Fill(x, R, module->size());
484 
485  if (R < muRxU.GetXaxis()->GetXmax()) {
486  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
487  muRxU_tmp->Fill(R, muon.getDot(*pmt));
488  }
489  }
490 
491  if (R < muRxT.GetXaxis()->GetXmax()) {
492  for (Int_t iy = 1; iy <= muRxT_tmp->GetYaxis()->GetNbins(); ++iy) {
493  muRxT_tmp->Fill(R, muRxT_tmp->GetYaxis()->GetBinCenter(iy), muRxT_tmp->GetYaxis()->GetBinWidth(iy));
494  }
495  }
496  }
497  }
498  }
499 
500 
501  //-------------------------------------------------------------------------
502  // normalisation of neutrino histograms
503  //-------------------------------------------------------------------------
504 
505  if (has_neutrino(*event)) {
506 
507  const double E = neutrino.getE();
508  const double x = log10(E);
509 
510  if (inst_vol.is_inside(vertex))
511  hits_per_E_in .Fill(x, NPE);
512  else
513  hits_per_E_out.Fill(x, NPE);
514 
515  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
516 
517  const double D = vertex.getDistance(*module);
518  const double ct0 = neutrino.getDot(module->getPosition());
519 
520  nuExD_tmp->Fill(x, D, module->size());
521  nuExc_tmp->Fill(x, ct0, module->size());
522  nuDxc_tmp->Fill(D, ct0, module->size());
523 
524  if (D < nuDxU.GetXaxis()->GetXmax()) {
525  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
526  nuDxU_tmp->Fill(D, neutrino.getDot(*pmt));
527  }
528  }
529 
530  if (D < nuDxT.GetXaxis()->GetXmax()) {
531  for (Int_t iy = 1; iy <= nuDxT_tmp->GetYaxis()->GetNbins(); ++iy) {
532  nuDxT_tmp->Fill(D, nuDxT_tmp->GetYaxis()->GetBinCenter(iy), nuDxT_tmp->GetYaxis()->GetBinWidth(iy));
533  }
534  }
535  }
536  }
537  } // end loop over events
538 
539  STATUS(endl);
540 
541 
542  //-------------------------------------------------------------------------
543  // normalisation of histograms
544  //-------------------------------------------------------------------------
545 
546  TH1D *job_sorted = makeSortedH1(&job, "hits_by_pdg"); // hits by track_in particles (PDG codes) per event, sorted
547  TH1D *trks_sorted = makeSortedH1(&trks, "trks_sorted"); // track_in particles (PDG codes) per event, sorted
548 
549  if (inputFile.getCounter() != 0) {
550 
551  const Double_t W = 1.0 / ((Double_t) inputFile.getCounter());
552 
553  for (TH1D* p : { &hits, &npe_per_pmt, &job, &trks, job_sorted, trks_sorted }) {
554  p->Scale(W);
555  }
556 
557  for (JManager_t::iterator i = npe_per_type.begin(); i != npe_per_type.end(); ++i) {
558  i->second->Scale(W);
559  }
560 
561  nuExD.Divide(nuExD_tmp);
562  nuExc.Divide(nuExc_tmp);
563  nuDxc.Divide(nuDxc_tmp);
564  nuDxU.Divide(nuDxU_tmp);
565  nuDxT.Divide(nuDxT_tmp);
566 
567  muExR.Divide(muExR_tmp);
568  muRxU.Divide(muRxU_tmp);
569  muRxT.Divide(muRxT_tmp);
570  }
571 
572 
573  //-------------------------------------------------------------------------
574  // output
575  //-------------------------------------------------------------------------
576 
577  TFile out(outputFile.c_str(), "recreate");
578 
579  out << hits << job << *job_sorted << trks << *trks_sorted << hits_per_E_in << hits_per_E_out << npe_per_type << npe_per_pmt;
580 
581  out << nuExD << nuExc << nuDxc << nuDxU << nuDxT;
582  out << muExR << muRxU << muRxT;
583 
584  out.Write();
585  out.Close();
586 }
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:33
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
JTrack3E getTrack(const Trk &track)
Get track.
#define STATUS(A)
Definition: JMessage.hh:63
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Detector data structure.
Definition: JDetector.hh:80
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
double getDot(const JAxis3D &axis) const
Get cosine angle of impact of Cherenkov light on PMT.
Definition: JVertex3D.hh:175
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
size_t getSize(T(&array)[N])
Get size of c-array.
Definition: JLangToolkit.hh:32
3D track with energy.
Definition: JTrack3E.hh:30
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Cylinder object.
Definition: JCylinder3D.hh:39
Detector file.
Definition: JHead.hh:196
double getE() const
Get energy.
Definition: JTrack3E.hh:93
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:47
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
int debug
debug level
Definition: JSirene.cc:63
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
Monte Carlo run header.
Definition: JHead.hh:1091
#define FATAL(A)
Definition: JMessage.hh:67
int id
track identifier
Definition: Trk.hh:15
JVertex3D getVertex() const
Get vertex of this track.
Definition: JTrack3D.hh:99
double getDistance(const JVector3D &pos) const
Get distance.
Definition: JAxis3D.hh:213
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:145
double getDot(const JAxis3D &axis) const
Get cosine angle of impact of Cherenkov light on PMT.
Definition: JTrack3D.hh:170
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62