JRadiation2D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <memory>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "TProfile2D.h"
 
#include "JTools/JQuantile.hh"
 
#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JConstants.hh"
 
#include "JSirene/JSeaWater.hh"
#include "JLang/JSharedPointer.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to histogram shower energy as a function of the muon energy.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;
 
  string  outputFile;
  int     numberOfPoints;
  int     debug;
 
  try {
 
    JParser<> zap("Example program to histogram shower energy.");
 
    zap['o'] = make_field(outputFile)      = "radiation.root";
    zap['n'] = make_field(numberOfPoints)  = 0;
    zap['d'] = make_field(debug)           = 0;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH2* h0 = new TH2D("total", NULL, 24, 2.0, 10.0, 20, -5.0, 0.0);              // master histogram
  TH2* h2 = (TH2*) h0->Clone("h2");
 
  NOTICE("Setting up radiation tables... " << flush);
 
  typedef pair<JSharedPointer<JRadiationInterface>, TH2*> tuple;
  typedef vector<tuple>                                   ntuple;
 
  ntuple radiation;
 
  const JRadiation hydrogen ( 1.0,  1.0, 40, 0.01, 0.1, 0.1);
  const JRadiation oxygen   ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
  const JRadiation chlorine (17.0, 35.0, 40, 0.01, 0.1, 0.1);
  const JRadiation sodium   (11.0, 23.0, 40, 0.01, 0.1, 0.1);
 
  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Oxygen   (new JRadiationFunction(oxygen,   300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Sodium   (new JRadiationFunction(sodium,   300, 0.2, 1.0e11));
 
  radiation.push_back(tuple(new JRadiationSource(11, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  EErad_t), (TH2D*) h0->Clone("[eerad O]" )));
  radiation.push_back(tuple(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t), (TH2D*) h0->Clone("[eerad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  EErad_t), (TH2D*) h0->Clone("[eerad H]" )));
  radiation.push_back(tuple(new JRadiationSource(14, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t), (TH2D*) h0->Clone("[eerad Na]" )));
  
  radiation.push_back(tuple(new JRadiationSource(21, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  Brems_t), (TH2D*) h0->Clone("[Brems O]" )));
  radiation.push_back(tuple(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t), (TH2D*) h0->Clone("[Brems Cl]")));
  radiation.push_back(tuple(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  Brems_t), (TH2D*) h0->Clone("[Brems H]" )));
  radiation.push_back(tuple(new JRadiationSource(24, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t), (TH2D*) h0->Clone("[Brems Na]" )));
 
  radiation.push_back(tuple(new JRadiationSource(31, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  GNrad_t), (TH2D*) h0->Clone("[gnrad O]" )));
  radiation.push_back(tuple(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t), (TH2D*) h0->Clone("[gnrad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  GNrad_t), (TH2D*) h0->Clone("[gnrad H]" )));
  radiation.push_back(tuple(new JRadiationSource(34, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t), (TH2D*) h0->Clone("[gnrad Na]" )));
 
  NOTICE("OK" << endl);
 
  for (int i = 1; i <= h0->GetXaxis()->GetNbins(); ++i) {
    
    const double x = h0->GetXaxis()->GetBinCenter(i);
    const double E = pow(10.0, x);
    const double W = 1.0 / (double) numberOfPoints;
 
    STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);
 
    for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {
 
      for (int n = 0; n != numberOfPoints; ++n) {
 
        const double li = p->first->getInverseInteractionLength(E);
        const double Es = p->first->getEnergyOfShower(E);
        const double y  = log10(Es/E);
 
        p->second->Fill(x, y, W);
 
        h2->Fill(x, y, li * W);
      }
    }
    STATUS(endl);
  }
 
  delete h0;
 
  out.Write();
  out.Close();
}