Program to histogram event-by-event data of shower light for making PDFs. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JConstants.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JIO/JFileStreamIO.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JTools/JHistogram1D_t.hh"
#include "JTools/JHistogramMap_t.hh"
#include "JTools/JTransformableMultiHistogram.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDF.hh"
#include "JSirene/JVisibleEnergyToolkit.hh"

Functions
int	main (int argc, char **argv)

Detailed Description

Program to histogram event-by-event data of shower light for making PDFs.

Author: lquinn

Definition in file JHistHDG.cc.

Function Documentation

◆ main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 56 of file JHistHDG.cc.

 {
   using namespace std;
   using namespace JPP;
  
   JMultipleFileScanner<Evt> inputFile;
   JLimit_t&   numberOfEvents = inputFile.getLimit();
   string      outputFile;
   string      detectorFile;
   bool        hadronicMode;
   int         debug;
  
   try {
  
     JParser<> zap("Program to histogram event-by-event data of shower light for making PDFs.");
  
     zap['f'] = make_field(inputFile);
     zap['o'] = make_field(outputFile);
     zap['n'] = make_field(numberOfEvents)         = JLimit::max();
     zap['a'] = make_field(detectorFile);
     zap['H'] = make_field(hadronicMode);
     zap['d'] = make_field(debug)                  = 1;
  
     zap(argc, argv);
   }
   catch(const exception &error) {
     FATAL(error.what() << endl);
   }
  
   if(hadronicMode) {
     NOTICE("Plotting hadronic showers." << endl);
   } else {
     NOTICE("Plotting EM showers." << endl);
   }
  
  
   JDetector detector;
  
   try {
     load(detectorFile, detector);
   }
   catch(const JException& error) {
     FATAL(error);
   }
  
   JHead head(JMultipleFileScanner<Head>(inputFile).getHeader());  // Monte Carlo header
  
   const Vec  offset = getOffset(head);
  
   NOTICE("Apply detector offset " << offset << endl);
  
   detector -= getPosition(offset);
  
   const JCylinder3D can(detector.begin(), detector.end());
  
   const JPMTRouter pmtRouter(detector);
  
   const double P_atm = NAMESPACE::getAmbientPressure();
   const double wmin  = getMinimalWavelength();
   const double wmax  = getMaximalWavelength();
  
   const JDispersion dispersion(P_atm);
  
   const double ng[] = { dispersion.getIndexOfRefractionGroup(wmax),
                         dispersion.getIndexOfRefractionGroup(wmin) };
  
   typedef JHistogram1D_t::abscissa_type                             abscissa_type;
  
   typedef JTransformableMultiHistogram<JHistogram1D_t, 
                                        JMAPLIST<JHistogramMap_t,
                                                 JHistogramMap_t,
                                                 JHistogramGridMap_t,
                                                 JHistogramGridMap_t>::maplist> JMultiHistogram_t;
  
   typedef JPDFTransformer<4, abscissa_type>                         JFunction4DTransformer_t;
  
   JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)
   JMultiHistogram_t h1;   // light from cascade
  
   h1.transformer.reset(new JFunction4DTransformer_t(21.5, 2, ng[0], 0.0,   JGeant(JGeanx(1.00, -2.2)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));
  
   set<double> C;  // cosine emission angle
  
   JQuadrature qeant(-1.0, 1.0, 20, JGeanx(1.00, -2.2));
  
   for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
     C.insert(i->getX());
   }
  
   C.insert(-1.01);
   C.insert(-1.00);
   C.insert( 1.00);
   C.insert( 1.01);
  
   const double R[] = {0.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 60.0, 70.0, 80.0};
   const double Dmax_m = 80.0;
  
   for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
  
     const double R_m = R[i];
     
     for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {
  
       const double cd = *c;
  
       const double grid  = 10.0 +  0.0 * R_m/100.0;                         // [deg]
       const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));        // azimuth angle unit step size
     
       const int    number_of_theta_points = max(2, (int) (180.0/(1.4 * grid)));
       const double theta_step             = PI / (number_of_theta_points + 1);
  
       for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
         
         const int    number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha));
         const double phi_step             = PI / (number_of_phi_points + 1);
         
         for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
     
           for (JMultiHistogram_t* p : { &h0, &h1 }) {
             (*p)[R_m][cd][theta][phi];
           }
         }
       }
     }
   }
  
   double buffer[] = { 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 7.5, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 85.0, 100.0 };
  
   for (JMultiHistogram_t::super_iterator i1 = h1.super_begin(); i1 != h1.super_end(); ++i1) {
     for (int j = 0; j != sizeof(buffer)/sizeof(buffer[0]); ++j) {
       i1.getValue()[buffer[j]];
     }
   }
  
   while (inputFile.hasNext()) {
  
     NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
  
     const Evt* event = inputFile.next();
  
     if (!has_neutrino(*event)) {
       WARNING("Event " << inputFile.getCounter() << " does not correspond to a neutrino interaction; skip.");
       continue;
     }
     
     if (!has_electron(*event) || count_hadrons(*event) == 0) {
       WARNING("No electron/hadrons found; skip.");
       continue;
     }
  
     const Trk& neutrino = get_neutrino(*event);
  
     const double Evis       = getVisibleEnergy      (*event, can);
     const Vec    EvisVector = getVisibleEnergyVector(*event, can);
  
     JVertex3D vertex = getVertex(neutrino);
     
     const JRotation3D R(getDirection(EvisVector));
     
     vertex.rotate(R);
  
     for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit)  {
  
       try {
  
         JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
  
         axis.transform(R, vertex.getPosition());
  
         if ((!hadronicMode && from_electron(*hit)) || (hadronicMode && from_hadron(*hit))) {
       
           const double t1    = vertex.getT() + axis.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
       
           const double D_m   = axis.getLength();
           const double cd    = axis.getZ() / D_m;
           const double theta = axis.getTheta();
           const double phi   = fabs(axis.getPhi());
           const double dt    = getTime(*hit) - t1;
           const double npe   = getNPE (*hit);
  
           if (D_m < Dmax_m) {
             h1.fill(D_m, cd, theta, phi, dt, npe/Evis);
           }
         }
       }
       catch(const exception& error) {
         FATAL(error.what() << endl);
       }
     }
  
     for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
  
       JPosition3D P = module->getPosition();
  
       P.rotate(R);
  
       P.sub(vertex);
  
       if (P.getLength() < h0.getXmax()) {
  
         for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
  
           JAxis3D axis = *pmt;
  
           axis.transform(R, vertex);
  
           h0.fill(axis.getLength(), axis.getZ()/axis.getLength(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
         }
       }
     }
   }
  
   double integral = 0;
   for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {
     integral+=i.getValue().getIntegral();
   }
   DEBUG("Integral:\t" << integral << endl);
  
   // output
  
   JFileStreamWriter out(outputFile.c_str());
  
   NOTICE("Storing, " << flush);
  
   for (const JMultiHistogram_t* p : { &h0, &h1 }) {
     out.store(*p);
   }
  
   out.close();
   NOTICE("done." << endl);
 }

Functions

Detailed Description

Function Documentation

◆ main()