MuonProfileBuilder.cc

Go to the documentation of this file.

#include "MuonProfileBuilder.h"

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

#include <det/Detector.h>

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

#include <utl/ErrorLogger.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/Particle.h>
#include <utl/UTMPoint.h>

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include <iostream>

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

using namespace fwk;
using namespace evt;
using namespace utl;

using namespace std;

namespace MuonProfileBuilder {

  MuonProfileBuilder::MuonProfileBuilder() :
    fMuonEnergyCut(0),
    fWriteRootFile(false),
    fRootFileName("dummyFile.root")
  {
  }

  MuonProfileBuilder::~MuonProfileBuilder()
  {
  }

  VModule::ResultFlag
  MuonProfileBuilder::Init()
  {
    Branch topB =
      CentralConfig::GetInstance()->GetTopBranch("MuonProfileBuilder");

    topB.GetChild("MuonEnergyCut").GetData(fMuonEnergyCut);
    topB.GetChild("WriteRootFile").GetData(fWriteRootFile);
    topB.GetChild("RootFileName").GetData(fRootFileName);

    return eSuccess;
  }


  VModule::ResultFlag
  MuonProfileBuilder::Run(evt::Event& event)
  {
    ShowerSimData& simShower = event.GetSimShower();

    const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
    const double zenith = (-simShower.GetDirection()).GetTheta(localCS); 

    double coreHeight = utl::UTMPoint(simShower.GetPosition(), ReferenceEllipsoid::GetWGS84()).GetHeight();
    const CoordinateSystemPtr coreCS = simShower.GetLocalCoordinateSystem();

    // assuming flat atmosphere!
    // for curved atmosphere, use slant profile model and remove cos(zenith) below
    const atm::Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
    const atm::ProfileResult& DepthVsHeight = theAtm.EvaluateDepthVsHeight();

    // unit is gram/cm2 for the fit of the histogram
    const double histLowerBound = 0.;
    const double histUpperBound = 2000.;
    const double binSize = 5.;
    const int nBins = (histUpperBound-histLowerBound)/binSize;

    TFile* rootFile(0);
    if (fWriteRootFile) {
      rootFile = new TFile(fRootFileName.c_str(),"RECREATE");
    }

    TH1D survivingMuons("surviving", "surviving muons", nBins, histLowerBound, histUpperBound);
    TH1D decayedMuons("decayed", "decayed muons", nBins, histLowerBound, histUpperBound);
    TH1D allMuons("all", "all muons", nBins, histLowerBound, histUpperBound);

    int i = 0;
    int type = Particle::eUndefined;

    for (utl::ShowerParticleIterator particleIt = simShower.GroundParticlesBegin();
         particleIt != simShower.GroundParticlesEnd(); ++particleIt) {

      type = particleIt->GetType();

      if (! (type == Particle::eMuon || type == Particle::eAntiMuon
             || type == Particle::eDecayedMuon || type == Particle::eDecayedAntiMuon))

        continue;

      if (!particleIt->HasProductionPoint())
        continue;

      if (particleIt->GetKineticEnergy() < fMuonEnergyCut)
        continue;

      double prod_height = particleIt->GetProductionPoint().GetZ(coreCS) + coreHeight;
      double x_prod = DepthVsHeight.Y(prod_height)/cos(zenith);
      double weight = particleIt->GetWeight();

      if (type == Particle::eMuon || type == Particle::eAntiMuon)
        survivingMuons.Fill(x_prod/(gram/cm2), weight);
      else if (type == Particle::eDecayedMuon || type == Particle::eDecayedAntiMuon)
        decayedMuons.Fill(x_prod/(gram/cm2), weight);

      allMuons.Fill(x_prod/(gram/cm2), weight);

      i++;
      if (i % 100000 == 0) {
        ostringstream msg;
        msg << "read " << i << " muons";
        INFO(msg);
      }

    }


    double histMax = allMuons.GetMaximum();
    double maxBin = allMuons.GetMaximumBin();
    int firstFitBin = 0;
    int lastFitBin = 0;
    double binContent;

    for (int i = 0; i < nBins; ++i) {
      binContent = allMuons.GetBinContent(i);
      if (binContent > .7*histMax && firstFitBin == 0)
        firstFitBin = i;
      if (binContent < .7*histMax && firstFitBin > 0 && lastFitBin == 0 && i > maxBin)
        lastFitBin = i-1;
    }

    double firstValue = allMuons.GetBinCenter(firstFitBin);
    double lastValue = allMuons.GetBinCenter(lastFitBin);

    TF1 func("func", "[0] + [2] * (x - [1])**2 + [3] * (x - [1])**3", firstValue, lastValue);
    func.SetParameters(histMax, allMuons.GetBinCenter(maxBin), -1, 0);
    double XmaxMu = allMuons.Fit("func", "RS")->Parameter(1);

    {
      ostringstream msg;
      msg << "fitted muon production profile between "
          << firstValue << " g/cm2 and "
          << lastValue << " g/cm2. "
          << "XmaxMu = " << XmaxMu << " g/cm2. " << endl;
      INFO(msg);
    }

    simShower.SetXmaxMu(XmaxMu*gram/cm2);

    {
      ostringstream msg;
      msg << "overwriting longitudinal profiles for muons";
      WARNING(msg);
    }

    if (fWriteRootFile) {
      survivingMuons.Write();
      decayedMuons.Write();
      allMuons.Write();
      if (rootFile) {
        rootFile->Close();
        delete rootFile;
      }
    }

    TabulatedFunction profileSurviving;
    TabulatedFunction profileDecayed;
    TabulatedFunction profileAll;

    for (int i = 0; i < nBins; ++i) {
      profileSurviving.PushBack(survivingMuons.GetBinCenter(i)*gram/cm2, survivingMuons.GetBin(i));
      profileDecayed.PushBack(decayedMuons.GetBinCenter(i)*gram/cm2, decayedMuons.GetBin(i));
      profileAll.PushBack(allMuons.GetBinCenter(i)*gram/cm2, allMuons.GetBin(i));
    }

    if (!simShower.HasLongitudinalProfile(ShowerSimData::eMuon))
      simShower.MakeLongitudinalProfile(ShowerSimData::eMuon);

    simShower.GetLongitudinalProfile(ShowerSimData::eMuon) = profileSurviving;

    if (!simShower.HasLongitudinalProfile(ShowerSimData::eMuonProduction))
      simShower.MakeLongitudinalProfile(ShowerSimData::eMuonProduction);

    simShower.GetLongitudinalProfile(ShowerSimData::eMuonProduction) = profileAll;

    return eSuccess;
  }


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

}