80{
   84 
   86  typedef JTriggeredFileScanner_t::multi_pointer_type       multi_pointer_type;
   87 
   88  JTriggeredFileScanner_t inputFile;
   91  size_t                  numberOfPrefits;
   94  double                  Emin_GeV;
   95  double                  NPE;
   96  int                     application;
   97  string                  option;
   98  bool                    isMuon;
   99  bool                    wrtNeutrino;
  100  double                  radius;
  103 
  104  try {
  105 
  106    JParser<> zap(
"Example program to histogram fit results.");
 
  107 
  117    zap[
'O'] = 
make_field(option)              = 
"E", 
"N", 
"LINE", 
"LOGE";
 
  120    zap[
'r'] = 
make_field(radius)              = numeric_limits<double>::max();     
 
  123 
  124    zap(argc, argv);
  125  }
  126  catch(const exception& error) {
  127    FATAL(error.what() << endl);
 
  128  }
  129  
  130  cout << "APPLICATION " << application << endl;
  131 
  133 
  134  try {
  136  }
  137  catch(const exception& error) {
  138    FATAL(error.what() << endl);
 
  139  }
  140 
  141  JVolume     volume(head, option != 
"LINE");
 
  145 
  146  cylinder.
add(offset);
 
  147 
  149 
  150  containment.add(origin);
  151 
  152  NOTICE(
"Offset: "      << offset      << endl);
 
  153  NOTICE(
"Origin: "      << origin      << endl);
 
  154  NOTICE(
"Cylinder:    " << cylinder    << endl);
 
  155  NOTICE(
"Containment: " << containment << endl);
 
  156 
  157  const double EMIN_GEV = 10e3;  
  158 
  160 
  161  TH1D job("job", NULL, 100, 0.5, 100.5);
  162 
  163  TH1D hn("hn",  "NDF (Number of Degrees of Freedom)", 250, -0.5, 249.5);  
  164  TH1D hq("hq",  "Fit Quality",  300,    0.0, 600.0);    
  165  TH1D hi("hi",  "Fit Index vs Best Resolution",  101, -0.5, 100.5);  
  166  TH1D hd("hd",  "Distance between Reco and True position",  2000,    0.0,  100.0);    
  167  TH1D hz("hz",  "Longitudinal Distance between Reco and MC lepton position",  100, -200.0, 200.0);    
  168  TH1D ht("ht",  "Time difference between Reco and MC lepton",  100, 0, 100.0);    
  169  TH1D h4D("h4D", "4D-distance between reco and true vertex", 2000, 0.0, 100.0); 
  170  
  171  TH1D ha("ha",  "Angle between reco direction and true direction",  1000, 0.0, 180.0);    
  172  
  173  TH1D e0("e0",  "True lepton energy",  100, volume.getXmin(), volume.getXmax());       
  174  TH1D e2("e2",  "Reconstructed energy",  100, volume.getXmin(), volume.getXmax());       
  175  TH1D er("er",  "Ratio of reconstructed energy and true energy",  100, -5.0, +5.0);       
  176  TH1D ea("ea", "Distribution of Energy error for all the events; E_{reco}-E_{true}", 100, -30, +30);
  177  TH1D ed3_5GeV("ed3_5GeV", "Distribution of Energy error for events with MC energy #in [3,5] GeV; E_{reco}-E_{true}", 100, -30, +30);
  178  TH1D ed8_11GeV("ed8_11GeV", "Distribution of Energy error for events with MC energy #in [8,11] GeV; E_{reco}-E_{true}", 100, -30, +30);
  179  TH1D ed15_21GeV("ed15_21GeV", "Distribution of Energy error for events with MC energy #in [15,21] GeV; E_{reco}-E_{true}", 100, -30, +30);
  180  TH2D ee("ee",  "; E_{true} [GeV]; E_{reco} [GeV]",
  181          100, volume.getXmin(), volume.getXmax(),
  182          100, volume.getXmin(), volume.getXmax());
  183  const Int_t    ny   = 28800;
  184  const Double_t ymin =     0.0;                 
  185  const Double_t ymax =   180.0;                 
  186  
  188 
  189  if (option.find('E') != string::npos) {
  190 
  191    for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
  192      X.push_back(x);
  193    }
  194 
  195  } else {
  196 
  198    
  199    for ( ; 
x <=  15.5; 
x +=  1.0) { X.push_back(x); }
 
  200    for ( ; 
x <=  25.5; 
x +=  2.0) { X.push_back(x); }
 
  201    for ( ; 
x <=  50.5; 
x +=  5.0) { X.push_back(x); }
 
  202    for ( ; 
x <= 100.5; 
x += 10.0) { X.push_back(x); }
 
  203    for ( ; 
x <= 250.5; 
x += 20.0) { X.push_back(x); }
 
  204  }
  205 
  206  TH2D     h2("h2", NULL, X.size() - 1, X.data(), ny,    ymin,    ymax);
  207  TH2D     h3("h3", NULL, X.size() - 1, X.data(), 2000,  0.0,     100.0);
  208  TProfile herrorE("herrorE", "Energy Error", X.size() - 1, X.data());
  209  TProfile hfracE( "hfracE", "Fractional Energy Error", X.size() - 1, X.data());
  210  TProfile he("he","Reconstruction Efficiency", X.size() - 1, X.data());
  211  TProfile htheta_nu_lep("htheta_nu_lep", "Angle between neutrino and primary lepton", X.size() - 1, X.data());
  212  TH1D hln1d("hln1d","Angle between neutrino and lepton",40,0,4);
  213  TH2D     hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
  214  TH2S     hmca("hmca", NULL, X.size() - 1, X.data(), 1000, 0, 180);
  215  TH2D hzenith("hzenith", "Reco Zenith vs MC Energy", 100, 0, 100, ny, ymin, ymax);
  216  TH2D hY("hY", "Reco Bjorken Y vs MC Bjorken Y", 40, 0, 1, 40, 0, 1);
  217  TH2D hby3_5GeV("hby3_5GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [3, 5] GeV", 40, 0, 1, 40, 0, 1);
  218  TH2D hby8_10GeV("hby8_10GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [8, 10] GeV", 40, 0, 1, 40, 0, 1);
  219 
  226  while (inputFile.hasNext()) {
  227 
  228    STATUS(
"event: " << setw(10) << inputFile.getCounter() << 
'\r'); 
DEBUG(endl);
 
  229 
  230    multi_pointer_type ps = inputFile.next();
  231 
  235 
  237 
  238    job.Fill(1.0);
  239 
  240    double Enu  = 0.0;
  241    double Elep = 0.0;
  242    double Emax = 0.0;
  243 
  245 
  247    
  250    
  251    }
  252    
  253    vector<Trk>::const_iterator lepton = event->mc_trks.end();  
  254 
  255    for (vector<Trk>::const_iterator mc_evt = event->mc_trks.begin(); mc_evt != event->mc_trks.end(); ++mc_evt) {
  256        
  258 
  259        Elep += mc_evt->E;
  260        
  261        if (mc_evt->E > Emax) {
  262          lepton = mc_evt;
  263          Emax = mc_evt->E;
  264        }
  265 
  266 
  267      } 
