Documentation/Tutorials/HasModels/EMComponent.cc

Go to the documentation of this file.

#include <tls/EMComponent.h>
#include <tls/MuonProfileUtilities.h>

#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Reader.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/UTMPoint.h>
#include <utl/Transformation.h>
#include <utl/Vector.h>
#include <utl/Point.h>

#include <iostream>
#include <cstdlib>
#include <sstream>
#include <string>
#include <cmath>

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

#include <TROOT.h>
#include <TObjectTable.h>
#include <TFile.h>
#include <TH2.h>

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

string gXMLConfigurationFileName;
vector<int> gTheta;
vector<int> gLgRMu;

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

// This is to create a coordinate system (in Malargue) since the framework is not there.
CoordinateSystemPtr GetLocalCS();

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

  double xMax = 5.0; // km
  double yMax = 5.0; // km
  int xBins = 200;
  int yBins = 200;

  TFile* f = TFile::Open("EMComponent.root","RECREATE");

  // The geomagnetic field that will be used to change reference frames:
  CoordinateSystemPtr localCs = GetLocalCS();
  const Vector gMagneticField = Vector(1,
                                       kPiOnTwo + tls::kInclination,
                                       kPiOnTwo - tls::kDeclination,
                                       localCs,
                                       Vector::kSpherical);

  Reader reader(gXMLConfigurationFileName.c_str());
  Branch branch = reader.GetTopBranch().GetFirstChild();
  while(branch) {
    EMComponent emComponent(branch);

    for(unsigned int irmu = 0; irmu != gLgRMu.size(); ++irmu)
      for(unsigned int ith = 0; ith != gTheta.size(); ++ith)
    {
      const double theta  = gTheta[ith]*degree;
      const double rmu    = pow(10.,gLgRMu[irmu]);
      double phi = 90*degree;

      const Vector axis(1,theta,phi,localCs, Vector::kSpherical);
      const CoordinateSystemPtr showerCs =
        CoordinateSystem::RotationY(theta,CoordinateSystem::RotationZ(phi, localCs));
      const CoordinateSystemPtr geomagneticCs =
        CoordinateSystem::RotationZ(gMagneticField.GetPhi(showerCs) - kPiOnTwo, showerCs);

      ostringstream label, title;
      char buf[256]; sprintf(buf,"%s_%i_%i",branch.GetName().c_str(),
                                 int(round(log10(rmu))),
                                 int(round(theta/degree)) );
      label << buf;
      title << branch.GetName() << ";x / km;y / km;S_{em} / S_{#mu}";
      TH2D* histo = new TH2D(label.str().c_str(),title.str().c_str(),xBins,-xMax,xMax, yBins, -yMax, yMax);
      histo->SetDirectory(0);

      for(int ix = 1; ix < xBins+1; ++ix)
        for(int iy = 1; iy < yBins+1; ++iy)
      {
        const double x = histo->GetXaxis()->GetBinCenter(ix);
        const double y = histo->GetYaxis()->GetBinCenter(iy);

        Point p(x, y, 0,showerCs);

        Point p_ground = p - Vector(0, 0, p.GetZ(localCs)/axis.GetZ(localCs), showerCs);
        const double xp = p_ground.GetX(localCs);
        const double yp = p_ground.GetY(localCs);


        const double semsmu = emComponent.SignalRatio(xp*km, yp*km, rmu, theta, phi);
        histo->SetBinContent(ix,iy,semsmu);
      }

      histo->SetTitle("Em correction on the shower plane");
      histo->Write();
    }

    branch = branch.GetNextSibling();
  }

  f->Close();

  return 0;
}

void processOptions(int argc, char* argv[])
{
  ostringstream msg;
  msg << "Usage: " << argv[0] << " -c <xml_config_file> -z <zenith> [-z <zenith> ...] -m <lgRMu> [-m <lgRMu> ... ]\n"
      << "Example: '" << argv[0] << " -z 60 -z 70 -m 0' will histogram the em components for the zenith angles 60 deg and 70 deg and an energy estimator lgRMu = 0.\n"
      << "Everything is saved in the file EMComponent.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<int> >(&gTheta),
     "Zenith angle (in degrees). Required.")
    ("lgrmu,m",
     po::value< vector<int> >(&gLgRMu),
     "lg(energy estimator). 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("lgrmu") ) {
    cerr << desc << endl;
    exit(EXIT_SUCCESS);
  }
}

CoordinateSystemPtr GetLocalCS()
{
  ReferenceEllipsoid wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  const CoordinateSystemPtr ecef = wgs84.GetECEF();

  UTMPoint malargueOrigin(6060000*m, 440000*m, 1400*m, 19, 'H', wgs84);
  return AugerBaseCoordinateSystem::Create(malargueOrigin, ecef);
}

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