Jpp  18.3.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEffectiveMass1D.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 
7 
8 #include "JDAQ/JDAQEventIO.hh"
9 
12 
13 #include "JAAnet/JAAnetToolkit.hh"
14 #include "JAAnet/JHeadToolkit.hh"
15 #include "JAAnet/JVolume.hh"
16 
20 #include "JSupport/JSupport.hh"
21 
22 #include "JDetector/JDetector.hh"
24 
25 #include "JPhysics/JConstants.hh"
26 #include "JPhysics/JGeane.hh"
27 
28 #include "JReconstruction/JEvt.hh"
30 
31 #include "JROOT/JManager.hh"
32 
33 #include "Jeep/JParser.hh"
34 #include "Jeep/JMessage.hh"
35 
36 #include "TROOT.h"
37 #include "TFile.h"
38 #include "TH1D.h"
39 
40 
41 namespace {
42 
43  /**
44  * Check if given event has a muon intersecting the given cylinder.
45  *
46  * \param event event
47  * \param cylinder cylinder
48  * \return true if event has intersecting muon; else false.
49  */
50  bool hasIntersectingMuon(const Evt& event,
51  const JGEOMETRY3D::JCylinder3D& cylinder)
52  {
53  using namespace std;
54  using namespace JPP;
55 
56  for (vector<Trk>::const_iterator i = event.mc_trks.cbegin(); i != event.mc_trks.cend(); ++i) {
57 
58  if (is_finalstate(*i) && is_muon(*i)) {
59 
60  const JCylinder3D::intersection_type& intersection = cylinder.getIntersection(getAxis(*i));
61 
62  const double Lmuon = gWater.getX(i->E, MASS_MUON / getSinThetaC());
63  const double Leff = (intersection.first < 0.0 ?
64  min(Lmuon, intersection.second) :
65  min(Lmuon, intersection.second) - intersection.first);
66 
67  if (Leff > 0.0) { return true; }
68  }
69  }
70 
71  return false;
72  }
73 }
74 
75 
76 /**
77  * \file
78  *
79  * Example program to histogram neutrino effective mass for triggered events inside the can volume.\n
80  * For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.\n
81  * The unit of the contents of the histogram is Mton.
82  * \author mdejong, bstrandberg
83  */
84 int main(int argc, char **argv)
85 {
86  using namespace std;
87  using namespace JPP;
88  using namespace KM3NETDAQ;
89 
90 
91  //------------------------------------------------------------------------------
92  // parse command line arguments, check that header can be found in the input file
93  //------------------------------------------------------------------------------
94 
95  JMultipleFileScanner_t inputFiles;
96  string outputFile;
97  double wall;
98  bool logx;
99  int numberOfBins;
100  string detectorFile;
101  int debug;
102 
103 
104  try {
105 
106  JParser<> zap("Example program to histogram neutrino effective mass for triggered events.");
107 
108  zap['f'] = make_field(inputFiles);
109  zap['o'] = make_field(outputFile) = "Meff.root";
110  zap['R'] = make_field(wall, "Addition margin around the volume in which the considered events must reside") = 0.0;
111  zap['X'] = make_field(logx, "Use logarithm of energy");
112  zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10;
113  zap['d'] = make_field(debug) = 1;
114  zap['a'] = make_field(detectorFile, "Detector file: if not provided, trigger fraction is not calculated") = "";
115 
116  zap(argc, argv);
117  }
118  catch(const exception &error) {
119  FATAL(error.what() << endl);
120  }
121 
122  //-------------------------------------------------------------------------
123  // Initialise the detector
124  //-------------------------------------------------------------------------
125 
127 
128  if (detectorFile != "") {
129 
130  try {
131  load(detectorFile, detector);
132  }
133  catch(const JException& error) {
134  FATAL(error);
135  }
136  }
137 
138 
139  //-------------------------------------------------------------------------
140  // Retrieve auxiliary header info
141  //-------------------------------------------------------------------------
142 
143  JCylinder3D canvol;
144  JCylinder3D instvol;
145 
146  vector<JVolume> volumes;
147 
148  double Xmin = numeric_limits<double>::max();
149  double Xmax = numeric_limits<double>::lowest();
150 
151  JEvtWeightFileScannerSet<> scanners(inputFiles);
152 
153  for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.cbegin(); scanner != scanners.cend(); ++scanner) {
154 
155  const JHead& header = scanner->getHeader();
156 
157  const JVolume volume(header, logx);
158 
159  if (volume.getXmin() < Xmin) { Xmin = volume.getXmin(); }
160  if (volume.getXmax() > Xmax) { Xmax = volume.getXmax(); }
161 
162  volumes.push_back(volume);
163 
164  // the can volume comes from the generator header and the coordinate system is consistent by construction
165  JCylinder3D can = get<JCylinder3D>(header);
166  canvol.addMargin(wall);
167 
168  if (can.getVolume() > canvol.getVolume()) { canvol = can; }
169 
170  // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced
171  // to the one used in the event generator
172  JPosition3D center = getPosition(get<Vec>(header));
173  detector -= center;
174 
175  JCylinder3D inst( detector.begin(), detector.end() );
176  instvol.addMargin(wall);
177 
178  if (inst.getVolume() > instvol.getVolume()) { instvol = inst; }
179 
180  detector += center;
181  }
182 
183 
184  //------------------------------------------------------------------------------
185  // initialise output
186  //------------------------------------------------------------------------------
187 
188  TFile out(outputFile.c_str(), "recreate");
189 
190  TH1::SetDefaultSumw2();
191 
192  JManager<int, TH1D> hM(new TH1D("hM[%]", "effective mass (can) [kg]" , numberOfBins, Xmin, Xmax));
193  JManager<int, TH1D> hm(new TH1D("hm[%]", "effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));
194  JManager<int, TH1D> hNM(new TH1D("hNM[%]", "normalization effective mass (can) [kg]" , numberOfBins, Xmin, Xmax));
195  JManager<int, TH1D> hNm(new TH1D("hNm[%]", "normalization effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));
196 
197 
198  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
199 
200  const size_t scannerIndex = distance(scanners.begin(), scanner);
201 
202  const JVolume& volume = volumes[scannerIndex];
203 
204  NOTICE("JEffectiveMass1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl );
205 
206  //------------------------------------------------------------------------------
207  // loop over generated events in the input files
208  //------------------------------------------------------------------------------
209 
210  NOTICE("Scanning generated events for file type " << scannerIndex << endl);
211 
212  while (scanner->hasNext()) {
213 
214  STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl);
215 
216  const Evt* event = scanner->next();
217 
218  const Trk& primary = get_primary(*event);
219 
220  const JVertex3D& vertex = getVertex(*event);
221 
222  const bool isValid1 = (canvol.is_inside(vertex) || hasIntersectingMuon(*event, canvol));
223  const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));
224 
225  // events are filled to hM if they are inside the can
226  hNM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? scanner->getWeight(*event) : 0.0);
227 
228  // events are filled to hm if they are inside the instrumented volume
229  hNm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? scanner->getWeight(*event) : 0.0);
230 
231  } // end loop over generated events
232  STATUS(endl);
233 
234  if (scanner->getCounter() == 0) {
235  FATAL("JEffectiveMass1D: generated events not stored in the input file, JTriggerEfficiency should be run without option -O");
236  }
237 
238  //------------------------------------------------------------------------------
239  // loop over triggered events in the input files
240  //------------------------------------------------------------------------------
241 
242  NOTICE("Scanning triggered events for file type " << scannerIndex << endl);
243 
244  JEvtWeightFileScanner<JTriggeredFileScanner<> > in(scanner->getFilelist());
245 
246  while (in.hasNext()) {
247 
248  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
249 
251 
252  const Evt* event = mp;
253 
254  const Trk& primary = get_primary(*event);
255 
256  const JVertex3D& vertex = getVertex(*event);
257 
258  const bool isValid1 = (canvol.is_inside(vertex) || hasIntersectingMuon(*event, canvol));
259  const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));
260 
261  hM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? in.getWeight(*event) : 0.0);
262  hm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? in.getWeight(*event) : 0.0);
263 
264  } // end loop over generated events
265  STATUS(endl);
266  }
267 
268 
269  //------------------------------------------------------------------------------
270  // convert 'generated' and 'triggered' histograms to effective mass
271  //------------------------------------------------------------------------------
272 
273  vector<int> pdgIDs = get_keys(hM);
274 
275  for (vector<int>::const_iterator i = pdgIDs.begin(); i != pdgIDs.end(); ++i) {
276 
277  const int pdgID = *i;
278 
279  hM[pdgID]->Divide(hNM[pdgID]);
280  hM[pdgID]->Scale(canvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
281  hm[pdgID]->Divide(hNm[pdgID]);
282  hm[pdgID]->Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
283  }
284 
285  out << hM << hm;
286 
287  out.Close();
288 }
JVertex3D getVertex(const Trk &track)
Get vertex.
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
Double_t getXmin() const
Get minimal abscissa value.
Definition: JVolume.hh:91
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Energy loss of muon.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
static const double MASS_MUON
muon mass [GeV]
Double_t getXmax() const
Get maximal abscissa value.
Definition: JVolume.hh:102
const Trk & get_primary(const Evt &evt)
Get primary.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
virtual double getX(const double E0, const double E1) const override
Get distance traveled by muon.
Definition: JGeane.hh:349
static const double DENSITY_SEA_WATER
Fixed environment values.
Dynamic ROOT object management.
intersection_type getIntersection(const JAxis3D &axis) const
Get intersection points of axis with cylinder.
Definition: JCylinder3D.hh:277
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
bool is_finalstate(const Trk &track)
Test whether given track corresponds to a final state particle.
Auxiliary class for histogramming of effective volume.
Definition: JVolume.hh:29
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:226
Monte Carlo run header.
Definition: JHead.hh:1234
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
JAxis3D getAxis(const Trk &track)
Get axis.
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
Physics constants.
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:172
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
int type
MC: particle type in PDG encoding.
Definition: Trk.hh:24
double getVolume() const
Get volume.
Definition: JCylinder3D.hh:185
Auxiliary base class for list of file names.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
std::vector< filescanner_type >::iterator iterator
Utility class to parse command line options.
Primary particle.
Definition: JHead.hh:1174
const JHead & getHeader() const
Get header.
Definition: JHead.hh:1270
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
Template event-weighter-associated file scanner.
do set_variable DETECTOR_TXT $WORKDIR detector
General purpose class for multiple pointers.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
std::vector< filescanner_type >::const_iterator const_iterator
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:69
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:49
const array_type< JKey_t > & get_keys(const std::map< JKey_t, JValue_t, JComparator_t, JAllocator_t > &data)
Method to create array of keys of map.
Definition: JVectorize.hh:139
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
The cylinder used for photon tracking.
Definition: JHead.hh:575