JRadiation.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "JTools/JConstants.hh"
#include "JTools/JQuantile.hh"

#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JGeane.hh"

#include "JSirene/JSeaWater.hh"
#include "JLang/JSharedPointer.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  static const std::string LENGTH = "length";  //!< calculation of interaction length
  static const std::string ENERGY = "energy";  //!< calculation of shower energy
  static const std::string RANGE  = "range";   //!< calculation of muon range
  static const std::string ELOSS  = "eloss";   //!< calculation of b parameter in dE/dx formula

  /**
   * Clone master histogram.
   *
   * \param  h1       pointer to master histogram
   * \param  buffer   name of clone
   * \return          pointer to cloned histogram
   */
  inline TH1* clone(TH1* h1, const char* buffer)
  {
    if (h1 != NULL) 
      return (TH1*) h1->Clone(buffer);
    else
      return NULL;
  }
}


/**
 * \file
 *
 * Example program to histogram radiation cross sections, shower energy, range and b(E)
 * as a function of the muon energy.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;

  string  outputFile;
  int     numberOfPoints;
  string  option;
  int     debug;

  try {

    JParser<> zap("Example program to histogram radiation cross sections, shower energy, range and b(E).");

    zap['o'] = make_field(outputFile)      = "radiation.root";
    zap['n'] = make_field(numberOfPoints)  = 0;
    zap['O'] = make_field(option)          = LENGTH, ENERGY, RANGE, ELOSS;
    zap['d'] = make_field(debug)           = 0;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPP;


  TFile out(outputFile.c_str(), "recreate");

  TH1* h0 = NULL;                                                               // master histogram

  if (option == LENGTH)  {
    h0 = new TH1D("total", NULL, 160, 2.0, 10.0);                               // interaction length [m]
  }

  if (option == ENERGY) {
    h0 = new TH1D("total", NULL,  24, 2.0, 10.0);                               // <E> [GeV]
  }

  NOTICE("Setting up radiation tables... " << flush);

  typedef pair<JSharedPointer<JRadiationInterface>, TH1*> 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);

  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));

  JRadiationSource::source_type EErad(&JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad);
  JRadiationSource::source_type Brems(&JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems);
  JRadiationSource::source_type GNrad(&JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad);

  radiation.push_back(tuple(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  EErad), clone(h0, "[eerad O]" )));
  radiation.push_back(tuple(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad), clone(h0, "[eerad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  EErad), clone(h0, "[eerad H]" )));
  
  radiation.push_back(tuple(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  Brems), clone(h0, "[Brems O]" )));
  radiation.push_back(tuple(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems), clone(h0, "[Brems Cl]")));
  radiation.push_back(tuple(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  Brems), clone(h0, "[Brems H]" )));

  radiation.push_back(tuple(new JRadiationSource(Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  GNrad), clone(h0, "[gnrad O]" )));
  radiation.push_back(tuple(new JRadiationSource(Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad), clone(h0, "[gnrad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  GNrad), clone(h0, "[gnrad H]" )));

  NOTICE("OK" << endl);

  
  if (option == LENGTH) {

    for (int i = 1; i <= h0->GetNbinsX(); ++i) {
    
      const double x = h0->GetBinCenter(i);
      const double E = pow(10.0, x); 
    
      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {

        const double li = p->first->getInverseInteractionLength(E);

        p->second->Fill(x, 1.0/li);
      }
    }
    STATUS(endl);
  }


  if (option == ENERGY) {

    for (int i = 1; i <= h0->GetNbinsX(); ++i) {
    
      const double x = h0->GetBinCenter(i);
      const double E = pow(10.0, x); 

      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {

        JQuantile Q;

        for (int n = 0; n != numberOfPoints; ++n) {
          Q.put(p->first->getEnergyOfShower(E));
        }
        
        p->second->SetBinContent(i, Q.getMean());
        //p->second->SetBinError  (i, Q.getSTDev());
      }
    }
    STATUS(endl);
  }


  if (option == RANGE) {
    
    TH1D* Ra = new TH1D("R[analytical]", NULL, 12, 2.0, 8.0);                                            // Range [m]
    TH1D* Rb = new TH1D("R[numerical]",  NULL, 12, 2.0, 8.0);                                            // Range [m]
    
    for (int i = 1; i <= Ra->GetNbinsX(); ++i) {

      const double x  = Ra->GetBinCenter(i);
      const double E0 = pow(10.0, x); 
      const double Z  = gWater(E0);
      
      Ra->SetBinContent(i, Z*1e-3);
    }
    
    for (int j = 1; j <= Rb->GetNbinsX(); ++j) {

      const double x  = Rb->GetBinCenter(j);
      const double E0 = pow(10.0, x); 
    
      STATUS("Energy: " << SCIENTIFIC(8,1) << E0 << '\r'); DEBUG(endl);

      JQuantile Q;

      for (int n = 0; n != numberOfPoints; ++n) {
      
        double E = E0;
        double z = 0.0;
      
        while (E >= 0.5) {
        
          const int N = radiation.size();
        
          double li[N];                                       // inverse interaction lengths
          double ls = 1.0e-5;                                 // minimal total inverse interaction length (100 km)^-1
        
          for (int i = 0; i != N; ++i) {
            ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
          }
        
          double dz = min(gRandom->Exp(1.0) / ls, gWater(E));

          E -= gWater.getA() * dz;
        
          double Es = 0.0;                                            // shower energy [GeV]
          double y  = gRandom->Uniform(ls);
        
          for (int i = 0; i != N; ++i) {
          
            y -= li[i];
          
            if (y < 0.0) {
              Es = radiation[i].first->getEnergyOfShower(E);          // shower energy [GeV]
              break;
            }
          }
        
          z += dz;
          E -= Es;
        }

        Q.put(z);
      }

      Rb->SetBinContent(j, Q.getMean() * 1e-3);
      //Rb->SetBinError  (j, Q.getSTDev()  * 1e-3);
    }
    STATUS(endl);
  }
  
    
  if (option == ELOSS) {

    TH1D* hb = new TH1D("hb", NULL, 12, 2.0, 8.0);                                            // b(E)

    for (int j = 1; j <= hb->GetNbinsX(); ++j) {

      const double x  = hb->GetBinCenter(j);
      const double E  = pow(10.0, x); 

      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      JQuantile Q;

      const int N = radiation.size();
    
      double li[N];                                       // inverse interaction lengths
      double ls = 1.0e-5;                                 // minimal total inverse interaction length (100 km)^-1
    
      for (int i = 0; i != N; ++i) {
        ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
      }
        
      for (int i = 0; i != numberOfPoints; ++i) {
      
        double Es = 0.0;                                            // shower energy [GeV]
        double y  = gRandom->Uniform(ls);
        
        for (int i = 0; i != N; ++i) {
        
          y -= li[i];
        
          if (y < 0.0) {
            Es = radiation[i].first->getEnergyOfShower(E);          // shower energy [GeV]
            break;
          }
        }

        Q.put(1e3*ls*Es/E);
      }

      hb->SetBinContent(j, Q.getMean());
      //hb->SetBinError  (j, Q.getSTDev());
    }
    STATUS(endl);
  }

  delete h0;

  out.Write();
  out.Close();
}