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 "JPhysics/JSeaWater.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<shared_ptr<JRadiationInterface>, TH2*> tuple;
  typedef vector<tuple>                               ntuple;
 
  ntuple radiation;
 
  const JRadiation hydrogen (JSeaWater::H .Z, JSeaWater::H .A, 40, 0.01, 0.1, 0.1);
  const JRadiation oxygen   (JSeaWater::O .Z, JSeaWater::O .A, 40, 0.01, 0.1, 0.1);
  const JRadiation chlorine (JSeaWater::Cl.Z, JSeaWater::Cl.A, 40, 0.01, 0.1, 0.1);
  const JRadiation sodium   (JSeaWater::Na.Z, JSeaWater::Na.A, 40, 0.01, 0.1, 0.1);
 
  shared_ptr<JRadiation> Hydrogen (make_shared<JRadiationFunction>(hydrogen, 300, 0.2, 1.0e11));
  shared_ptr<JRadiation> Oxygen   (make_shared<JRadiationFunction>(oxygen,   300, 0.2, 1.0e11));
  shared_ptr<JRadiation> Chlorine (make_shared<JRadiationFunction>(chlorine, 300, 0.2, 1.0e11));
  shared_ptr<JRadiation> Sodium   (make_shared<JRadiationFunction>(sodium,   300, 0.2, 1.0e11));
 
  radiation.push_back(tuple(make_shared<JRadiationSource>(11, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  JRadiation::EErad_t), (TH2D*) h0->Clone("[eerad O]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::EErad_t), (TH2D*) h0->Clone("[eerad Cl]")));
  radiation.push_back(tuple(make_shared<JRadiationSource>(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  JRadiation::EErad_t), (TH2D*) h0->Clone("[eerad H]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(14, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::EErad_t), (TH2D*) h0->Clone("[eerad Na]")));
  
  radiation.push_back(tuple(make_shared<JRadiationSource>(21, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  JRadiation::Brems_t), (TH2D*) h0->Clone("[Brems O]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::Brems_t), (TH2D*) h0->Clone("[Brems Cl]")));
  radiation.push_back(tuple(make_shared<JRadiationSource>(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  JRadiation::Brems_t), (TH2D*) h0->Clone("[Brems H]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(24, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::Brems_t), (TH2D*) h0->Clone("[Brems Na]")));
 
  radiation.push_back(tuple(make_shared<JRadiationSource>(31, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  JRadiation::GNrad_t), (TH2D*) h0->Clone("[gnrad O]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), JRadiation::GNrad_t), (TH2D*) h0->Clone("[gnrad Cl]")));
  radiation.push_back(tuple(make_shared<JRadiationSource>(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  JRadiation::GNrad_t), (TH2D*) h0->Clone("[gnrad H]" )));
  radiation.push_back(tuple(make_shared<JRadiationSource>(34, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), JRadiation::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();
}