else if(isMuon && 
is_muon(*mc_evt)){
 
  268        
  269        Elep += mc_evt->E;
  270        if (mc_evt->E > Emax) {
  271          lepton = mc_evt;
  272          Emax = mc_evt->E;
  273        }
  274      }
  275    }
  276 
  277    if (lepton == event->mc_trks.end()) {
  278      continue;
  279    }
  280 
  281    job.Fill(3.0);
  282 
  283    double true_BjY = (Enu - Elep) / Enu;
  284    
  285    
  286 
  288    
  289    if (option.find('E') != string::npos){
  290 
  291      if(wrtNeutrino == 
true && !isMuon) 
x = volume.getX(neutrino.
E);
 
  292      else if(!wrtNeutrino && !isMuon)   
x = volume.getX(lepton->E);
 
  293    
  294    } else{
  296    }
  297    
  298    
  299    Double_t W = 0.0;
  300 
  301    if (!evt->empty()) {
  302 
  304 
  305      if (evt->begin() == __end) {
  306        continue;
  307      }
  308      
  309      job.Fill(4.0);
  310 
  311      
  312      if (numberOfPrefits > 0 ) {       
  313 
  314        JEvt::iterator __q = __end;
  315 
  316        advance(__end = evt->begin(), min(numberOfPrefits, (
size_t) 
distance(evt->begin(), __q)));
 
  317        
  319 
  320      } else {
  321 
  323      }
  324 
  325      JEvt::iterator best = evt->begin();  
  326 
  330 
  332 
  333      if(!wrtNeutrino && !isMuon){
  335      } else if(wrtNeutrino == true && !isMuon){
  339      }
  340 
  341      JEvt::iterator fit_with_smallest_error = best;
  342 
  344        fit_with_smallest_error  = position(evt->begin(), __end);
  346        fit_with_smallest_error  = energy(evt->begin(), __end);
  347      } else {
  348        fit_with_smallest_error  = pointing(evt->begin(), __end);
  349      }
  350      best = fit_with_smallest_error;
  351      
  352      const Double_t beta  = pointing.getAngle(*best);
  353      const double   Efit  = best->getE();
  354 
  356      
  357        
  358        bool ok = (Efit >= Emin_GeV);
  359      
  360        if (ok) {
  361 
  362          W = 1.0;
  363 
  364          job.Fill(5.0);
  365 
  366          hn.Fill((Double_t) best->getNDF());
  367          hq.Fill(best->getQ());
  368          hi.Fill((Double_t) 
distance(evt->begin(), fit_with_smallest_error));
 
  369          Q.put(beta);
  370 
  372 
  373          if(wrtNeutrino){
  375          } else {
  377          }
  378 
  380 
  382          
  384        
  385          
  387          tb.
sub(converter.putTime());
 
  388          double time_true = ta.
getT();
 
  389          double time_reco=tb.
getT();
 
  390          double distance_m;
  391          
  392          if(!isMuon && !wrtNeutrino){
  396            hd.Fill(fabs(distance_m));
  397            ht.Fill(fabs(time_reco - time_true));
  398            Qp.put(fabs(distance_m));
  399            Qt.put(fabs(time_reco-time_true));
  401            h4D.Fill(distance4d);
  402            Q4d.put(distance4d);
  403            h3.Fill(x, fabs(distance_m));
  404          } else {
  407            hd.Fill(fabs(distance_m));
  408            ht.Fill(fabs(time_reco - time_true));
  409            Qp.put(fabs(distance_m));
  410            Qt.put(fabs(time_reco-time_true));
  412            h4D.Fill(distance4d);
  413            Q4d.put(distance4d);
  414            h3.Fill(x, fabs(distance_m));
  415          }
  416        
  418          
  419            job.Fill(6.0);
  420 
  422          
  424 
  426 
  427              job.Fill(7.0);
  428 
  430            }
  431          }
  432        
  433          if (best->getE() >= EMIN_GEV) {
  436          
  438          }
  440 
  441            hY.Fill(true_BjY, best->getW(5));
  442 
  443            if(Efit > 2.5 && Efit < 6) hby3_5GeV.Fill(true_BjY, best->getW(5));
  444            else if(Efit > 7.5 && Efit < 11) hby8_10GeV.Fill(true_BjY, best->getW(5));
  445 
  446          }
  447          if(wrtNeutrino){
  448 
  449            e0.Fill(volume.getX(Enu,  true));
  450            er.Fill(volume.getX(Efit) - volume.getX(Enu));
  451            ee.Fill(volume.getX(Enu), volume.getX(Efit));
  452            ea.Fill(Efit - Enu);
  453            hzenith.Fill(Enu, zenith);
  454            QE.put(log10(Efit/Enu));
  455            herrorE.Fill(x, volume.getX(Efit) - volume.getX(Enu));
  456            hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Enu))/volume.getX(Enu));
  457            
  458            if(Enu >= 3 && Enu <= 5) ed3_5GeV.Fill(Efit - Enu);
  459            else if(Enu >= 8 && Enu <= 11) ed8_11GeV.Fill(Efit - Enu);
  460            else if(Enu >= 15 && Enu <= 21) ed15_21GeV.Fill(Efit - Enu);
  461 
  462          } else {
  463 
  464            e0.Fill(volume.getX(Elep,  true));
  465            er.Fill(volume.getX(Efit) - volume.getX(Elep));
  466            ee.Fill(volume.getX(Elep), volume.getX(Efit));
  467            ea.Fill(Efit - Elep);
  468            hzenith.Fill(Elep, zenith);
  469            QE.put(log10(Efit/Elep));
  470            herrorE.Fill(x, volume.getX(Efit) - volume.getX(Elep));
  471            hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Elep))/volume.getX(Elep));      
  472 
  473            if(Elep >= 3 && Elep <= 5) ed3_5GeV.Fill(Efit - Elep);
  474            else if(Elep >= 8 && Elep <= 11) ed8_11GeV.Fill(Efit - Elep);
  475            else if(Elep >= 15 && Elep <= 21) ed15_21GeV.Fill(Efit - Elep);
  476 
  477          }
  478          e2.Fill(volume.getX(Efit, true));
  479          h2.Fill(x, beta);
  480        }
  481      }
  482    }
  483    he.Fill(x, W);}
  484 
  486 
  487  NOTICE(
"Number of events input           " << setw(8) << right << job.GetBinContent(1) << endl);
 
  488  if(!isMuon){
  489    NOTICE(
"Number of events with electron   " << setw(8) << right << job.GetBinContent(3) << endl);
 
  490  } else{
  491    NOTICE(
"Number of events with muon     " << setw(8) << right << job.GetBinContent(3) << endl);
 
  492  }
  493  NOTICE(
"Number of events with fit        " << setw(8) << right << job.GetBinContent(4) << endl);
 
  494  NOTICE(
"Number of events selected        " << setw(8) << right << job.GetBinContent(5) << endl);
 
  495  NOTICE(
"Number of events with neutrino   " << setw(8) << right << job.GetBinContent(6) << endl);
 
  496  NOTICE(
"Number of events contained       " << setw(8) << right << job.GetBinContent(7) << endl);
 
  497 
  498  if (Q.getCount() != 0) {
  499    NOTICE(
"Median, 70, 80 and 90% quantiles of space angle [deg] " << 
FIXED     (6,3) << Q.getQuantile(0.5) << 
" "<< Q.getQuantile(0.7) << 
" "<< Q.getQuantile(0.8)<< 
" "<< Q.getQuantile(0.9)<< endl);
 
  500  }
  501  if (QE.getCount() != 0) {
  502      NOTICE(
"Median, 70%, 80% and 90% quantiles of energy ratio log(Ereco/Etrue) " << 
FIXED     (6,3) << QE.getQuantile(0.5) << 
" "<< QE.getQuantile(0.7) << 
" "<< QE.getQuantile(0.8)<< 
" "<< QE.getQuantile(0.9) <<  endl);
 
  503  }
  504  if(Qp.getCount() != 0){
  505        NOTICE(
"Median distance to vertex:      " << setw(8) << Qp.getQuantile(0.5) << 
" "<< Qp.getQuantile(0.7)<<
" "<< Qp.getQuantile(0.8)<<  
" "<<Qp.getQuantile(0.9)<<endl);
 
  506        NOTICE(
"Median time to vertex:      " << setw(8)  << Qt.getQuantile(0.5) << 
" "<< Qt.getQuantile(0.7)<<
" "<< Qt.getQuantile(0.8)<<  
" "<<Qt.getQuantile(0.9)<<endl);
 
  507        NOTICE(
"Median 4D-distance to vertex:      " << setw(8)  << Q4d.getQuantile(0.5) << 
" "<< Q4d.getQuantile(0.7)<<
" "<< Q4d.getQuantile(0.8)<<  
" "<<Q4d.getQuantile(0.9)<<endl);
 
  508  }
  509      out.Write();
  510      out.Close();
  511}
