Jpp 20.0.0-72-g597b30bc9
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JEffectiveMassCountingOffline.cc File Reference

Example program to histogram neutrino effective mass and trigger efficiency for triggered events by counting events inside the can volume. More...

#include <string>
#include <iostream>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JAAnet/JVolume.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSirene/JVisibleEnergyToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JROOT/JManager.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JAbstractHistogram.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to histogram neutrino effective mass and trigger efficiency for triggered events by counting events inside the can volume.


All generated events should therefore be available in the input file(s).
This implies that option JSirene -N 0 should be used and option JTriggerEfficiency -O should not be used.
Alternatively, the option -C can be used to specify a cylinder as the can.
The unit of the contents of the histogram is $Mton$ or $km^{3}$.

Author
mdejong, bstrandberg, bjung

Definition in file JEffectiveMassCountingOffline.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 84 of file JEffectiveMassCountingOffline.cc.

85{
86 using namespace std;
87 using namespace JPP;
88 typedef JAbstractHistogram<double> JAbstractHistogram_t;
89
90 JMultipleFileScanner_t generatorFiles;
91 JMultipleFileScanner_t triggerFiles;
92 string outputFile;
93 double margin;
94 double epsilon;
95 bool logx;
96 int numberOfBins;
97 JCylinder3D cylinder;
98 std::string option;
99 JAbstractHistogram_t X;
100 int debug;
101
102
103 try {
104
105 JParser<> zap("Example program to histogram neutrino effective mass and trigger effifiency for triggered events with offline files. Visible energy will be computed inside the user cylinder.");
106
107 zap['f'] = make_field(generatorFiles);
108 zap['F'] = make_field(triggerFiles);
109 zap['o'] = make_field(outputFile)
110 = "Meff.root";
111 zap['m'] = make_field(margin, "Margin around the volume in which the considered events must reside")
112 = 0.0;
113 zap['e'] = make_field(epsilon, "Minimal amount of visible energy assigned to events") = 1e-8;
114 zap['X'] = make_field(logx, "Use logarithm of energy");
115 zap['x'] = make_field(X, "Histogram binning")
116 = JAbstractHistogram_t(100, 1.0, 100.0);
117 zap['C'] = make_field(cylinder, "User cylinder");
118 zap['O'] = make_field(option, "Result option")
119 = Mass_t, Volume_t;
120 zap['d'] = make_field(debug)
121 = 2;
122
123 zap(argc, argv);
124 }
125 catch(const exception &error) {
126 FATAL(error.what() << endl);
127 }
128
129 try{
130
131 //------------------------------------------------------------------------------
132 // Histogram definitions
133 //------------------------------------------------------------------------------
134
135 setLongprint(cout);
136
137 double Xmin = numeric_limits<double>::max();
138 double Xmax = numeric_limits<double>::lowest();
139
140 JEvtWeightFileScannerSet<> scanners1(generatorFiles);
141
142 numberOfBins = X.getNumberOfBins();
143 Xmin = X.getLowerLimit();
144 Xmax = X.getUpperLimit();
145
146 JManager<int, TH1D> hM0 (new TH1D("hM0[%]", NULL, numberOfBins, Xmin, Xmax));
147 JManager<int, TH1D> hNM0(new TH1D("hNM0[%]", NULL, numberOfBins, Xmin, Xmax));
148 JManager<int, TH1D> hnT0(new TH1D("hnT0[%]", NULL, numberOfBins, Xmin, Xmax));
149 JManager<int, TH1D> hT0 (new TH1D("hT0[%]", NULL, numberOfBins, Xmin, Xmax));
150
151 JManager<int, TH1D> hM1 (new TH1D("hM1[%]", NULL, numberOfBins, Xmin, Xmax));
152 JManager<int, TH1D> hNM1(new TH1D("hNM1[%]", NULL, numberOfBins, Xmin, Xmax));
153 JManager<int, TH1D> hnT1(new TH1D("hnT1[%]", NULL, numberOfBins, Xmin, Xmax));
154 JManager<int, TH1D> hT1 (new TH1D("hT1[%]", NULL, numberOfBins, Xmin, Xmax));
155
156 JManager<int, TH1D> hM2 (new TH1D("hM2[%]", NULL, numberOfBins, Xmin, Xmax));
157 JManager<int, TH1D> hNM2(new TH1D("hNM2[%]", NULL, numberOfBins, Xmin, Xmax));
158 JManager<int, TH1D> hnT2(new TH1D("hnT2[%]", NULL, numberOfBins, Xmin, Xmax));
159 JManager<int, TH1D> hT2 (new TH1D("hT2[%]", NULL, numberOfBins, Xmin, Xmax));
160
161 JManager<int, TH1D> hM3 (new TH1D("hM3[%]", NULL, numberOfBins, Xmin, Xmax));
162 JManager<int, TH1D> hNM3(new TH1D("hNM3[%]", NULL, numberOfBins, Xmin, Xmax));
163 JManager<int, TH1D> hnT3(new TH1D("hnT3[%]", NULL, numberOfBins, Xmin, Xmax));
164 JManager<int, TH1D> hT3 (new TH1D("hT3[%]", NULL, numberOfBins, Xmin, Xmax));
165
166 JCylinder3D can = cylinder;
167 can.addMargin(margin);
168
169 NOTICE("Cylinder, including margin: " << can << endl);
170 //------------------------------------------------------------------------------
171 // Loop over generator files
172 //------------------------------------------------------------------------------
173
174 for (JEvtWeightFileScannerSet<>::iterator scanner = scanners1.begin(); scanner != scanners1.end(); ++scanner) {
175
176 const JHead header = scanner->getHeader();
177
178 NOTICE("Scanning file type " << scanner->getName() << ": " << endl << header << endl);
179
180
181
182 for (JMultipleFileScanner<Evt> in(scanner->getFilelist()); in.hasNext(); ) {
183
184 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
185
186 Evt* event = in.next();
187
188 if (!has_primary(*event)) {
189 WARNING(endl << "Could not find primary for the following event:" << endl << *event << endl);
190 continue;
191 }
192
193 for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
194 track->clearusr(); // Prevent visible energy from being calculated based on `mc_usr_keys::energy_lost_in_can`
195 }
196
197 const Trk& primary = get_primary(*event);
198 const JVertex3D& vertex = getVertex (primary);
199
200
201
202 double Evis = getVisibleEnergy (*event, can);
203 double EvisLL = (has_leading_lepton(*event) ?
204 getVisibleEnergyLeadingLepton(*event, can) :
205 0.0);
206 double EvisHad = Evis - EvisLL;
207
208 if (Evis == 0) Evis = epsilon;
209 if (EvisLL == 0) EvisLL = epsilon;
210 if (EvisHad == 0) EvisHad = epsilon;
211
212 const double x0 = logx ? log10(primary.E) : primary.E;
213 const double x1 = logx ? log10(Evis) : Evis;
214 const double x2 = logx ? log10(EvisLL) : EvisLL;
215 const double x3 = logx ? log10(EvisHad) : EvisHad;
216
217 const int interactionID0 = getInteractionID(*event, false, false);
218 const int interactionID1 = getInteractionID(*event, true, false);
219
220 hnT0[interactionID0]->Fill(x0);
221 hnT1[interactionID0]->Fill(x1);
222 hnT2[interactionID0]->Fill(x2);
223 hnT3[interactionID0]->Fill(x3);
224
225 hnT0[interactionID1]->Fill(x0);
226 hnT1[interactionID1]->Fill(x1);
227 hnT2[interactionID1]->Fill(x2);
228 hnT3[interactionID1]->Fill(x3);
229
230 if (can.is_inside(vertex)) {
231 hNM0[interactionID0]->Fill(x0);
232 hNM1[interactionID0]->Fill(x1);
233 hNM2[interactionID0]->Fill(x2);
234 hNM3[interactionID0]->Fill(x3);
235
236 hNM0[interactionID1]->Fill(x0);
237 hNM1[interactionID1]->Fill(x1);
238 hNM2[interactionID1]->Fill(x2);
239 hNM3[interactionID1]->Fill(x3);
240 }
241 if (abs(primary.type) == TRACK_TYPE_NUTAU) {
242
243 const int interactionID2 = getInteractionID(*event, false, true);
244 const int interactionID3 = getInteractionID(*event, true, true);
245
246 hnT0[interactionID2]->Fill(x0);
247 hnT1[interactionID2]->Fill(x1);
248 hnT2[interactionID2]->Fill(x2);
249 hnT3[interactionID2]->Fill(x3);
250
251 hnT0[interactionID3]->Fill(x0);
252 hnT1[interactionID3]->Fill(x1);
253 hnT2[interactionID3]->Fill(x2);
254 hnT3[interactionID3]->Fill(x3);
255
256 if(can.is_inside(vertex)){
257 hNM0[interactionID2]->Fill(x0);
258 hNM1[interactionID2]->Fill(x1);
259 hNM2[interactionID2]->Fill(x2);
260 hNM3[interactionID2]->Fill(x3);
261
262 hNM0[interactionID3]->Fill(x0);
263 hNM1[interactionID3]->Fill(x1);
264 hNM2[interactionID3]->Fill(x2);
265 hNM3[interactionID3]->Fill(x3);
266 }
267 }
268 }
269 STATUS(endl);
270 }
271
272
273 //------------------------------------------------------------------------------
274 // Loop over trigger files
275 //------------------------------------------------------------------------------
276 JEvtWeightFileScannerSet<> scanners2(triggerFiles);
277
278 for (JEvtWeightFileScannerSet<>::iterator scanner = scanners2.begin(); scanner != scanners2.end(); ++scanner) {
279
280 const JHead header = scanner->getHeader();
281
282 NOTICE("Scanning file type " << scanner->getName() << ": " << endl << header << endl);
283
284 for (JMultipleFileScanner<Evt> in(scanner->getFilelist()); in.hasNext(); ) {
285
286 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
287
288 Evt* event = in.next();
289
290 if (!has_primary(*event)) {
291 WARNING(endl << "Could not find primary for the following event:" << endl << *event << endl);
292 continue;
293 }
294
295 for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
296 track->clearusr(); // Prevent visible energy from being calculated based on `mc_usr_keys::energy_lost_in_can`
297 }
298
299 const Trk& primary = get_primary(*event);
300
301 double Evis = getVisibleEnergy (*event, can);
302 double EvisLL = (has_leading_lepton(*event) ?
303 getVisibleEnergyLeadingLepton(*event, can) :
304 0.0);
305 double EvisHad = Evis - EvisLL;
306
307 if (Evis == 0) Evis = epsilon;
308 if (EvisLL == 0) EvisLL = epsilon;
309 if (EvisHad == 0) EvisHad = epsilon;
310
311 const double x0 = logx ? log10(primary.E) : primary.E;
312 const double x1 = logx ? log10(Evis) : Evis;
313 const double x2 = logx ? log10(EvisLL) : EvisLL;
314 const double x3 = logx ? log10(EvisHad) : EvisHad;
315
316 const int interactionID0 = getInteractionID(*event, false, false);
317 const int interactionID1 = getInteractionID(*event, true, false);
318
319 hM0[interactionID0]->Fill(x0);
320 hM1[interactionID0]->Fill(x1);
321 hM2[interactionID0]->Fill(x2);
322 hM3[interactionID0]->Fill(x3);
323
324 hM0[interactionID1]->Fill(x0);
325 hM1[interactionID1]->Fill(x1);
326 hM2[interactionID1]->Fill(x2);
327 hM3[interactionID1]->Fill(x3);
328
329 if (abs(primary.type) == TRACK_TYPE_NUTAU) {
330
331 const int interactionID2 = getInteractionID(*event, false, true);
332 const int interactionID3 = getInteractionID(*event, true, true);
333
334 hM0[interactionID2]->Fill(x0);
335 hM1[interactionID2]->Fill(x1);
336 hM2[interactionID2]->Fill(x2);
337 hM3[interactionID2]->Fill(x3);
338
339 hM0[interactionID3]->Fill(x0);
340 hM1[interactionID3]->Fill(x1);
341 hM2[interactionID3]->Fill(x2);
342 hM3[interactionID3]->Fill(x3);
343 }
344 }
345 STATUS(endl);
346 }
347
348
349 //------------------------------------------------------------------------------
350 // convert 'generated' and 'triggered' histograms to effective mass and trigger efficiency
351 //------------------------------------------------------------------------------
352 double scale = 0.0;
353
354 if (option == Mass_t) {
355 scale = can.getVolume() * DENSITY_SEA_WATER * 1e-6; // Mton
356 } else if (option == Volume_t) {
357 scale = can.getVolume() * 1.0e-9; // km^3
358 }
359
360 for (int type : get_keys(hM0)) {
361 //Trigger Efficiency
362 hT0.insert(make_pair(type, (TH1D*) divide(*hM0[type], *hnT0[type], MAKE_STRING("hT0[" << type << "]"))));
363 hT1.insert(make_pair(type, (TH1D*) divide(*hM1[type], *hnT1[type], MAKE_STRING("hT1[" << type << "]"))));
364 hT2.insert(make_pair(type, (TH1D*) divide(*hM2[type], *hnT2[type], MAKE_STRING("hT2[" << type << "]"))));
365 hT3.insert(make_pair(type, (TH1D*) divide(*hM3[type], *hnT3[type], MAKE_STRING("hT3[" << type << "]"))));
366
367 //EffMass
368 hM0[type]->Scale(scale);
369 hM1[type]->Scale(scale);
370 hM2[type]->Scale(scale);
371 hM3[type]->Scale(scale);
372
373 hM0[type]->Divide(hNM0[type]);
374 hM1[type]->Divide(hNM1[type]);
375 hM2[type]->Divide(hNM2[type]);
376 hM3[type]->Divide(hNM3[type]);
377
378 }
379
380 TFile out(outputFile.c_str(), "recreate");
381
382 out << hM0 << hM1 << hM2 << hM3 << hT0 << hT1 << hT2 << hT3;
383
384 out.Close();
385 }
386 catch(const JException& error) {
387 FATAL(error << endl);
388 }
389}
scale
void scale(vector< double > &v, double c)
scale vector content
Definition JBeaconSimulator.cc:73
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
setLongprint
void setLongprint(std::ostream &out)
Set long print option.
Definition JManip.hh:132
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
WARNING
#define WARNING(A)
Definition JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_STRING
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
numberOfBins
int numberOfBins
number of bins for average CDF integral of optical module
Definition JSirene.cc:73
JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236
JAANET::JHead::getHeader
const JHead & getHeader() const
Get header.
Definition JHead.hh:1270
JLANG::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306
JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition JMultipleFileScanner.hh:409
JAANET::has_primary
bool has_primary(const Evt &evt)
Auxiliary function to check if an event contains a primary track.
Definition JAAnetToolkit.hh:647
JAANET::getVertex
JVertex3D getVertex(const Trk &track)
Get vertex.
Definition JAAnetToolkit.hh:257
JAANET::has_leading_lepton
bool has_leading_lepton(const Evt &event)
Auxiliary function to check if an event contains the leading lepton of a neutrino interaction.
Definition JAAnetToolkit.hh:689
JAANET::get_primary
const Trk & get_primary(const Evt &evt)
Auxiliary function to retrieve the primary track of an event.
Definition JAAnetToolkit.hh:668
JAANET::getInteractionID
int getInteractionID(const Evt &event, const bool useScatteringType=true, const bool useDecayType=true)
Retrieve interaction identifier.
Definition JAAnetToolkit.hh:1260
JLANG::get_keys
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
JPHYSICS::DENSITY_SEA_WATER
static const double DENSITY_SEA_WATER
Fixed environment values.
Definition JPhysics/JConstants.hh:44
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:48
JSIRENE::getVisibleEnergyLeadingLepton
double getVisibleEnergyLeadingLepton(const Trk &, const JCylinder3D &)
JSIRENE::getVisibleEnergy
double getVisibleEnergy(const Trk &, const JCylinder3D &)
Get the visible energy of a track.
Definition JVisibleEnergyToolkit.hh:73
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
JAANET::can
The cylinder used for photon tracking.
Definition JHead.hh:575
JAANET::primary
Primary particle.
Definition JHead.hh:1174
JAANET::primary::type
int type
Particle type.
Definition JHead.hh:1204
JSUPPORT::JEvtWeightFileScannerSet
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
Definition JEvtWeightFileScannerSet.hh:58
JSUPPORT::JMultipleFileScanner_t
Auxiliary base class for list of file names.
Definition JMultipleFileScanner.hh:51
JTOOLS::JAbstractHistogram
Simple data structure for histogram binning.
Definition JAbstractHistogram.hh:26
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
Generated by doxygen 1.12.0