Program to histogram event-by-event data of muon 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/JPDFTransformer.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JConstants.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JIO/JFileStreamIO.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.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 "JSirene/JVisibleEnergyToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

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

Detailed Description

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

Author: mdejong

Definition in file JHistHDF.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 103 of file JHistHDF.cc.

 {
   using namespace std;
   using namespace JPP;
   using namespace KM3NET;  
 
   JMultipleFileScanner<Evt> inputFile;
   JLimit_t&   numberOfEvents = inputFile.getLimit();
   string      outputFile;
   string      detectorFile;
   int         debug;
 
   try {
 
     JParser<> zap("Program to histogram event-by-event data of muon 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['d'] = make_field(debug)               = 1;
 
     zap(argc, argv);
   }
   catch(const exception &error) {
     FATAL(error.what() << endl);
   }
 
 
   JDetector detector;
 
   try {
     load(detectorFile, detector);
   }
   catch(const JException& error) {
     FATAL(error);
   }
 
   JHead head(JMultipleFileScanner<Head>(inputFile).getHeader());  // Monte Carlo header
 
   const bool KM3Sim = is_km3sim(head);
   
   if (KM3Sim) {
     detector -= JCenter3D(detector.begin(), detector.end());      // For KM3Sim, you get it from the detector file
   } else {
     detector -= get<JPosition3D>(head);                           // For KM3, the MC origin is stored in the header
   }
 
   const JCylinder3D can = get<JCylinder3D>(head);
 
   const JPMTRouter pmtRouter(detector);
 
   const double P_atm = NAMESPACE::getAmbientPressure();
   const double wmax  = getMaximalWavelength();
     
   const JDispersion dispersion(P_atm);
 
   typedef JHistogram1D_t::abscissa_type                             abscissa_type;
 
   typedef JTransformableMultiHistogram<JHistogram1D_t, 
                                        JMAPLIST<JHistogramMap_t,
                                                 JHistogramGridMap_t,
                                                 JHistogramGridMap_t>::maplist> JMultiHistogram_t;
 
   typedef JPDFTransformer<3, abscissa_type>                         JFunction3DTransformer_t;
 
   JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)
   JMultiHistogram_t h1;   // light from muon
   JMultiHistogram_t h2;   // light from showers
 
   const double cmin = dispersion.getKmin (wmax);
 
   h1.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
   h2.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
 
 
   // configure bins
 
   const double R[] = { 0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 70.0, 80.0, 90.0, 100.0, 170.0, 250.0 };
 
   for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
 
     const double R_m = R[i];
     
     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, &h2 }) {
           (*p)[R_m][theta][phi];
         }
       }
     }
   }
 
   for (JMultiHistogram_t::super_iterator 
          i1 = h1.super_begin(),
          i2 = h2.super_begin(); i1 != h1.super_end(); ++i1, ++i2) {
 
     for (double x : { 0.0, 1.0, 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, 125.0, 150.0, 200.0, 500.0} ) {
       i1.getValue()[x];
       i2.getValue()[x];
     }
   }
 
 
   while (inputFile.hasNext()) {
 
     NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
 
     const Evt* event = inputFile.next();
 
     for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
 
       if (is_muon(*track)) {
 
         // track parameters
 
         const double Evis       = getVisibleEnergy      (*event, can);
         const Vec    EvisVector = getVisibleEnergyVector(*event, can);
 
         JVertex3D vertex = getVertex(*track);
 
         const JRotation3D R(getDirection(EvisVector));
 
         vertex.rotate(R);
 
         JRange<double> Z(0.0, gWater(track->E));   // range of muon
 
         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);
 
             const double t1 = vertex.getT() + (axis.getZ() + axis.getX()*getTanThetaC()) / getSpeedOfLight();
         
             const double R_m   = axis.getX();
             const double theta = axis.getTheta();
             const double phi   = fabs(axis.getPhi());
             const double dt    = getTime(*hit) - t1;
             const double npe   = getNPE (*hit);
 
             const int classID = (KM3Sim ? classify_km3sim(*track, *hit) : classify_km3(*track, *hit));
             
             switch (classID) {
 
             case 1:
               h1.fill(R_m, theta, phi, dt, npe);              // light from muon
               break;
 
             case 2:
               h2.fill(R_m, theta, phi, dt, npe / max(1.0, Evis)); // light from shower
               break;
             }
           }
           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 (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= h0.getXmax()) {
 
             for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
           
               JAxis3D axis = *pmt;
 
               axis.transform(R, vertex);
 
               if (Z(axis.getZ() - axis.getX()/getTanThetaC())) {
                 h0.fill(axis.getX(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
               }
             }
           }
         }
       }
     }
   }
   NOTICE(endl);
 
 
   // output
 
   JFileStreamWriter out(outputFile.c_str());
 
   for (const JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
     out.store(*p);
   }
 
   out.close();
 }

Functions

Detailed Description

Function Documentation