#define DEBUG(A)
Message macros.
 
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
 
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
 
Data structure for circle in two dimensions.
 
Data structure for position in two dimensions.
 
Data structure for vector in two dimensions.
 
double getLengthSquared() const
Get length squared.
 
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
 
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
 
JCylinder3D & add(const JVector3D &pos)
Add position.
 
Data structure for position in three dimensions.
 
const JPosition3D & getPosition() const
Get position.
 
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
 
JTime & add(const JTime &value)
Addition operator.
 
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
 
JTime & sub(const JTime &value)
Subtraction operator.
 
Data structure for vector in three dimensions.
 
JVector3D & add(const JVector3D &vector)
Add vector.
 
double getLength() const
Get length.
 
double getLengthSquared() const
Get length squared.
 
double getTheta() const
Get theta angle.
 
Utility class to parse command line options.
 
double getMaximum(const double E) const
Get depth of shower maximum.
 
Auxiliary class to evaluate atmospheric muon hypothesis.
 
Auxiliary class to compare fit results with respect to a reference direction (e.g....
 
Auxiliary class to compare fit results with respect to a reference position (e.g. true muon).
 
const_iterator< T > end() const
Get end of data.
 
const_iterator< T > begin() const
Get begin of data.
 
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
 
static const int JSHOWERPOSITIONFIT
 
static const int JSHOWERDIRECTIONPREFIT
 
static const int JSHOWERPOINTSIMPLEX
 
static const int JSHOWERCOMPLETEFIT
 
static const int JSHOWER_BJORKEN_Y
 
static const int JSHOWERPREFIT
 
static const int JSHOWERENERGYPREFIT
 
Vec getOrigin(const JHead &header)
Get origin.
 
JDirection3D getDirection(const Vec &dir)
Get direction.
 
JTrack3E getTrack(const Trk &track)
Get track.
 
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
 
JCylinder3D getCylinder(const JHead &header)
Get cylinder corresponding to the positions of generated tracks (i.e.
 
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
 
JPosition3D getPosition(const Vec &pos)
Get position.
 
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
 
Vec getOffset(const JHead &header)
Get offset.
 
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
 
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
 
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
 
T pow(const T &x, const double y)
Power .
 
static const double PI
Mathematical constants.
 
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
 
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
 
const double getInverseSpeedOfLight()
Get inverse speed of light.
 
const double getSpeedOfLight()
Get speed of light.
 
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
 
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
 
Head getCommonHeader(const JMultipleFileScanner_t &file_list)
Get common Monte Carlo header.
 
KM3NeT DAQ data structures and auxiliaries.
 
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
 
Auxiliary data structure for floating point format specification.
 
Auxiliary class for histogramming of effective volume.
 
Auxiliary class to test history.
 
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
 
General purpose sorter of fit results.
 
Auxiliary class for defining the range of iterations of objects.
 
const JLimit & getLimit() const
Get limit.
 
static counter_type max()
Get maximum counter value.
 
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
 
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
 
double E
Energy [GeV] (either MC truth or reconstructed)
 
Reconstruction type dependent comparison of track quality.