Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JShowerMCEvt.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 
7 #include "JAAnet/JHead.hh"
8 #include "JAAnet/JHeadToolkit.hh"
9 #include "JDAQ/JDAQEventIO.hh"
10 #include "JDAQ/JDAQTimesliceIO.hh"
17 #include "JSupport/JSupport.hh"
18 #include "JSupport/JMeta.hh"
19 #include "JPhysics/JGeanz.hh"
20 
21 #include "JReconstruction/JEvt.hh"
23 
24 #include "Jeep/JParser.hh"
25 #include "Jeep/JMessage.hh"
26 
27 
28 /**
29  * \file
30  *
31  * Auxiliary program to store Monte Carlo information from a neutrino or the primary electron in JFIT::JEvt format.
32  * By default the neutrino is considered, a bool variable allows to select the primary electron instead.
33  *
34  * \author mdejong, adomi
35  */
36 int main(int argc, char **argv)
37 {
38  using namespace std;
39  using namespace JPP;
40  using namespace KM3NETDAQ;
41 
43 
44  JTriggeredFileScanner<> inputFile;
46  JLimit_t& numberOfEvents = inputFile.getLimit();
47  int debug;
48  bool take_electron;
49 
50  try {
51 
52  JParser<> zap("Auxiliary program to store Monte Carlo true muon in format for subsequent fits.");
53 
54  zap['f'] = make_field(inputFile);
55  zap['o'] = make_field(outputFile) = "mcevt.root";
56  zap['n'] = make_field(numberOfEvents) = JLimit::max();
57  zap['d'] = make_field(debug) = 2;
58  zap['e'] = make_field(take_electron);
59 
60  zap(argc, argv);
61  }
62  catch(const exception& error) {
63  FATAL(error.what() << endl);
64  }
65 
66 
67  const JVector3D center = get<JPosition3D>(getHeader(inputFile));
68 
69  outputFile.open();
70 
71  outputFile.put(JMeta(argc, argv));
72 
73  while (inputFile.hasNext()) {
74 
75  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
76 
77  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
78 
79  const JDAQEvent* tev = ps;
80  const Evt* event = ps;
81  JEvt out;
82 
83  const JTimeConverter converter(*event, *tev);
84 
86 
87  if(!take_electron) fermion = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_neutrino);
88  else fermion = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_electron);
89 
90  if (fermion != event->mc_trks.end()) {
91 
92  // if fermion = neutrino => fermion position = neutrino interaction vertex
93  double time = converter.putTime();
94  JVector3D position = getPosition(*fermion);
95 
96  // if fermion = electron => fermion position and time = EM shower brightest point
97  if(take_electron){
98  JVector3D electron_dir = getDirection(*fermion);
99  double shower_elongation = geanz.getMaximum(fermion->E);
100  electron_dir *= shower_elongation;
101  position.add(electron_dir);
102  time += (shower_elongation * getInverseSpeedOfLight());
103  }
104 
105  JTrack3D td(position, getDirection(*fermion), time);
106  JTrack3E ta(td, fermion->E);
107 
108  ta.add(center);
109 
110  out.push_back(getFit(JMCEVT, ta, 0, 0.0, ta.getE()));
111 
112  }
113 
114  outputFile.put(out);
115  }
116  STATUS(endl);
117 
119 
120  io >> outputFile;
121 
122  outputFile.close();
123 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1493
ROOT TTree parameter settings.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
double putTime() const
Get Monte Carlo minus DAQ/trigger hit time.
#define STATUS(A)
Definition: JMessage.hh:63
Recording of objects on file according a format that follows from the file name extension.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
string outputFile
3D track with energy.
Definition: JTrack3E.hh:29
JTime & add(const JTime &value)
Addition operator.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Type list.
Definition: JTypeList.hh:22
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
double getE() const
Get energy.
Definition: JTrack3E.hh:92
Data structure for vector in three dimensions.
Definition: JVector3D.hh:33
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
ROOT I/O of application specific meta data.
int debug
debug level
Definition: JSirene.cc:61
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=0)
Get fit.
Definition: JEvtToolkit.hh:126
static const int JMCEVT
Data structure for set of track fit results.
Definition: JEvt.hh:294
double getMaximum(const double E) const
Get depth of shower maximum.
Definition: JGeanz.hh:162
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
JDirection3D getDirection(const Vec &v)
Get direction.
Auxiliary class to convert DAQ/trigger hit time to/from Monte Carlo hit time.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
Longitudinal emission profile EM-shower.
General purpose class for multiple pointers.
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:141
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])
Definition: Main.cpp:15