Jpp 20.0.0-72-g597b30bc9
the software that should make you happy
Loading...
Searching...
No Matches
JEffectiveMassCountingOffline.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3
4#include "TROOT.h"
5#include "TFile.h"
6#include "TH1D.h"
7
8#include "km3net-dataformat/offline/Head.hh"
9#include "km3net-dataformat/offline/Evt.hh"
10
11#include "JGeometry3D/JCylinder3D.hh"
12#include "JGeometry3D/JPosition3D.hh"
13
14#include "JAAnet/JVolume.hh"
15#include "JAAnet/JHeadToolkit.hh"
16#include "JAAnet/JParticleTypes.hh"
17#include "JAAnet/JAAnetToolkit.hh"
18
19#include "JSirene/JVisibleEnergyToolkit.hh"
20
21#include "JSupport/JSupport.hh"
22#include "JSupport/JMultipleFileScanner.hh"
23#include "JSupport/JEvtWeightFileScannerSet.hh"
24#include "JSupport/JMonteCarloFileSupportkit.hh"
25
26#include "JROOT/JManager.hh"
27
28#include "Jeep/JParser.hh"
29#include "Jeep/JMessage.hh"
30
31#include "JTools/JAbstractHistogram.hh"
32
33namespace {
34 const char* const Mass_t = "Mass";
35 const char* const Volume_t = "Volume";
36
37 /**
38 * Divide correlated histograms.
39 *
40 * \param h1 first histogram (numerator)
41 * \param h0 second histogram (denominator)
42 * \param name divided histogram name
43 * \return divided histogram
44 */
45 TH1* divide(const TH1& h1, const TH1& h0, const std::string name = "")
46 {
47 using namespace std;
48 using namespace JPP;
49
50 TH1* h2 = (TH1*) h1.Clone(name.c_str());
51
52 for (Int_t i = 0; i <= h1.GetNbinsX()+1; ++i) {
53
54 const Double_t y1 = h1.GetBinContent(i);
55 const Double_t y0 = h0.GetBinContent(i);
56
57 if (y0 != 0.0) {
58
59 const Double_t y3 = y1 / y0;
60 const Double_t w3 = sqrt(y1 * (y0 - y1) / (y0*y0*y0));
61
62 h2->SetBinContent(i, y3);
63 h2->SetBinError (i, w3);
64
65 } else {
66 ERROR("Bin " << h0.GetName() << "[" << i << "] empty." << endl);
67 }
68 }
69
70 return h2;
71 }
72}
73
74/**
75 * \file
76 *
77 * Example program to histogram neutrino effective mass and trigger efficiency for triggered events by counting events inside the can volume.\n
78 * All generated events should therefore be available in the input file(s).\n
79 * This implies that option <tt>JSirene -N 0</tt> should be used and option <tt>JTriggerEfficiency -O</tt> should not be used.\n
80 * Alternatively, the option <tt>-C</tt> can be used to specify a cylinder as the can.\n
81 * The unit of the contents of the histogram is \f$Mton\f$ or \f$km^{3}\f$.
82 * \author mdejong, bstrandberg, bjung
83 */
84int main(int argc, char **argv)
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}
JAAnetToolkit.hh
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
scale
void scale(vector< double > &v, double c)
scale vector content
Definition JBeaconSimulator.cc:73
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
main
int main(int argc, char **argv)
Definition JEffectiveMassCountingOffline.cc:84
JManager.hh
Dynamic ROOT object management.
setLongprint
void setLongprint(std::ostream &out)
Set long print option.
Definition JManip.hh:132
JMessage.hh
General purpose messaging.
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
JMultipleFileScanner.hh
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JParticleTypes.hh
Definition of particle types.
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
JSupport.hh
ROOT TTree parameter settings of various packages.
JVisibleEnergyToolkit.hh
Auxiliary methods for evaluating visible energies.
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