Documentation/Tutorials/HasModels/MuonProfile.cc

Go to the documentation of this file.

#include <det/Detector.h>

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

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

#include <string>
#include <sstream>
#include <vector>

using namespace tls;
using namespace utl;
using namespace std;

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

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

namespace po = boost::program_options;

string gXMLConfigurationFileName;
Vector gMagneticField;
vector<double> gTheta;
vector<double> gPhi;

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

  const double xMax = 3000*m;
  const double yMax = 3000*m;
  const int xBins = 200;
  const int yBins = 200;

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

  // Get the configurations for the muon maps:
  Reader reader(gXMLConfigurationFileName.c_str());
  Branch topBranch = reader.GetTopBranch();
  Branch b = topBranch.GetFirstChild();
  vector<MuonProfile> muonProfiles;
  vector<string> muonProfileNames;
  while (b)
  {
    muonProfileNames.push_back(b.GetName());
    muonProfiles.push_back(MuonProfile(b));
    b = b.GetNextSibling();
  }
    
  // The geomagnetic field that will be used to change reference frames:
  CoordinateSystemPtr localCs = GetLocalCS();
  gMagneticField = Vector(1,
                          kPiOnTwo + tls::kInclination,
                          kPiOnTwo - tls::kDeclination,
                          localCs,
                          Vector::kSpherical);

  // Loop for all values of zenith and azimuth, filling histograms according to the muon maps
  for(unsigned int i = 0; i != gTheta.size(); ++i)
    for(unsigned int j = 0; j != gPhi.size(); ++j)
  {
    const double theta  = gTheta[i]*degree;
    const double phi    = gPhi[j]*degree;
    const double energy = 10.*EeV;

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

    for (unsigned int iType=0;iType<muonProfiles.size();++iType)
    {
      const string& name = muonProfileNames[iType];
      MuonProfile& mp = muonProfiles[iType];
      
      ostringstream label;
      label << name << "NMu" << "_" << theta/deg << "_" << phi/deg;
      boost::shared_ptr<TH2D> histoNMu(new TH2D(label.str().c_str(),name.c_str(),
        xBins,-xMax,xMax, yBins, -yMax, yMax));

      label.str("");
      label << name << "ThMu" << "_" << theta/deg << "_" << phi/deg;
      boost::shared_ptr<TH2D> histoThMu(new TH2D(label.str().c_str(),name.c_str(),
        xBins,-xMax,xMax, yBins, -yMax, yMax));

      label.str("");
      label << name << "EnMu" << "_" << theta/deg << "_" << phi/deg;
      boost::shared_ptr<TH2D> histoEnMu(new TH2D(label.str().c_str(),name.c_str(),
        xBins,-xMax,xMax, yBins, -yMax, yMax));

      cout << "Processing " << name << " " << theta/deg << " " << phi/deg << endl;
      
      for(int i = 1; i <= xBins+1; ++i) for(int j = 1; j <= yBins+1; ++j)
      {
        const double x = histoNMu->GetXaxis()->GetBinCenter(i);
        const double y = histoNMu->GetYaxis()->GetBinCenter(j);
        Point p(x, y, 0, geomagneticCs);
        Point p_ground = p - Vector(0, 0, p.GetZ(localCs)/axis.GetZ(localCs), geomagneticCs);
        const double xp = p_ground.GetX(localCs);
        const double yp = p_ground.GetY(localCs);
       
        histoNMu->Fill(x, y, mp.NMuon(xp, yp, theta, phi, energy));
        histoThMu->Fill(x, y, mp.ThetaMuon(xp, yp, theta, phi, energy));
        histoEnMu->Fill(x, y, mp.EnergyMuon(xp, yp, theta, phi, energy));
      }

      f.cd();
      histoNMu->Write();
      histoThMu->Write();
      histoEnMu->Write();
    } // loop over types
  } // loop over theta, phi combinatins

  f.cd();
  f.Close();

  return 0;
}

void processOptions(int argc, char* argv[])
{
  ostringstream msg;
  msg << "Usage: " << argv[0] << " -c <xml_config_file> -z <zenith> [-z <zenith>...] -a <azimuth> [-a <azimuth>...]\n"
      << "Example: '" << argv[0] << " -c MuonProfile.xml -z 60 -z 70 -a 0 -a 90' will histogram the maps for the following directions: (60 deg,0), (60 deg,90 deg), (70 deg,0) and (79 deg,90 deg).\n"
      << "Everything is saved in the file MuonProfile.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.")
    ("azimuth,a",
     po::value< vector<double> >(&gPhi),
     "Azimuth angle (in degrees). 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("azimuth")) {
    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: