SdSimMuonNumberFitter.cc

Go to the documentation of this file.

#include "SdSimMuonNumberFitter.h"

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/UTMPoint.h>
#include <utl/Reader.h>
#include <utl/Particle.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>

#include <sdet/MuonProfile.h>
#include <utl/HASUtilities.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/Header.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <sevt/StationRecData.h>

#include <sstream>
#include <fstream>
#include <iomanip>
#include <iostream>

#include <TFile.h>
#include <TMath.h>
#include <TH2D.h>
#include <TMinuit.h>
#include <TCanvas.h>
#include <TVirtualPad.h>

#include <boost/tuple/tuple.hpp>

using namespace SdSimMuonNumberFitterNS;
using namespace fwk;
using namespace utl;
using namespace HASUtilities;
using namespace sdet;
using namespace evt;
using namespace sevt;
using namespace std;


typedef boost::tuple<double, double> Tuple2D;

//ofstream gDebugLog("debug.log");

static int kHASRecoStatus = 10;

MuonProfile* fgMuonProfile = nullptr;

// square of the perpendicular distance from the axis
inline
double
RPerp2(const Vector& axis, const Vector& station)
{
  const double scal = axis*station;
  return station*station - scal*scal;
}


// perpendicular distance from the axis
inline
double
RPerp(const Vector& axis, const Vector& station)
{
  return sqrt(RPerp2(axis, station));
}


template<typename G>
inline
string
ToString(const G& g, const CoordinateSystemPtr cs)
{
  ostringstream os;
  os << '('
     << g.GetX(cs)/meter << ", "
     << g.GetY(cs)/meter << ", "
     << g.GetZ(cs)/meter << ") [m]";
  return os.str();
}


inline
CoordinateSystemPtr
GetShowerCoordinateSystem(const double theta, const double phi,
  const Point& core)
{
  return CoordinateSystem::RotationY(theta,
                                     CoordinateSystem::RotationZ(phi,
                                                                 LocalCoordinateSystem::Create(core)));
}


// project point into core plane in the direction of the shower origin
inline
Tuple2D
ProjectIntoCorePlane(const Point& pos, const Point& showerOrigin, const CoordinateSystemPtr& coreCS)
{
  const double dz = pos.GetZ(coreCS);
  const Vector stationToOrigin = showerOrigin-pos;
  const Point corePlanePos = stationToOrigin/stationToOrigin.GetZ(coreCS)*dz + pos;
  return Tuple2D(corePlanePos.GetX(coreCS),corePlanePos.GetY(coreCS));
}


inline
Tuple2D
ProjectIntoCorePlane2(const Point& pos, const Point& showerOrigin, const CoordinateSystemPtr& coreCS)
{
  const Vector axis = showerOrigin - Point(0,0,0, coreCS);
  const double dz = pos.GetZ(coreCS);
  const Point pos2 = pos - axis/axis.GetZ(coreCS)*dz;
  return Tuple2D(pos2.GetX(coreCS), pos2.GetY(coreCS));
}


// get local shower inclination at position "pos"
inline
double
LocalCosTheta(const Point& pos, const Point& showerOrigin)
{
  const Vector stationToOrigin = showerOrigin - pos;
  return stationToOrigin.GetCosTheta(LocalCoordinateSystem::Create(pos));
}


void
Likelihood::AddStation(const int id, const Point& pos, const Point& showerOrigin, const CoordinateSystemPtr& coreCS, const unsigned int nmuon)
{
  // projection into flat muon density plane along shower axis
  // axis points upwards to the shower "source", we safely assume that axis is normalised
  Tuple2D sPosProjXY = ProjectIntoCorePlane(pos, showerOrigin, coreCS);
  const double sTheta = acos(LocalCosTheta(pos, showerOrigin));
  fStationList.push_back(InternalStationData(id, sPosProjXY.get<0>(), sPosProjXY.get<1>(), sTheta,nmuon));
}


void
Likelihood::MinuitFnc(int& nPar, double* const grad, double& value, double* const par, const int flag)
{
  // p[0] N19
  // p[1] x_core
  // p[2] y_core

  if (par[0] <= 0) {
    value = DBL_MAX;
    return;
  }

  double logLikelihood = 0;
  for (InternalStationList::const_iterator sIt = fStationList.begin();
       sIt != fStationList.end(); ++sIt) {
    const double& x = sIt->fX;
    const double& y = sIt->fY;
    //const double& th = sIt->fTheta;
    const int& nmu = sIt->fNMuon;
    double nmuExpected = 0;
    try {
      nmuExpected = par[0]
        * fgMuonProfile->NMuon(x - par[1], y - par[2], fTheta, fPhi, 10*EeV);
    } catch (OutOfBoundException& e) {
      nmuExpected = 1e-9;
    }

    double p = TMath::Poisson(nmu, nmuExpected);
    if (p > 0)
      logLikelihood += log(p);
    else
      logLikelihood += -DBL_MAX_EXP;
  }

  if (logLikelihood < -DBL_MAX) {
    ERROR("logLikelihood = -infinity!");
    value = DBL_MAX;
    return;
  }

  if (logLikelihood != logLikelihood) {
    ERROR("Likelihood = nan!");
    value = DBL_MAX;
    return;
  }

  value = -logLikelihood;
}


static Likelihood gLikelihood;


SdSimMuonNumberFitter::SdSimMuonNumberFitter()
{ }


SdSimMuonNumberFitter::~SdSimMuonNumberFitter()
{
  delete fgMuonProfile;
}


VModule::ResultFlag
SdSimMuonNumberFitter::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();

  Branch topB = cc->GetTopBranch("SdSimMuonNumberFitter");

  topB.GetChild("Verbosity").GetData(fVerbosity);

  // muon profiles are only defined up to about 4000 m distance
  topB.GetChild("MaxDistanceZeroSignalStations").GetData(fSilentRadius);

  string muonProfile;
  topB.GetChild("MuonProfile").GetData(muonProfile);
  Branch muonProfileBranch = cc->GetTopBranch("MuonProfile").GetChild(muonProfile.c_str());
  fgMuonProfile = new sdet::MuonProfile(muonProfileBranch);

  topB.GetChild("EnergyConstant").GetData(fEnergyConstant);

  topB.GetChild("Gamma").GetData(fGamma);

  fFixedCore = bool(topB.GetChild("FixedCore"));

  if (fVerbosity >= ErrorLogger::eInfo) {
    ostringstream msg;
    msg << "\nParameters:\n"
        << "  Verbosity                       : " << fVerbosity << "\n"
        << "  MaxDistanceZeroSignalStations   : " << fSilentRadius << " m\n"
        << "  MuonProfile                     : " << muonProfileBranch.GetName()
        << " Theta limits [" << fgMuonProfile->GetThetaMin()/degree << " Deg, "
        << fgMuonProfile->GetThetaMax()/degree << " Deg]" << "\n"
        << "  EnergyConstant(EnergyConversion): " << fEnergyConstant << " eV\n"
        << "  Gamma         (EnergyConversion): " << fGamma << "\n"
        << "  Fixed core                      : " << (fFixedCore ? "yes\n" : "no\n");
    INFO(msg);
  }

  return eSuccess;
}


VModule::ResultFlag
SdSimMuonNumberFitter::Run(evt::Event & event)
{
  // nothing to do if there is no SimShower or no SEvent
  if (!event.HasSimShower())
    return eContinueLoop;
  if (!event.HasSEvent())
    return eContinueLoop;

  sevt::SEvent& sEvent = event.GetSEvent();
  const evt::ShowerSimData& showerSimData = event.GetSimShower();
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  if (!event.HasRecShower())
    event.MakeRecShower();
  if (!event.GetRecShower().HasSRecShower())
    event.GetRecShower().MakeSRecShower();
  ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();

  // Get true shower core, axis, core time
  const Point& simCore = showerSimData.GetPosition();
  const TimeStamp& simTime = showerSimData.GetTimeStamp();
  const Vector& showerAxis = -showerSimData.GetDirection();
  fSimCoreCS = showerSimData.GetLocalCoordinateSystem();
  
  const double energyMC = showerSimData.GetEnergy();
  const double thetaMC = showerSimData.GetZenith();
  const double phiMC = showerSimData.GetAzimuth();

  // Build up station lists
  vector<StationData> candidateList;
  vector<StationData> silentList;
  for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {
    const sdet::Station& dStation = sDetector.GetStation(*sIt);
    const sevt::StationSimData& sSim = sIt->GetSimData();

    StationData sd;
    sd.fId = sIt->GetId();
    sd.fPos = dStation.GetPosition();
    sd.fNMuon = 0;
    sd.fTime = TimeStamp(0,0);

    bool noParticleData = true;
    for (sevt::StationSimData::ConstParticleIterator it = sSim.ParticlesBegin();
         it != sSim.ParticlesEnd(); ++it) {
      noParticleData = false;
      if (abs(it->GetType()) == Particle::eMuon)
        ++sd.fNMuon;
      const TimeStamp pTime = simTime + it->GetTime();
      if (!sd.fTime)
        sd.fTime = pTime;
      if (pTime < sd.fTime)
        sd.fTime = pTime;
    }
    if (noParticleData && sIt->IsCandidate() && sIt->HasRecData()) {
       const sevt::StationRecData& sRec = sIt->GetRecData();
      // if we don't have the particle information, assume that the muon number is encoded in the tank signal value
      sd.fNMuon = (unsigned int)(sRec.GetTotalSignal());
      sd.fTime = sRec.GetSignalStartTime();
    }

    if (sd.fNMuon > 0)
      candidateList.push_back(sd);
    else if (RPerp(showerAxis, sd.fPos - simCore) < fSilentRadius)
      silentList.push_back(sd);
  }

  if (fVerbosity > ErrorLogger::eDebug) {
    for (unsigned int i = 0; i < candidateList.size(); ++i) {
      const StationData& sd = candidateList[i];
      cerr << "Station " << sd.fId << ": " << sd.fNMuon << endl;
    }
  }

  if(fVerbosity > ErrorLogger::eInfo) {
    ostringstream msg;
    msg << "\tSEvent Id = " << sEvent.GetHeader().GetId() << "\n"
        << "\tThetaMC " << showerSimData.GetZenith()/deg << " [deg]" << "\n"
        << "\t# total number of stations     = " << sEvent.GetNumberOfStations() << "\n"
        << "\t# number of candidate stations = " << candidateList.size() << "\n"
        << "\t# number of silent    stations = " << silentList.size() << endl;
    INFO(msg);
  }

  // not worth it
  if (candidateList.size() < 3) {
    INFO("Not enough stations.");
    showersrec.SetS1000(0, 0);
    return eSuccess;
  }

  const Point barycenter = FindMuonDensityBaryCenter(candidateList);

  // once a barycenter, create: local coordinate system, calculate shower origin, barytime
  const CoordinateSystemPtr baryCS = LocalCoordinateSystem::Create(barycenter);
  Point showerOrigin = barycenter + showerAxis*AverageCurvatureRadius(thetaMC);
  const TimeStamp coreTime = simTime -
    TimeInterval(showerAxis*(barycenter - simCore)/kSpeedOfLight);

  for (vector<StationData>::const_iterator sIt = candidateList.begin();
       sIt != candidateList.end(); ++sIt) {
    sevt::Station& station = sEvent.GetStation(sIt->fId);
    if (!station.HasRecData())
      station.MakeRecData();
    StationRecData& sRec = station.GetRecData();

    const double signal = sIt->fNMuon;

    const TimeStamp sTime = sIt->fTime + (coreTime - simTime);

    sRec.SetTotalSignal(signal, sqrt(signal));
    sRec.SetSignalStartTime(sTime);
    sRec.SetT50(0,0);
  }

  // warn, if event is out of range
  const bool notInMuonProfileRange  =
    thetaMC < fgMuonProfile->GetThetaMin() || fgMuonProfile->GetThetaMax() < thetaMC;

  if (notInMuonProfileRange) {
    if (fVerbosity > ErrorLogger::eInfo) {
      ostringstream msg;
      msg << "Theta = " << thetaMC/degree << " Deg was not in valid range of\n ";
      if (notInMuonProfileRange)
        msg << " MuonProfile";
      msg << "\n";
      INFO(msg);
    }
  }

  // energy reco
  Point core;
  Vector coreErr;
  double n19 = 0;
  double n19Err = 0;
  double logLikelihoodMin = 0;
  gLikelihood.Clear();
  for (vector<StationData>::const_iterator sIt=candidateList.begin();
       sIt != candidateList.end(); ++sIt) {
    gLikelihood.AddStation(sIt->fId, sIt->fPos, showerOrigin, baryCS, sIt->fNMuon);
  }
  for (vector<StationData>::const_iterator sIt = silentList.begin();
       sIt != silentList.end(); ++sIt) {
    gLikelihood.AddStation(sIt->fId, sIt->fPos, showerOrigin, baryCS, 0);
  }
  gLikelihood.SetTheta(thetaMC);
  gLikelihood.SetPhi(phiMC);

  n19 = 1;
  bool success = EnergyFit(baryCS, core, coreErr, n19, n19Err, logLikelihoodMin, true);
  if (!success) {
    ERROR("Failed N19 estimate fit (fixed core to barycenter)");
    showersrec.SetS1000(0, 0);
    return eSuccess;
  }
  const double n19estimate = n19;

  if (200 <= fVerbosity) { // draw loglikelihood debug stuff
    char buf[1024];
    int phii = int(phiMC/degree);
    if (phii < 0)
      phii += 360;
    sprintf(buf, "lh_map_%09i_%i_%i.root", sEvent.GetHeader().GetId(), int(thetaMC/degree), phii);
    TFile f(buf, "RECREATE");
    TCanvas* const c = new TCanvas;
    c->Divide(1, 2);
    c->cd(1);
    DrawMuonDensity(gPad, n19estimate, thetaMC, phiMC, "est");
    c->cd(2);
    DrawLogLikelihood(gPad, n19estimate, "est");
    c->Write();
    f.Close();
    delete c;
  }

  success = EnergyFit(baryCS, core, coreErr, n19, n19Err, logLikelihoodMin, false);
  if (!success) {
    showersrec.SetS1000(n19estimate, 0);
    ERROR("Failed N19 + core fit");

    if (100 <= fVerbosity && fVerbosity < 200) { // draw loglikelihood debug stuff
      char buf[1024];
      int phii = int(phiMC/degree);
      if (phii < 0)
        phii += 360;
      sprintf(buf, "lh_map_%09i_%i_%i.root", sEvent.GetHeader().GetId(), int(thetaMC/degree), phii);
      TFile f(buf, "RECREATE");
      TCanvas* const c = new TCanvas;
      c->Divide(1, 2);
      c->cd(1);
      DrawMuonDensity(gPad, n19estimate, thetaMC, phiMC, "est");
      c->cd(2);
      DrawLogLikelihood(gPad, n19estimate, "est");
      c->Write();
      f.Close();
      delete c;
    }

    return eSuccess;
  }
  showerOrigin = core + showerAxis*AverageCurvatureRadius(thetaMC);
  
  // Convert N19 to energy
  const double Energy = fEnergyConstant * pow(n19,fGamma);
  const double EnergyError = n19Err * fEnergyConstant*fGamma*pow(n19, fGamma-1);

  const double diffToMC = Energy/energyMC - 1;
  const double diffToMCSigma = (Energy - energyMC) / EnergyError;

  if (fVerbosity > ErrorLogger::eInfo) {
    ostringstream msg;
    msg << "Result:\n"
      << "\tEnergy - XML Calibration   : " << Energy/EeV << " +/- " << EnergyError/EeV
      << " [EeV]\n"
      << "\tMonte-Carlo Energy         : " << energyMC/EeV << " [EeV]"
      << " ( " << setprecision(3) << diffToMC*100.
      << " % | " << setprecision(3) << diffToMCSigma << " std.dev. )" << endl;
    INFO(msg);
  }

  // fill SRecShower
  if (!event.HasRecShower())
    event.MakeRecShower();
  if (!event.GetRecShower().HasSRecShower())
    event.GetRecShower().MakeSRecShower();

  const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);
  const CoordinateSystemPtr showerPlaneCS = GetShowerCoordinateSystem(thetaMC, phiMC, core);

  double chi2 = 0;
  unsigned int ndof = 0;
  for (vector<StationData>::const_iterator sIt = candidateList.begin();
       sIt != candidateList.end(); ++sIt) {
    const Point& sPos = sIt->fPos;

    const Tuple2D sPosProjXY = ProjectIntoCorePlane(sPos, showerOrigin, coreCS);
    
    double muons = 0;
    try {
      muons = n19
        * fgMuonProfile->NMuon(sPosProjXY.get<0>(),
                               sPosProjXY.get<1>(), thetaMC, phiMC, 10*EeV);
    } catch (OutOfBoundException& e) {
      continue;
    }

    const double signalPred = muons;
    sevt::Station& station = sEvent.GetStation(sIt->fId);
    StationRecData& sRec = station.GetRecData();

    const double signal = sIt->fNMuon;
    const double sigma = sqrt(signalPred);
    const double res = (signal - signalPred) / sigma;

    const double rho = RPerp(showerAxis, sPos - core);
    const double drho = 0; // TODO calculate this, see LDFFinder

    chi2 += (res*res) / signalPred; // sigma^2 = mean for poisson statistics
    ++ndof;

    sRec.SetTotalSignal(signal, sigma);
    sRec.SetAzimuthShowerPlane(sPos.GetPhi(showerPlaneCS));
    sRec.SetLDFResidual((signal - signalPred)/sigma);
    sRec.SetSPDistance(rho, drho);
  }

  showersrec.SetBarycenter(barycenter);
  showersrec.SetCorePosition(core);
  showersrec.SetCoreError(coreErr);
  showersrec.SetCoreTime(coreTime, TimeInterval(0));
  showersrec.SetAxis(showerAxis);
  showersrec.SetCurvature(0, 0);
  showersrec.SetEnergy(Energy, EnergyError);
  showersrec.SetS1000(n19, n19Err);
  showersrec.SetLDFChi2(chi2, ndof);
  showersrec.SetLDFLikelihood(logLikelihoodMin);
  // LDFRecStatus variable > 10 for HAS
  showersrec.SetLDFRecStage(ShowerSRecData::kLDF + kHASRecoStatus);

  return eSuccess;
}


VModule::ResultFlag
SdSimMuonNumberFitter::Finish()
{
  return eSuccess;
}


utl::Point
SdSimMuonNumberFitter::FindMuonDensityBaryCenter(const vector<StationData>& stationList)
  const
{
  const det::Detector& detector = det::Detector::GetInstance();

  // calculate barycenter and bary-time
  const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();

  // the absolute points in time and space
  const Point siteOrigin(0,0,0, siteCS);

  // weighted sums
  Vector barySum(0,0,0, siteCS);
  double weightSum = 0;
  int nStations = 0;
  for (vector<StationData>::const_iterator it = stationList.begin();
       it != stationList.end(); ++it) {
    const StationData& sd = *it;
    const double weight = sqrt(double(sd.fNMuon));
    weightSum += weight;
    barySum += weight * (sd.fPos - siteOrigin);
    ++nStations;
  }

  if (!weightSum) {
    ERROR("sum of weights = 0");
    return Point();
  }

  barySum /= weightSum;
  Point barycenter = siteOrigin + barySum;

  ostringstream msg;
  msg << "\nbarycenter = " << ToString(barycenter, fSimCoreCS) << " (showerGroundCS)" << "\n";
  INFO(msg);

  return barycenter;
}


void
SdSimMuonNumberFitter::N19FitFnc(int& nPar, double* const grad, double& value, double* const par, const int flag)
{
  gLikelihood.MinuitFnc(nPar, grad, value, par, flag);
}


bool
SdSimMuonNumberFitter::EnergyFit(const CoordinateSystemPtr& localCS,
                                 Point& core, Vector& coreErr, double& n19, double& n19Err,
                                 double& logLikelihoodMin, const bool fixedCore)
  const
{
  TMinuit minuit(3);

  bool badFit = false;

  int errFlag = 0;
  double argList[10];
  argList[0] = -1;

  minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
  minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
  minuit.SetPrintLevel(-1);

  minuit.SetFCN(SdSimMuonNumberFitter::N19FitFnc);
  minuit.mnparm(0, "N19"  , n19,         0.2, 0.0, 0.0, errFlag);
  minuit.mnparm(1, "xCore", 0.0, 100.0*meter, 0.0, 0.0, errFlag);
  minuit.mnparm(2, "yCore", 0.0, 100.0*meter, 0.0, 0.0, errFlag);

  //Change in FCN of one-standard-deviation error
  argList[0] = 0.5;
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

  /* Minimisation strategy: Try MIGRAD minimisation first - which provides also uncertainty estimates.
     If MIGRAD fails, the more robust SIMPLEX algorithm is used to find the minimum, which does NOT
     provide uncertainty estimates. In that case, we run MINOS to the find the uncertainties.
     If SIMPLEX fails too, mark the reconstruction as bad.
  */
  if (fixedCore) {
    minuit.FixParameter(1);
    minuit.FixParameter(2);
  }
  argList[0] = 500; argList[1] = 1;
  minuit.mnexcm("MIGRAD", argList, 2, errFlag);
  if (errFlag) {
    ostringstream msg;
    msg << "minuit MIGRAD error: " << errFlag;
    ERROR(msg);
    minuit.mnexcm("SIMPLEX", argList , 2, errFlag);
    if (errFlag) {
      ostringstream msg;
      msg << "minuit SIMPLEX error: " << errFlag;
      ERROR(msg);
      badFit = true;
    } else {
      minuit.mnexcm("MINOS", argList , 1, errFlag);
      if (errFlag) {
        ostringstream msg;
        msg << "minuit MINOS error: " << errFlag;
        ERROR(msg);
        badFit = true;
      }
    }
  }

  // Get results
  double xShift = 0;
  double xShiftErr = 0;
  double yShift = 0;
  double yShiftErr = 0;
  minuit.GetParameter(0, n19, n19Err);
  minuit.GetParameter(1, xShift, xShiftErr);
  minuit.GetParameter(2, yShift, yShiftErr);

  if (n19 <= 0 || n19Err <= 0 || TMath::IsNaN(n19) || TMath::IsNaN(n19Err)) {
    ERROR("Fit failed: N19 or N19Err are bogus.");
    badFit = true;
  }

  if (!fFixedCore) {
    core = utl::Point(xShift, yShift, 0, localCS);
    coreErr = utl::Vector(xShiftErr, yShiftErr, 0, localCS);
  }

  logLikelihoodMin = -minuit.fAmin; // minus sign here, because we minimise -LL in Minuit

  return !badFit;
}


