BetheBlochFit_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file   BetheBlochFit_module.cc
// \brief  Bethe-Bloch fit to muondEdx output
// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"

#include <iostream>

#include "TF1.h"
#include "TFile.h"
#include "TGeoManager.h"
#include "TGraph.h"
#include "TH2.h"
#include "TProfile.h"

#include "CalibrationUtils/CalibUtil.h"
#include "Calibration/func/BetheBloch.h"
#include "Calibration/func/EnergyLossVsDistance.h"

#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"


namespace calib
{
  class BetheBlochFit: public art::EDAnalyzer
  {
  public:
    explicit BetheBlochFit(const fhicl::ParameterSet& pset);
    virtual ~BetheBlochFit();

    void analyze(const art::Event& evt) override;

    void beginJob() override;
    void endJob() override;
    void beginRun(const art::Run& run) override;
    void endRun(const art::Run& run) override;

  protected:
    TH1* ProfileMPV(TH2* h) const;

    double fWindowStart;
    double fWindowLength;

    TH2* fdEdx;
    TH2* fdEdx_true;
  };

  //......................................................................
  BetheBlochFit::BetheBlochFit(const fhicl::ParameterSet& pset)
    : EDAnalyzer(pset)
  {
    // TODO: get these from pset
    fWindowStart = 100.;
    fWindowLength = 100.;
  }

  //......................................................................
  BetheBlochFit::~BetheBlochFit()
  {
  }

  //......................................................................
  void BetheBlochFit::analyze(const art::Event& evt)
  {
  }

  //......................................................................
  void BetheBlochFit::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fdEdx = tfs->make<TH2D>("dEdx", ";Distance to track end (cm);PECorr / cm",
                            50, 0, 500, 100, 0, 200);
    fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm",
                                 50, 0, 500, 100, 0, 5);
  }

  //......................................................................
  void BetheBlochFit::beginRun(const art::Run& run)
  {
  }

  //......................................................................
  void BetheBlochFit::endRun(const art::Run& run)
  {
    art::Handle<TH2D> dE, dEtrue;
    run.getByLabel("muondedx", "dEdx", dE);
    run.getByLabel("muondedx", "dEdxTrue", dEtrue);

    fdEdx->Add(&*dE);
    fdEdx_true->Add(&*dEtrue);
  }

  //......................................................................
  TH1* BetheBlochFit::ProfileMPV(TH2* h) const
  {
    TH1* ret = new TH1F("", "", h->GetNbinsX(),
                        h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
    for(int x = 0; x < h->GetNbinsX()+2; ++x){
      double maxVal = -1e10;
      double maxy = 0;
      for(int y = 1; y < h->GetNbinsY()+1; ++y){
        const double val = h->GetBinContent(x, y);
        if(val > maxVal){
          maxVal = val;
          maxy = h->GetYaxis()->GetBinCenter(y);
        }
      }
      TH1* proj = h->ProjectionY("", x, x);
      proj->Fit("gaus", "0Q", "", maxy-0.2, maxy+0.2);
      TF1* func = proj->GetFunction("gaus");
      if(func){
        ret->SetBinContent(x, func->GetParameter(1));
        ret->SetBinError(x, func->GetParError(1));
      }
    }

    return ret;
  }

  //......................................................................
  void BetheBlochFit::endJob()
  {
    art::ServiceHandle<geo::Geometry> geom;

    // This differs from soup by about 0.2%
    // Representative of the whole detector
    TGeoVolume* block = geom->ROOTGeoManager()->FindVolumeFast("vBlock");
    const double avgDensity = 1e3*block->Weight()/block->Capacity();

    //    const double avgDensity =
    //      geom->ROOTGeoManager()->GetMaterial("Soup")->GetDensity();

    const double scintDensity =
      geom->ROOTGeoManager()->GetMaterial("Scintillator")->GetDensity();

    std::cout << "Detector average (soup) density is " << avgDensity << std::endl;
    std::cout << "Scintillator density is " << scintDensity << std::endl;

    TProfile* prof = fdEdx->ProfileX();

    IBetheBloch* bb_ana = new BetheBlochAnalytic(kPolyethylene);
    IBetheBloch* bb_tab = new BetheBlochTables(kPolyethylene);

    EnergyLossVsDistance el_ana(bb_ana, scintDensity, avgDensity);
    EnergyLossVsDistance el_tab(bb_tab, scintDensity, avgDensity);

    TF1* f_ana = el_ana.GetTF1();
    prof->Fit(f_ana, "", "", 50, 500);
    TF1* f_tab = el_tab.GetTF1();
    prof->Fit(f_tab, "", "", 50, 500);

    // Average seen by muon travelling along z-axis
    const double Lz = (geom->Plane(0)->Cell(4)->HalfD()+geom->Plane(1)->Cell(4)->HalfD())/2;

    const double pePerMeV = f_ana->GetParameter(0);

    std::cout << "1 MeV = " << pePerMeV << " PE -> 1 PE = " << 1/pePerMeV << " MeV" << std::endl;

    const double MeVPerCmPerMIP = bb_ana->MIP()*scintDensity;

    std::cout << "1 MIP = " << MeVPerCmPerMIP*Lz << " MeV [MeVScale]" << std::endl;

    const double pePerMIPPerCm = pePerMeV*MeVPerCmPerMIP;

    std::cout << "1 MIP = " << pePerMIPPerCm*Lz << " PE -> 1 PE = " << 1/pePerMIPPerCm/Lz << " MIP [MIPScale]" << std::endl;

    // Subvert the file service
    art::ServiceHandle<art::TFileService> tfs;
    tfs->file().cd("bb");
    f_ana->Draw("c");
    f_ana->GetHistogram()->Write("fit_ana");
    f_tab->Draw("c");
    f_tab->GetHistogram()->Write("fit_tab");

    el_ana.GetEnergyLossGraph()->Write("el_ana");
    el_tab.GetEnergyLossGraph()->Write("el_tab");

    bb_tab->GetdEdxGraph(scintDensity)->Write("dEdx_curve");
    bb_ana->GetdEdxGraph(scintDensity)->Write("dEdx_curve_ana");
    bb_ana->GetMPVGraph(scintDensity, Lz)->Write("mpv_curve_ana");
    bb_ana->GetMPVGraph(scintDensity, Lz*2)->Write("mpv_curve_ana2");
    bb_ana->GetMPVGraph(scintDensity, Lz/2)->Write("mpv_curve_ana3");

    ProfileMPV(fdEdx)->Write("prof_mpv");
  }

  DEFINE_ART_MODULE(BetheBlochFit)

} // namespace