DedxDistribution.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  DedxDistribution.cxx                                                                              
/// \brief This is a helper class for ParticleIDAlg that provides a tidy structure
///  in which to hold the dE/dx histograms and provides access methods to the data
///
///
///                                                                                                         
/// \version $Id: DedxDistribution.cxx,v 1.0 2013-10-04 10:00:00 absmith Exp $                             
/// \author  Alex Smith (smith@physics.umn.edu)                  
////////////////////////////////////////////////////////////////////////

#include "messagefacility/MessageLogger/MessageLogger.h"
#include "ShowerLID/DedxDistribution.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/WeightedHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Utilities/func/MathUtil.h"
#include "Utilities/func/EnvExpand.h"

#include "TGeoManager.h"
#include "TGeoNode.h"

#include <iostream>
#include "TH1.h"

namespace slid {

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   

  DedxDistribution::DedxDistribution()
    : fEmptyLongHist(0), fEmptyTransHist(0)
  {
  }


  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////               

  DedxDistribution::~DedxDistribution() {
      
    for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
      for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
        //
        // delete longitudinal dE/dx histograms
        //
        for (int iPlane = 0; iPlane < kNumLongitudinalPlane; iPlane++) {
          TH1F* h = fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane];
          if(h != fEmptyLongHist) delete h;
        }
        //
        // delete transverse dE/dx histograms
        //

