LIDTrainingDedx.cxx

Go to the documentation of this file.

#include "TChain.h"
#include "TFile.h"
#include "TH1D.h"

#include <iostream>
#include <fstream>
#include <string>
#include "ShowerLID/train/LIDTrainingDedx.h"
#include <math.h>

namespace slidt{

  /////////////////////////////////////////////////////////////////////////////////////////

  void ActivateBranches()
  {

    // if(fActivated)
    //   return;
          
    fChain->SetBranchAddress("showerTrueNuCCNC",        &showerTrueNuCCNC);
    fChain->SetBranchAddress("showerTrueNuMode",        &showerTrueNuMode);
    fChain->SetBranchAddress("showerTrueNuPdg",         &showerTrueNuPdg);
    fChain->SetBranchAddress("showerNPlane",            &showerNPlane);
    fChain->SetBranchAddress("showerStart[3]",          showerStart);
    fChain->SetBranchAddress("showerStop[3]",           showerStop);
    fChain->SetBranchAddress("showerEnergy",            &showerEnergy);
    fChain->SetBranchAddress("showerIsFid",             &showerIsFid);
    fChain->SetBranchAddress("showerPlaneDedx[200]",    showerPlaneDedx);
    fChain->SetBranchAddress("showerTransCellDedx[20]", showerTransCellDedx);
    fChain->SetBranchAddress("showerTrueDang",          &showerTrueDang);
    fChain->SetBranchAddress("showerTruePdg",           &showerTruePdg);
    fChain->SetBranchAddress("showerTrueP4[4]",         showerTrueP4);
    fChain->SetBranchAddress("showerDir[3]",            showerDir);
    
    fActivated = true;

  }// end of ActivateSetBranchAddresses


  /////////////////////////////////////////////////////////////////////////////////////////
  