void
SdSimMuonNumberFitter::DrawLogLikelihood(TVirtualPad* const pad, const double n19, const char* const comment)
  const
{
  int nPar = 3;
  double par[3];
  double value = 0;
  par[0] = n19;
  int flag = 0;
  double* grad = nullptr;

  INFO("Draw 1");
  TH2D* const h1 = new TH2D(TString("lh_map1") += comment, "", 100, -25e3, 25e3, 100, -25e3, 25e3);
  for (int ix = 1; ix < h1->GetNbinsX()+1; ++ix)
    for(int iy = 1; iy < h1->GetNbinsY()+1; ++iy) {
      par[1] = h1->GetXaxis()->GetBinCenter(ix);
      par[2] = h1->GetYaxis()->GetBinCenter(iy);
      gLikelihood.MinuitFnc(nPar, grad, value, par, flag);
      h1->SetBinContent(ix,iy,value);
    }

  INFO("Draw 2");
  TH2D* const h2 = new TH2D(TString("lh_map2") += comment, "", 100, -5e3, 5e3, 100, -5e3, 5e3);
  for (int ix = 1; ix < h2->GetNbinsX()+1; ++ix)
    for(int iy = 1; iy < h2->GetNbinsY()+1; ++iy) {
      par[1] = h2->GetXaxis()->GetBinCenter(ix);
      par[2] = h2->GetYaxis()->GetBinCenter(iy);
      gLikelihood.MinuitFnc(nPar, grad, value, par, flag);
      h2->SetBinContent(ix,iy,value);
    }
  pad->Divide(2);
  pad->cd(1);
  h1->Draw("colz");
  pad->cd(2);
  h2->Draw("colz");
}


void
SdSimMuonNumberFitter::DrawMuonDensity(TVirtualPad* const pad, const double n19, const double theta, const double phi, const char* const comment)
  const
{
  INFO("Draw 1");
  TH2D* const h1 = new TH2D(TString("nmu_map1") += comment, "", 50, -25e3, 25e3, 50, -25e3, 25e3);
  for (int ix = 1; ix < h1->GetNbinsX()+1; ++ix)
    for(int iy = 1; iy < h1->GetNbinsY()+1; ++iy) {
      const double x = h1->GetXaxis()->GetBinCenter(ix);
      const double y = h1->GetYaxis()->GetBinCenter(iy);
      double muons = 0.0;
      try {
        muons = n19
          *fgMuonProfile->NMuon(x,y, theta, phi, 10*EeV);
      } catch (OutOfBoundException& e)
      { }
    h1->SetBinContent(ix, iy, muons);
  }

  INFO("Draw 2");
  TH2D* const h2 = new TH2D(TString("nmu_map2") += comment, "", 50, -5e3, 5e3, 50, -5e3, 5e3);
  for (int ix = 1; ix < h2->GetNbinsX()+1; ++ix)
    for(int iy = 1; iy < h2->GetNbinsY()+1; ++iy) {
      const double x = h2->GetXaxis()->GetBinCenter(ix);
      const double y = h2->GetYaxis()->GetBinCenter(iy);
      double muons = 0.0;
      try {
        muons = n19
          * fgMuonProfile->NMuon(x,y, theta, phi, 10*EeV);
      } catch (OutOfBoundException& e)
      { }
    h2->SetBinContent(ix, iy, muons);
  }
  pad->Divide(2);
  pad->cd(1);
  h1->Draw("colz");
  pad->cd(2);
  h2->Draw("colz");
}