        for (int iPlane = 0; iPlane < kNumTransversePlane; iPlane++) {
          TH1* h = fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iPlane];
          if(h != fEmptyTransHist) delete h;
        }
      }
    }

    delete fEmptyLongHist;
    delete fEmptyTransHist;
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   

  void DedxDistribution::Initialize(const fhicl::ParameterSet& pset,
                                    std::string fileNameFclTag) {

    std::string libPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
    std::string thistname = pset.get< std::string >(fileNameFclTag);
    /// FD dE/dx input histogram files
    std::string fileName = libPath + thistname;
    LoadHistFromFile(fileName);
    return;
  }

  void DedxDistribution::LoadHistFromFile(std::string filePath) {

    TFile* fDistFile = new TFile(filePath.c_str(), "READ");
    LoadDedxHistograms(fDistFile);
    fDistFile->Close();
    delete fDistFile;
  }

  ////////////////////////////////////////////////////////////////////////
  TH1F* DedxDistribution::TH1DToTH1F(TH1D* h, bool trans)
  {
    // A surprisingly large fraction (~1/2) of these histograms are empty. No
    // need to have multiple seperate copies of identical empty histograms.
    if(h->GetEntries() == 0){
      // If we already have an appropriate histogram cached, that's what we'll
      // return.
      TH1F* ret = 0;
      if(trans && fEmptyTransHist) ret = fEmptyTransHist;
      if(!trans && fEmptyLongHist) ret = fEmptyLongHist;

      // Check the binning of the cached histogram in case we ever try to have
      // variable binning within the longitudunal or transverse histogram sets.
      if(ret && 
         (h->GetNbinsX() != ret->GetNbinsX() ||
          h->GetXaxis()->GetXmin() != ret->GetXaxis()->GetXmin() ||
          h->GetXaxis()->GetXmax() != ret->GetXaxis()->GetXmax())){
        std::cout << "Error in DedxDistribution::TH1DToTH1F. "
                  << "Cached empty " << (trans ? "transverse" : "longitudinal")
                  << " histogram binning does not match loaded histogram. "
                  << ret->GetNbinsX() << " bins from "
                  << ret->GetXaxis()->GetXmin() << " to "
                  << ret->GetXaxis()->GetXmax() << " vs "
                  << h->GetNbinsX() << " bins from "
                  << h->GetXaxis()->GetXmin() << " to "
                  << h->GetXaxis()->GetXmax() << std::endl;
        abort();
      }

      if(ret){
        delete h;
        return ret;
      }
    }

    // Make a histogram the same size, but with floats. Don't copy the title
    // which is usually long and descriptive and takes up a surprising amount
    // of storage.
    TH1F* ret = new TH1F(h->GetName(),
                         "",
                         h->GetNbinsX(),
                         h->GetXaxis()->GetXmin(),
                         h->GetXaxis()->GetXmax());

    // Copy the bin contents across
    for(int i = 0; i < h->GetNbinsX() + 2; ++i){
      ret->SetBinContent(i, h->GetBinContent(i));
    }

    // We also need to set the number of entries which is checked by some code
    // that uses these histograms.
    ret->SetEntries(h->GetEntries());
    ret->SetDirectory(0);

    delete h;

    // If the histogram is empty and we got here, which means the cache is also
    // empty, set the cache entry to this histogram.
    if(ret->GetEntries() == 0){
      if(trans)
        fEmptyTransHist = ret;
      else
        fEmptyLongHist = ret;
    }

    return ret;
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////     

  bool DedxDistribution::LoadDedxHistograms(TFile* fDistFile) {

    //The detector is divided into 4 regions in XY, these regions extend in z. 
    //Region 1 is the top left quadrant in XY when looking at the front of the detector
    //The training histograms are then divided into 11 energy regions from 0 to 5 GeV.
    //In the training each shower is put into a bin based on its reco energy
    //The first and last energy bin are 0-0.25 and 4.75-5 GeV.  All bins inbetween are 0.5 GeV wide.
    //For each region and each energy there is a de/dx histogram for each plane from 0 to 200
    //the histograms being loaded below are for longitudinal dedx.
    //These training histograms are also broken up by particle type (electron, muon, gamma, proton, pi0, pion, neutron)

    //NOTE: Jianming is exploring using the shower length as a training variable.
    //In order to do this the variables below are initialized and then store the
    //shower length at each energy in each region.  This is done by finding the 
    //last plane in which there was energy deposited.
    //These variables are not currently being used in the production version of the algorithm
    for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
      for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
        fHtExpPlane[iXYRegion][iEnergyBin] = 0;
      }
    }
    //
    //  Load histograms for longitudinal dE/dx
    //
    for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
      for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
        for (int iPlane = 0; iPlane < kNumLongitudinalPlane; iPlane++) {
          char cht[200];
          sprintf(cht, "ht_%d_%d_%d", iXYRegion, iEnergyBin, iPlane);
          //convert TH1D to TH1F in order to save memory
          fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane] = TH1DToTH1F((TH1D*)fDistFile->Get(cht), false);

          if (fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetEntries() > 0 && fHtExpPlane[iXYRegion][iEnergyBin] == 0) {
            if (fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetBinContent(1) / (float) fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iPlane]->GetEntries() > 0.997) {
              fHtExpPlane[iXYRegion][iEnergyBin] = iPlane;
            }
          }
        }
      }
    }

    //
    //  Load histograms for transverse dE/dx
    //
    for (int iXYRegion = 0; iXYRegion < kNumXYRegion; iXYRegion++) {
      for (int iEnergyBin = 0; iEnergyBin < kNumEnergyBin; iEnergyBin++) {
        for (int iPlane = 0; iPlane < kNumTransversePlane; iPlane++) {
          char cht[200];
          sprintf(cht, "httrans_%d_%d_%d", iXYRegion, iEnergyBin, iPlane);
          //convert TH1D to TH1F in order to save memory.
          fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iPlane] = TH1DToTH1F((TH1D*)fDistFile->Get(cht), true);
        }
      }
    }
    return true;
  }



  ///  \brief Transverse plane hit dE/dx.

  TH1F* DedxDistribution::GetLongitudinalPlaneHitDedxHist(int iXYRegion, int iEnergyBin, int iLongitudinalPlane) {
    return fLongitudinalPlaneHitDedx[iXYRegion][iEnergyBin][iLongitudinalPlane];
  }
  ///  \brief Transverse plane hit dE/dx.

  TH1F* DedxDistribution::GetTransversePlaneCellDedxHist(int iXYRegion, int iEnergyBin, int iTransversePlane) {
    return fTransversePlaneCellDedx[iXYRegion][iEnergyBin][iTransversePlane];
  }





}// End of namespace slid