AbsCalib_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       AbsCalib
// Module Type: analyzer
// File:        AbsCalib_module.cc
//
// Generated at Fri May 10 14:12:00 2013 by Abigail Waldron using artmod
// from cetpkgsupport v1_02_00.
////////////////////////////////////////////////////////////////////////

#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/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include <TH2D.h>
#include <TH1D.h>
#include <TF1.h>
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include <iostream>

namespace calib{
  class AbsCalib;
}

namespace calib{
  class AbsCalib : public art::EDAnalyzer {
  public:
    explicit AbsCalib(fhicl::ParameterSet const & p);
    virtual ~AbsCalib();
    
    void analyze(art::Event const & e) override;

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


  private:

    double fitMEU(TH1D*& proj, double& sigma, bool isTruth=false); 
    
    // Declare member data here.
    double fWindowStart;
    double fWindowLength;

    double fMEU, fMEU_minus, fMEU_plus, fMEU_true;
    double fMEU_error, fMEU_error_minus, fMEU_error_plus, fMEU_error_true;

    TH2D* fdEdx;
    TH2D* fdEdx_minus;
    TH2D* fdEdx_plus;
    TH2D* fdEdx_true;
    TH1D* fWindow;
    TH1D* fWindow_minus;
    TH1D* fWindow_plus;
    TH1D* fWindow_true;
  };
  
  
  AbsCalib::AbsCalib(fhicl::ParameterSet const & p)
    : EDAnalyzer(p)
  {
    fWindowStart = 100.;
    fWindowLength = 100.;
  }
  
  AbsCalib::~AbsCalib()
  {
  }
  
  void AbsCalib::analyze(art::Event const & e)
  {
  }
  
  void AbsCalib::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fdEdx = tfs->make<TH2D>("dEdx","; ; ;",
                            50,0.,500.,200,0.,200.);
    fdEdx_minus = tfs->make<TH2D>("dEdx_minus","; ; ;",
                            50,0.,500.,200,0.,200.);
    fdEdx_plus = tfs->make<TH2D>("dEdx_plus","; ; ;",
                            50,0.,500.,200,0.,200.);
    fdEdx_true = tfs->make<TH2D>("dEdx_true","; ; ;",
                                 50,0.,500.,200,0.,5.);
    fWindow = tfs->make<TH1D>("window","window",
                              200,0.,200.);
    fWindow_minus = tfs->make<TH1D>("window_minus","window_minus",
                              200,0.,200.);
    fWindow_plus = tfs->make<TH1D>("window_plus","window_plus",
                              200,0.,200.);
    fWindow_true = tfs->make<TH1D>("window_true","window_true",
                                   200,0.,5.);
  }
  
  void AbsCalib::endJob()
  {
    // For now calculate MEU per job, maybe later move to run
    fMEU = this->fitMEU(fWindow,fMEU_error);
    fMEU_minus = this->fitMEU(fWindow_minus,fMEU_error_minus);
    fMEU_plus = this->fitMEU(fWindow_plus,fMEU_error_plus);
    fMEU_true = this->fitMEU(fWindow_true,fMEU_error_true,true);


    double sys_err = fabs(fMEU_plus - fMEU);
    sys_err = (fabs(fMEU - fMEU_minus) > sys_err ) ? fabs(fMEU - fMEU_minus) : sys_err;

    std::cout << "MEU = " << fMEU << " +/-" << fMEU_error << "(stat) PECorr/cm"
              << " +/-" << sys_err << "(sys) PECorr/cm"
              << std::endl;


    std::cout << "True MEU = " << fMEU_true << " +/-" << fMEU_error_true << "(stat) MeV/cm"
              << std::endl;
 

  }
  
  double AbsCalib::fitMEU(TH1D*& proj, double& sigma, bool isTruth){
    
    Double_t parameters[3];
    parameters[0] = proj->GetBinContent(proj->GetMaximumBin());
    parameters[1] = proj->GetBinCenter(proj->GetMaximumBin());
    parameters[2] = proj->GetRMS();


    // only care about finding peak position, use gaussian truncated about mode  
    Double_t bound_size = (isTruth) ? 0.2 : 0.5; // narrower for truth, *not* MC reco... don't mess this up.

    Double_t lowerBound = std::max(0.,parameters[1] - bound_size*proj->GetRMS());
    Double_t upperBound = parameters[1] + bound_size*proj->GetRMS();


    TF1* fitfunction = new TF1("fitfunction","[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
    fitfunction->SetParameters(parameters);


    // fitting                                                          
    proj->Fit(fitfunction,"RQ");


    // retrieve values
    sigma = fitfunction->GetParError(1);
    return fitfunction->GetParameter(1);
  }


  void AbsCalib::beginRun(art::Run const & r)
  {
  }
  
  void AbsCalib::endRun(art::Run const & r)
  {
    art::Handle<TH2D> tmp, tmpMinus, tmpPlus, tmpTrue;
    r.getByLabel("muondedx","dEdx",tmp);
    r.getByLabel("muondedx","dEdxMinus",tmpMinus);
    r.getByLabel("muondedx","dEdxPlus",tmpPlus);
    r.getByLabel("muondedx","dEdxTrue",tmpTrue);
    
    for( int i=0; i<fdEdx->GetNbinsX(); ++i ){
      for( int j=0; j<fdEdx->GetNbinsY(); ++j ){
        double thing = fdEdx->GetBinContent(i+1,j+1);
        double thingMinus = fdEdx_minus->GetBinContent(i+1,j+1);
        double thingPlus = fdEdx_plus->GetBinContent(i+1,j+1);
        double thingTrue = fdEdx_true->GetBinContent(i+1,j+1);
        thing += tmp->GetBinContent(i+1,j+1);   
        thingMinus += tmpMinus->GetBinContent(i+1,j+1); 
        thingPlus += tmpPlus->GetBinContent(i+1,j+1);   
        thingTrue += tmpTrue->GetBinContent(i+1,j+1);   
        fdEdx->SetBinContent(i+1,j+1,thing);
        fdEdx_minus->SetBinContent(i+1,j+1,thingMinus);
        fdEdx_plus->SetBinContent(i+1,j+1,thingPlus);
        fdEdx_true->SetBinContent(i+1,j+1,thingTrue);   
      }
    }

    // project into track window
    int startBin = fdEdx->GetXaxis()->FindBin(fWindowStart);
    int endBin = 1 + fdEdx->GetXaxis()->FindBin(fWindowStart+fWindowLength);
    fWindow->Add((TH1D*)tmp->ProjectionY("window",startBin,endBin));
    fWindow_minus->Add((TH1D*)tmpMinus->ProjectionY("window_minus",startBin,endBin));
    fWindow_plus->Add((TH1D*)tmpPlus->ProjectionY("window_plus",startBin,endBin));
    fWindow_true->Add((TH1D*)tmpTrue->ProjectionY("window_true",startBin,endBin));
  }
}// end namspace
  
DEFINE_ART_MODULE(calib::AbsCalib)