Jpp
JORCAShowerPositionPDF.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <cmath>
3 
4 #include "evt/Evt.hh"
5 
6 #include "JDAQ/JDAQEvent.hh"
7 
8 #include "Jeep/JParser.hh"
9 #include "Jeep/JMessage.hh"
10 
11 #include "JFit/JPoint4D.hh"
12 
13 #include "JLang/JTypeList.hh"
14 
15 #include "JDAQ/JDAQEvent.hh"
16 
17 #include "JSupport/JSupport.hh"
19 #include "JSupport/JTreeScanner.hh"
20 
22 #include "JTrigger/JHitL0.hh"
23 #include "JTrigger/JHitL1.hh"
24 #include "JTrigger/JBuildL0.hh"
25 
26 #include "JDetector/JDetector.hh"
28 #include "JDetector/JPMTRouter.hh"
30 
33 
34 #include <TFile.h>
35 #include <TH2D.h>
36 
37 int main(int argc, char **argv){
38 
39  using namespace JPP;
40 
41  typedef JTriggeredFileScanner<> JTriggeredFileScanner_t;
42  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
43 
44  JTriggeredFileScanner_t inputFile;
45  string outputFile;
46  string detectorFile;
47  JRange<double> t_res_ns;
48  JRange<double> D_m;
49  JRange<double> cosT;
50  bool PMTweight; // weight PDF to take into account PMT growing number
51  // with respect to vertex
52  int debug;
53 
54  try{
55  JParser<> zap;
56 
57  zap['f'] = make_field(inputFile);
58  zap['o'] = make_field(outputFile) = "JORCAShowerPositionFit_PDF.root";
59  zap['a'] = make_field(detectorFile);
60  zap['T'] = make_field(t_res_ns) = JRange<double>(-100, 100); // [ns]
61  zap['D'] = make_field(D_m) = JRange<double>(0, 100); // [m]
62  zap['C'] = make_field(cosT) = JRange<double>(-1, 1);
63  zap['d'] = make_field(debug) = 2;
64  zap['P'] = make_field(PMTweight);
65 
66  if (zap.read(argc, argv) != 0) return 1;
67  }
68  catch(const exception& error){
69  FATAL(error.what() << endl);
70  }
71 
72  cout.tie(&cerr);
73 
74  JDetector detector;
75 
76  try{
77  load(detectorFile, detector);
78  }
79  catch(const JException& error){
80  FATAL(error);
81  }
82 
83  const JPMTRouter pmtRouter(detector);
84  const JModuleRouter moduleRouter(detector);
85 
86  typedef vector<JHitL0> JDataL0_t;
87  const JBuildL0<JHitL0> buildL0;
88 
89  TFile *out = new TFile(outputFile.c_str(), "RECREATE");
90  if ( !out->IsOpen() ) printf("Can't open outputfile\n");
91 
92  TH2D hPDF2Dist("hPDF2Dist", "PDF for ORCA Shower Position Fit; #delta t [ns]; D [m]",
93  2000, t_res_ns.getLowerLimit(), t_res_ns.getUpperLimit(),
94  700, D_m.getLowerLimit(), D_m.getUpperLimit());
95 
96  while(inputFile.hasNext()){
97 
98  multi_pointer_type ps = inputFile.next();
99 
100  JDAQEvent* tev = ps;
101  Evt* event = ps;
102  JTimeConverter converter(*event, *tev);
103 
104  // Get MC truth neutrino and shower
105  Trk cascade;
106  cascade.pos = event->mc_trks[1].pos;
107  cascade.pos -= Vec(0.0, 0.0, -117.164);
108 
109  for (unsigned int i = 1; i < event->mc_trks.size(); i++) {
110  cascade.E += event->mc_trks[i].E;
111  cascade.dir += event->mc_trks[i].dir * event->mc_trks[i].E;
112  }
113  cascade.dir.normalize();
114 
115  JPosition3D vx = getPosition(cascade.pos);
116  double t_vx = converter.putTime(event->mc_trks[1].t);
117  JDirection3D dir = getDirection(cascade.dir);
118 
119  JDataL0_t dataL0;
120  buildL0(JDAQTimeslice(*tev, true), moduleRouter, back_inserter(dataL0));
121 
122  //**************** PDF with real interaction vertex ***************
123  const JRotation3D R(dir);
124  vx.rotate(R);
125 
126  for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
127 
128  JHitL0 hit(*i);
129  hit.transform(R, vx);
130 
131  const double D = hit.getLength();
132  const double t_expected = t_vx + D * getInverseSpeedOfLight() * getIndexOfRefraction();
133  const double dt = hit.getT() - t_expected;
134  const JDirection3D photonDir(hit.getPosition());
135  const double ct = photonDir.getDot(hit.getDirection());
136 
137  if((D >= D_m.getLowerLimit() && D <= D_m.getUpperLimit()) &&
138  (dt >= t_res_ns.getLowerLimit() && dt <= t_res_ns.getUpperLimit()) &&
139  (ct >= cosT.getLowerLimit() && ct <= cosT.getUpperLimit())){
140 
141  double w = 1;
142 
143  if(PMTweight){
144  w = 1 / (D*D);
145  }
146 
147  hPDF2Dist.Fill(dt, D, w);
148  }
149  }
150  }
151  STATUS(endl);
152 
153  out->cd();
154  hPDF2Dist.Write();
155  out->Close();
156 
157 }
main
int main(int argc, char **argv)
Definition: JORCAShowerPositionPDF.cc:37
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:34
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:708
JMessage.hh
JPosition3D.hh
JDirection3D.hh
JPoint4D.hh
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JTimeConverter.hh
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:36
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
JTriggeredFileScanner.hh
JTreeScanner.hh
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JBuildL0.hh
JTOOLS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
debug
int debug
debug level
Definition: JSirene.cc:59
JTOOLS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
JModuleRouter.hh
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JParser.hh
JDetectorToolkit.hh
JHitL1.hh
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
JPMTRouter.hh
JTypeList.hh
JDAQEvent.hh
JDetector.hh
JHitL0.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:1791