  void SaveHists(const int inPdg, 
                 const std::string outDir, 
                 const int*  inNuMode)
  {
    std::string partname, mode = "";
    
    switch(inPdg){
    case 11:
      partname = "electron";
      break;
    case 22:
      partname = "gamma";
      break;
    case 13:
      partname = "muon";
      break;
    case 111:
      partname = "pi0";
      break;
    case 211:
      partname = "pion";
      break;
    case 2212:
      partname = "proton";
      break;
    case 2112:
      partname = "neutron";
      break;
    default:
      std::cerr<<"The specified PDG particle is not included in PID.\n";
      std::abort();
    }

    if( inNuMode != NULL ){
      switch(*inNuMode){
      case 0:
        mode = "_qe";
        break;
      case 1:
        mode = "_res";
        break;
      case 2:
        mode = "_dis";
        break;
      case 3:
        mode = "_coh";
        break;
      default:
        std::cerr<<"Invalid neutrino interaction mode. Valid choices are 0,1,2,3.\n";
        std::abort();
      }
    }

    std::string outFileName = outDir + partname + mode + "_dedxhists.root";
    std::cout<<"output file name is "<<outFileName<<std::endl;
    TFile fout(outFileName.c_str(), "RECREATE");

    for( int idet=0; idet < kNumXYRegion; ++idet){
      for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
        for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
          hlongdedx[idet][iebin][iplane]->Write();
        }
      }
    }
    for( int idet=0; idet < kNumXYRegion; ++idet){
      for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
        for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
          htransdedx[idet][iebin][iplane]->Write();
        }
      }
    }
    fout.Close();
  }// End of SaveHists


  /////////////////////////////////////////////////////////////////////////////////////////
  
  void Reset()
  {
    for( int idet=0; idet < kNumXYRegion; ++idet){
      for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
        for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
          htransdedx[idet][iebin][iplane]->Reset();
        }
      }
    }
    for( int idet=0; idet < kNumXYRegion; ++idet){
      for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
        for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
          hlongdedx[idet][iebin][iplane]->Reset();
        }
      }
    }
    
    fChain->Reset();
  }// Reset all histograms

  /////////////////////////////////////////////////////////////////////////////////////////

  void InitializeHists()
  {
    if( fInitialized )
      return;

    for( int idet=0; idet < kNumXYRegion; ++idet){
      for( int iebin=0; iebin < kNumEnergyBin; ++iebin){

        for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
          char hname[100], htitle[200];
          sprintf(hname,"hlongdedx_det%d_ebin%d_plane%d", idet, iebin, iplane);
          sprintf(htitle,"Longitudinal dE/dx: Det Region %d, Energy Bin %d, Plane %d;dE/dX (GeV/cm);Events", idet, iebin, iplane);
          hlongdedx[idet][iebin][iplane] =  new TH1D(hname, htitle, 50, 0, 0.050);        
        }// end loop over longitudinal planes

        for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
          char hname[100], htitle[200];
          sprintf(hname,"htransdedx_det%d_ebin%d_plane%d", idet, iebin, iplane);
          sprintf(htitle,"Transverse dE/dx: Det Region %d, Energy Bin %d, Plane %d;dE/dX (GeV/cm);Events", idet, iebin, iplane);
          float hmax = 0.01;
          if( iplane==0 )
            hmax = 0.01;
          if( iplane==1 )
            hmax = 0.01;
          if( iplane==2 )
            hmax = 0.002;
          if( iplane>=3 && iplane<=5 )
            hmax = 0.001;
          if( iplane>=6 && iplane<=9 )
            hmax = 0.0005;
          if( iplane>=10 && iplane<=14 )
            hmax = 0.00010;
          if( iplane>=15 )
            hmax = 0.00005;
          htransdedx[idet][iebin][iplane] =  new TH1D(hname, htitle, 500, 0, hmax);       

        }// end loop over transverse planes

      }// end loop over energy bins
    }// end loop over detregions
    fInitialized = true; 
    return;
  }// end of InitializeHists

  /////////////////////////////////////////////////////////////////////////////////////////

  void FillHists( const char*   inFnames, 
                  const char*   inTreename,
                  const char*   outDir,
                  const int     inPdg,
                  const float   inDangLower,
                  const float   inDangUpper,
                  const int*    inNuPdg,
                  const int*    inNuMode,
                  const int*    inNuCC,
                  const double* eTrueFrac)
  {

    fChain = new TChain(inTreename);
    fChain->Add(inFnames);

    ActivateBranches();
    InitializeHists();
    
    int nentries = fChain->GetEntries();
    std::cout<<"Number of entries "<<nentries<<std::endl;

    for(int e = 0; e < nentries; ++e){
      
      fChain->GetEntry(e);
      
      // skip to next event if any of the required cuts are not
      // passed

      if( abs(showerTruePdg) != inPdg          ||
          showerTrueDang     <= inDangLower    ||
          showerTrueDang     >=  inDangUpper    ||
          showerIsFid        != 1               )
        continue;
      
      // the following three cuts don't apply to all particles.
      // see what's been provided.
      if( inNuPdg != NULL ){
        if( showerTrueNuPdg != *inNuPdg )
          continue;
      }
      
      if( inNuMode != NULL ){
        if( showerTrueNuMode != *inNuMode )
          continue;
      }
        
      if( inNuCC != NULL ){
        if( showerTrueNuCCNC != *inNuCC )
          continue;
      }

      if( eTrueFrac != NULL ){
        if( showerEnergy/showerTrueP4[3] <= *eTrueFrac )
          continue;
      }
      
      // check which energy bin the shower falls in
      static const double* eBins = kEnergyBins;
      int nEBins = kNumEnergyBin;
      int showerEBin = -1;
      for(int ibin = 0; ibin <= nEBins; ++ibin){
        if(showerEnergy > eBins[ibin] && showerEnergy < eBins[ibin+1]){
          showerEBin = ibin;
          if( showerEBin == 2)
          break;
        }
      }// end loop over energy bins
      
      if( showerEBin == -1)
        continue;

      // figure out which of the four detector regions  we are in
      int detRegion = -1;
      if( (showerStart[0]+showerStop[0])/2.0>0&&(showerStart[1]+showerStop[1])/2.0>0 )
        detRegion = 0;
      if( (showerStart[0]+showerStop[0])/2.0>0&&(showerStart[1]+showerStop[1])/2.0<0 )
        detRegion = 1;
      if( (showerStart[0]+showerStop[0])/2.0<0&&(showerStart[1]+showerStop[1])/2.0>0 )
        detRegion = 2;
      if( (showerStart[0]+showerStop[0])/2.0<0&&(showerStart[1]+showerStop[1])/2.0<0 )
        detRegion = 3;
      
      int nplanes = std::min(kNumLongitudinalPlane, showerNPlane);

      for( int iplane =0; iplane < nplanes; ++iplane)
        hlongdedx[detRegion][showerEBin][iplane]->Fill( showerPlaneDedx[int(iplane*fabs(showerDir[2]))] );
      
      for( int iplane =0; iplane < kNumTransversePlane; ++iplane)
        htransdedx[detRegion][showerEBin][iplane]->Fill( showerTransCellDedx[iplane] );
      
    }// end loop over entries
    
    std::cout<<"the output directory is "<<outDir<<std::endl;
    std::string outDirStr(outDir);
    
    SaveHists( inPdg, outDirStr, inNuMode);
    Reset();
  }// end of FillHists
  
}// end of namespace