TestParticleCorrections_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    TestParticleCorrections_module.cc
// \brief   Make plots of threshold effect from PCHits files
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

//C++ and C:
#include <cmath>
#include <vector>

//Framework:
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"

//NovaSoft:
#include "Calibrator/Calibrator.h"
#include "CalibrationUtils/CalibUtil.h"
#include "CalibrationDataProducts/PCHit.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "RecoBase/RecoHit.h"

// ROOT
#include "TFile.h"
#include "TGraph.h"
#include "TH2.h"

namespace calib
{
  /// Make plots of threshold effect from PCHits files
  class TestParticleCorrections: public art::EDAnalyzer
  {
  public:
    explicit TestParticleCorrections(const fhicl::ParameterSet& pset);
    virtual ~TestParticleCorrections();

    virtual void beginJob() override;

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

    virtual void reconfigure(const fhicl::ParameterSet& pset);

  protected:
    TH1 *fE, *fMu, *fProt, *fPi;
    TH1 *fEPre, *fMuPre, *fProtPre, *fPiPre;

    TH1* fCell;

    TGraph *fGE, *fGMu, *fGProt, *fGPi;
  };

  //---------------------------------------------------------------------
  TestParticleCorrections::TestParticleCorrections(const fhicl::ParameterSet& pset)
    : EDAnalyzer(pset)
  {
    reconfigure(pset);
  }

  //---------------------------------------------------------------------
  TestParticleCorrections::~TestParticleCorrections()
  {
  }

  //---------------------------------------------------------------------
  void TestParticleCorrections::reconfigure(const fhicl::ParameterSet& pset)
  {
  }

  //---------------------------------------------------------------------
  void TestParticleCorrections::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;

    fE = tfs->make<TH1F>("e", "True e^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
    fMu = tfs->make<TH1F>("mu", "True #mu^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
    fProt = tfs->make<TH1F>("prot", "True p;Energy estimate (GeV);Slices", 400, 0, 4);
    fPi = tfs->make<TH1F>("pi", "True #pi^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);


    fEPre = tfs->make<TH1F>("e_pre", "True e^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
    fMuPre = tfs->make<TH1F>("mu_pre", "True #mu^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
    fProtPre = tfs->make<TH1F>("prot_pre", "True p;Energy estimate (GeV);Slices", 400, 0, 4);
    fPiPre = tfs->make<TH1F>("pi_pre", "True #pi^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);

    fE->SetLineColor(kRed);
    fMu->SetLineColor(kBlue);
    fProt->SetLineColor(kBlack);
    fPi->SetLineColor(kMagenta);

    fEPre->SetLineColor(kRed);
    fMuPre->SetLineColor(kBlue);
    fProtPre->SetLineColor(kBlack);
    fPiPre->SetLineColor(kMagenta);

    fCell = tfs->make<TH1F>("cell", ";Cell reco/true", 100, 0, 2);

    TFile* f = new TFile("pcorr_curves.root");
    fGE = (TGraph*)f->Get("e");
    fGMu = (TGraph*)f->Get("mu");
    fGProt = (TGraph*)f->Get("p");
    fGPi = (TGraph*)f->Get("pi");
  }

  //---------------------------------------------------------------------
  double corr(TGraph* g, double thresh)
  {
    for(int i = 0; i < g->GetN(); ++i){
      double x, y;
      g->GetPoint(i, x, y);

      if(x > thresh) return 100/y;
    }

    return 1; // Uh oh
  }

  //---------------------------------------------------------------------
  void TestParticleCorrections::analyze(const art::Event& evt)
  {
    const double hack = 4095/2170.; // why?

    art::ServiceHandle<cheat::BackTracker> bt;
    art::ServiceHandle<calib::Calibrator> cal;

    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel("slicer", slices);

    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      const art::Ptr<rb::Cluster> slice(slices, sliceIdx);
      if(slice->IsNoise()) continue;

      for(unsigned int iCell = 0; iCell < slice->NCell(); ++iCell){
        const double reco = slice->RecoHit( iCell).GeV()*hack;
        const std::vector<sim::FLSHit> flss = bt->HitToFLSHit(slice->Cell(iCell));
        double truth = 0;
        for(auto f: flss) truth += f.GetEdep();

        fCell->Fill(reco/truth);
      }

      const double est = slice->TotalGeV();

      const std::vector<const sim::Particle*> parts = bt->HitsToParticle(slice->AllCells());

      if(parts.empty()) continue;
      const sim::Particle* part = parts[0];

      rb::CellHit dummy;
      dummy.SetPlane(slice->MinPlane());
      dummy.SetCell(slice->MinCell(geo::kX));
      dummy.SetADC(0, 25);
      dummy.SetPE(cal->GetPE(&dummy));

      dummy.SetMC();

      const rb::RecoHit rhit = cal->MakeRecoHit(dummy, slice->MeanX()); // symmetric

      // Crash when not calibrated. Always should be in MC
      const double thresh = rhit.GeV()*1000; // MeV
      //      std::cout << thresh << std::endl;

      // From gdml, correct to total density of soup
      const double dead = (.85*.859+.15*1.49)/(.85*.859);

      switch(part->PdgCode()){
      case -11:
      case +11:
        fE->Fill(est*corr(fGE, thresh)*dead*hack);
        fEPre->Fill(est*hack);
        break;
      case -13:
      case +13:
        fMu->Fill(est*corr(fGMu, thresh)*dead*hack);
        fMuPre->Fill(est*hack);
        break;
      case 2212:
        fProt->Fill(est*corr(fGProt, thresh)*dead*hack);
        fProtPre->Fill(est*hack);
        break;
      case -211:
      case +211:
        fPi->Fill(est*corr(fGPi, thresh)*dead*hack);
        fPiPre->Fill(est*hack);
        break;
      }
    }
  }

  DEFINE_ART_MODULE(calib::TestParticleCorrections)

} //end namespace