Documentation/Tutorials/HasModels/TankResponse.cc

Go to the documentation of this file.

#include <tls/VTankResponse.h>
#include <tls/TankResponseFactory.h>
#include <tls/TankResponseUtilities.h>

#include <boost/program_options.hpp>
#include <boost/shared_ptr.hpp>

#include <utl/Reader.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>

#include <iostream>
#include <cmath>
#include <string>

#include <TFile.h>
#include <TH1D.h>

namespace po = boost::program_options;
using namespace tls;
using namespace utl;
using namespace std;

string gXMLConfigurationFileName;
vector<double> gTheta;
int gMuonMin, gMuonMax;

// The following function is just to set the global variables and is at the end.
void processOptions(int argc, char* argv[]);

int main(int argc, char* argv[]) {

  processOptions(argc, argv);

  //TH1::AddDirectory(kFALSE); // this breaks everything in USC!
  TFile f("TankResponse.root","RECREATE");

  // Get the configuration. Which implementation is used is determined there.
  Reader reader(gXMLConfigurationFileName.c_str());
  Branch branch = reader.GetTopBranch().GetFirstChild();
  vector<Branch> branchList;
  while (branch) { branchList.push_back(branch); branch = branch.GetNextSibling(); }

  // calculate minimum and maximum signal to plot
  double minMu = gMuonMin-3*sqrt((double)gMuonMin);
  if (minMu<0) minMu = 0;
  double maxMu = gMuonMax+3*sqrt((double)gMuonMax);
  
  unsigned short nPoints = 400;

  for (unsigned short iz=0;iz<gTheta.size();++iz) {
    const double zenith = gTheta[iz];
    const double smin = minMu*TankMeanTrackLength(zenith*degree);
    const double smax = maxMu*TankMeanTrackLength(zenith*degree);
    
    for (int muons=gMuonMin;muons<=gMuonMax;++muons)
      for (unsigned short ib=0;ib<branchList.size();++ib)
    {
      Branch& branch = branchList[ib];
      
      printf("=============================================\n");
      printf("Using %s",branch.GetName().c_str());
      printf(" Theta = %.1f Deg, Muons = %i\n",zenith,muons);
      printf("=============================================\n");
  
      ostringstream name;
      name << branch.GetName() << "_" << int(zenith) << "_" << muons;

      try {
        VTankResponse& tankResponse = TankResponseFactory::GetInstance().GetTankResponse(branch);
    
        boost::shared_ptr<TH1D> h(new TH1D(name.str().c_str(),
          ";S / VEM;dp/dS / VEM^{-1}",nPoints,smin,smax));
    
        ostringstream name2;
        name2 << branch.GetName() << "_int_" << int(zenith) << "_" << muons;
        boost::shared_ptr<TH1D> h2(new TH1D(name2.str().c_str(),
          ";S_{min} / VEM;P(S > S_{min})",
          nPoints,smin,smax));
    
        for (int ix=1;ix<h->GetNbinsX()+1;++ix) {
          double s = h->GetXaxis()->GetBinCenter(ix);
          double p = tankResponse.PDF(s, zenith*degree, 1000*meter, muons);
          double p2 = tankResponse.CDF(s, zenith*degree, 1000*meter, muons);
          h->SetBinContent(ix,p);
          h2->SetBinContent(ix,p2);
          if (p<0 || p>1) {
            cerr << "WARNING: s = " << s << " VEM --> p = " << p << endl;
            if (p<0) p = -p;
          }
          if (p2<0 || p2>1) {
            cerr << "WARNING: s = " << s << " VEM --> p2 = " << p2 << endl;
            if (p2>1) p2 = 1.;
          }
        }
    
        cout << " Average signal              "
            << tankResponse.Mean(zenith*degree, 0, static_cast<int>(muons)) << "\n"
            << " PDF for S < 3.0 VEM "
            << 1.-tankResponse.CDF(3.0,zenith*degree, 1000*meter, static_cast<int>(muons))
            << endl;
    
        f.cd();
        h->Write();
        h2->Write();
      } catch (NonExistentComponentException& e) {
        ERROR(e.GetMessage());
      }
    } // loop over tank responses, no. of muons
  } // loop over zenith angles
  return 0;
}

void processOptions(int argc, char* argv[])
{
  ostringstream msg;
  msg << "Usage: " << argv[0] << " -c <xml_config_file> -z <zenith> [-z <zenith>...] --min <muons> --max <muons>\n"
      << "Example: '" << argv[0] << " -z 60 -z 80 --min 1 --max 10' will histogram the tank\n"
      << "response probability distributions at theta = 60 deg and theta = 80 deg, for 1 up\n"
      << "up to 10 muons. Everything is saved in the file testTankResponse.root.\n"
      << "\nAll options are";
  po::options_description desc(msg.str());
  desc.add_options()
    ("help,h",
     "write this message")
    ("config,c",
     po::value< string > (&gXMLConfigurationFileName),
     "path to configuration file. Required.")
    ("zenith,z",
     po::value< vector<double> >(&gTheta),
     "Zenith angle (in degrees). Required.")
    ("min",
     po::value<int>(&gMuonMin),
     "Minimal number of expected muons. Required.")
    ("max",
     po::value<int>(&gMuonMax),
     "Maximal number of expected muons. Required.")
    ;


  po::variables_map vm;
  try {
    po::store(po::parse_command_line(argc, argv, desc), vm);
    po::notify(vm);
  }
  catch (po::error& er) {
    cerr << "Command line error : " << er.what() << '\n'
         << desc << endl;
    exit(EXIT_FAILURE);
  }

  if (vm.count("help") || !vm.count("zenith") || !vm.count("min") || !vm.count("max")) {
    cerr << desc << endl;
    exit(EXIT_SUCCESS);
  }
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: