eventRates.h

Go to the documentation of this file.

//////////////////////////////////////////////////////////
// Created by L. Loiacono
// Modified by L. Fields
//////////////////////////////////////////////////////////

#ifndef eventRates_h
#define eventRates_h

// C++
#include <string>
#include <vector>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>

//ROOT
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TRandom3.h>

// Dk2nu
//#ifndef __CINT__
#include "dk2nu/tree/dk2nu.h"
//#endif // ifndef __ CINT__

class eventRates {
public :

   TChain          *fChain;   //!pointer to the analyzed TTree or TChain
   //TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain
   
   Double_t        fTotalPOT; //total pot used for all files
   std::string     ffilename; //filename for saving histograms

   std::string detectorname;
   double detx;     // detector location
   double dety;
   double detz;

   TRandom3 *rand3;

   // Declaration of leaf types
 //LBNEDataNtp_t   *data;
   Int_t           run;
   Int_t           evtno;
   Int_t           protonN;
   Float_t         beamHWidth;
   Float_t         beamVWidth;
   Float_t         beamX;
   Float_t         beamY;
   Float_t         protonX;
   Float_t         protonY;
   Float_t         protonZ;
   Float_t         protonPx;
   Float_t         protonPy;
   Float_t         protonPz;
   Float_t         nuTarZ;
   Float_t         hornCurrent;
   Float_t         Ndxdz;
   Float_t         Ndydz;
   Float_t         Npz;
   Float_t         Nenergy;
   Float_t         NdxdzNear[5];
   Float_t         NdydzNear[5];
   Float_t         NenergyN[5];
   Double_t        NWtNear[5];
   Float_t         NdxdzFar[3];
   Float_t         NdydzFar[3];
   Float_t         NenergyF[3];
   Double_t        NWtFar[3];
   Int_t           Norig;
   Int_t           Ndecay;
   Int_t           Ntype;
   Float_t         Vx;
   Float_t         Vy;
   Float_t         Vz;
   Float_t         pdPx;
   Float_t         pdPy;
   Float_t         pdPz;
   Float_t         ppdxdz;
   Float_t         ppdydz;
   Float_t         pppz;
   Float_t         ppenergy;
   Float_t         ppmedium;
   Int_t           ptype;
   Int_t           ptrkid;
   Float_t         ppvx;
   Float_t         ppvy;
   Float_t         ppvz;
   Float_t         muparpx;
   Float_t         muparpy;
   Float_t         muparpz;
   Float_t         mupare;
   Float_t         Necm;
   Double_t        Nimpwt;
   Float_t         xpoint;
   Float_t         ypoint;
   Float_t         zpoint;
   Float_t         tvx;
   Float_t         tvy;
   Float_t         tvz;
   Float_t         tpx;
   Float_t         tpy;
   Float_t         tpz;
   Int_t           tptype;
   Int_t           tgen;
 //map<int,TrackPoint_t> fTrkPtMap;

   bool isDk2nu = false;

   // List of branches
   bsim::Dk2Nu* dk2nu;

   TBranch        *b_data_run;   //!
   TBranch        *b_data_evtno;   //!
   TBranch        *b_data_protonN;   //!
   TBranch        *b_data_beamHWidth;   //!
   TBranch        *b_data_beamVWidth;   //!
   TBranch        *b_data_beamX;   //!
   TBranch        *b_data_beamY;   //!
   TBranch        *b_data_protonX;   //!
   TBranch        *b_data_protonY;   //!
   TBranch        *b_data_protonZ;   //!
   TBranch        *b_data_protonPx;   //!
   TBranch        *b_data_protonPy;   //!
   TBranch        *b_data_protonPz;   //!
   TBranch        *b_data_nuTarZ;   //!
   TBranch        *b_data_hornCurrent;   //!
   TBranch        *b_data_Ndxdz;   //!
   TBranch        *b_data_Ndydz;   //!
   TBranch        *b_data_Npz;   //!
   TBranch        *b_data_Nenergy;   //!
   TBranch        *b_data_NdxdzNear;   //!
   TBranch        *b_data_NdydzNear;   //!
   TBranch        *b_data_NenergyN;   //!
   TBranch        *b_data_NWtNear;   //!
   TBranch        *b_data_NdxdzFar;   //!
   TBranch        *b_data_NdydzFar;   //!
   TBranch        *b_data_NenergyF;   //!
   TBranch        *b_data_NWtFar;   //!
   TBranch        *b_data_Norig;   //!
   TBranch        *b_data_Ndecay;   //!
   TBranch        *b_data_Ntype;   //!
   TBranch        *b_data_Vx;   //!
   TBranch        *b_data_Vy;   //!
   TBranch        *b_data_Vz;   //!
   TBranch        *b_data_pdPx;   //!
   TBranch        *b_data_pdPy;   //!
   TBranch        *b_data_pdPz;   //!
   TBranch        *b_data_ppdxdz;   //!
   TBranch        *b_data_ppdydz;   //!
   TBranch        *b_data_pppz;   //!
   TBranch        *b_data_ppenergy;   //!
   TBranch        *b_data_ppmedium;   //!
   TBranch        *b_data_ptype;   //!
   TBranch        *b_data_ptrkid;   //!
   TBranch        *b_data_ppvx;   //!
   TBranch        *b_data_ppvy;   //!
   TBranch        *b_data_ppvz;   //!
   TBranch        *b_data_muparpx;   //!
   TBranch        *b_data_muparpy;   //!
   TBranch        *b_data_muparpz;   //!
   TBranch        *b_data_mupare;   //!
   TBranch        *b_data_Necm;   //!
   TBranch        *b_data_Nimpwt;   //!
   TBranch        *b_data_xpoint;   //!
   TBranch        *b_data_ypoint;   //!
   TBranch        *b_data_zpoint;   //!
   TBranch        *b_data_tvx;   //!
   TBranch        *b_data_tvy;   //!
   TBranch        *b_data_tvz;   //!
   TBranch        *b_data_tpx;   //!
   TBranch        *b_data_tpy;   //!
   TBranch        *b_data_tpz;   //!
   TBranch        *b_data_tptype;   //!
   TBranch        *b_data_tgen;   //!

   eventRates(std::string input_user, std::string version, std::string macro, 
              std::string current, std::string location, std::string physics_list, 
              int n_files, int start_index, double pot_per_file, bool on_grid,
              bool dk2nu);

   virtual ~eventRates();
   virtual Int_t    Cut(Long64_t entry);
   virtual Int_t    GetEntry(Long64_t entry);
   virtual Long64_t LoadTree(Long64_t entry);
   virtual void     Init(TTree *tree);
   virtual void     Loop();
   virtual Bool_t   Notify();
   virtual void     Show(Long64_t entry = -1);

   std::string GetPOTAsString(const double dpot);
   void SetTitles(TH1* h, 
                  const std::string &xtitle = "", 
                  const std::string &ytitle = "");
   
   double GetWeight(const std::vector<double> xdet,
                            double& nu_wght, 
                            double& nu_energy);

   double GetXSec(  double nu_type, 
                    double nu_energy,
                    std::string current);

   void ReadXSecsFromFiles(  );

   double GetOscillatedNeutrinoType(double ntype, double E);


 private:

    std::ifstream fdat_file[6];
    int fnbins;
    int fnlines;
    double f_e_arr[1500][6][2]; // energy bins; neutrino type; CC/NC
    double f_xsec_arr[1500][6][2]; // energy bins; neutrino type; CC/NC
   
};

#endif

#ifdef eventRates_cxx

eventRates::eventRates(std::string input_user, std::string version, std::string macro, std::string current, std::string location, std::string physics_list, int n_files,int start_index, double potperfile, bool on_grid,bool dk2nu)
{
  // simulation = G4PBeam (default) or Fluka
  // macro = Nominal, etc
  // location = LBNEFD (default), LBNEND, etc 
  // physics_list = QGSP_BERT (default), QGSP, FTFP_BERT, etc
  // n_files = 250 (default)
  // potperfile = 100000 (default)

  isDk2nu = dk2nu;

  // Read the cross-sections from files
  ReadXSecsFromFiles();

  // Otherwise Set detector coordinates based on requested detector location
  if(location=="LBNEND") {
    detx = 0.0;
    dety = 0.0;
    detz = 57400.0;
  }
  
  else if(location=="LBNEFD") {
    detx = 0.0;
    dety = 0.0;
    detz = 129700000.0;
  }

  detectorname = location;

  std::string input_flux_dir =  "/pnfs/dune/scratch/users/"+input_user+"/fluxfiles/g4lbne/"+version+"/"+physics_list+"/"+macro+"/"+current+"/flux/";
  if(input_user=="beam")
    input_flux_dir = "/dune/data/beam/fluxfiles/g4lbne/"+version+"/"+physics_list+"/"+macro+"/"+current+"/flux/";
  if(on_grid) {
    input_flux_dir = getenv("_CONDOR_SCRATCH_DIR");
    input_flux_dir = input_flux_dir + "/";
  }

  std::vector<std::string> fFileVec;

  
  for(int i = start_index; i< start_index+n_files; i++) {
    
    // convert index into zero-padded string
    std::ostringstream ss_5;
    std::ostringstream ss_3;
    ss_5 << std::setw( 5 ) << std::setfill( '0' ) << i;
    ss_3 << std::setw( 3 ) << std::setfill( '0' ) << i;
    std::string index_string_5 = ss_5.str();
    std::string index_string_3 = ss_3.str();
    
    std::string flux_file = input_flux_dir + "g4lbne_"+version+"_"+physics_list+"_"+macro+"_"+current+"_"+index_string_3+".root";

    if(isDk2nu)
      flux_file = input_flux_dir + "g4lbne_"+version+"_"+physics_list+"_"+macro+"_"+current+"_dk2nu_"+index_string_5+".root";

    // check that the file exists and is a valid root file
    TFile f(flux_file.c_str());
    if (!f.IsZombie())
      fFileVec.push_back(flux_file.c_str());
  }

  //set number of pot per file !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //
  fTotalPOT  = potperfile*(double)fFileVec.size();

  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //set the filename prefix for saving histogram plots 
  //
  std::string start_index_string;
  std::string end_index_string;
  
  std::ostringstream start_convert;
  start_convert << start_index;
  start_index_string = start_convert.str();
 
  std::ostringstream end_convert;
  end_convert << start_index + n_files;
  end_index_string = end_convert.str();
  
   // if someone is making files from someone else's flux files
  // he/she probably wants to write to his own directory
  // instead of someone else's
  std::string output_flux_dir = input_flux_dir;
  if(!on_grid) {
    std::string output_user = getenv("USER");
    output_flux_dir = "/dune/data/users/"+output_user+"/fluxfiles/g4lbne/"+version+"/"+physics_list+"/"+macro+"/"+current+"/flux/";
    system(("mkdir -p "+output_flux_dir).c_str());
  }

  ffilename = output_flux_dir + "histos_g4lbne_"+version+"_"+physics_list+"_"+macro+"_"+current;

  if(n_files==1)
    ffilename = output_flux_dir + "histos_g4lbne_"+version+"_"+physics_list+"_"+macro+"_"+current+"_"+start_index_string;

   //
   //????????????????????????????
   //????????????????????????????
  std::string ntuple_name = "nudata";
  if(isDk2nu)
    ntuple_name = "dk2nuTree";

  fChain = new TChain(ntuple_name.c_str());
   for(std::vector<std::string>::const_iterator sit = fFileVec.begin(); sit != fFileVec.end(); ++sit)
   {
     fChain -> Add(sit -> c_str());
   }

   Init(fChain);

   // initialize random numbers used for oscillation calculations
   rand3 = new TRandom3(0); 

}

eventRates::~eventRates()
{
   if (!fChain) return;
   delete fChain->GetCurrentFile();
}

Int_t eventRates::GetEntry(Long64_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Long64_t eventRates::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
   if (!fChain) return -5;
   Long64_t centry = fChain->LoadTree(entry);
   if (centry < 0) return centry;
   if (!fChain->InheritsFrom(TChain::Class()))  return centry;
   TChain *chain = (TChain*)fChain;
   if (chain->GetTreeNumber() != fCurrent) {
      fCurrent = chain->GetTreeNumber();
      Notify();
   }
   return centry;
}

void eventRates::Init(TTree *tree)
{
   // The Init() function is called when the selector needs to initialize
   // a new tree or chain. Typically here the branch addresses and branch
   // pointers of the tree will be set.
   // It is normally not necessary to make changes to the generated
   // code, but the routine can be extended by the user if needed.
   // Init() will be called many times when running on PROOF
   // (once per file to be processed).

   // Set branch addresses and branch pointers

   fCurrent = -1;
   
   

   if(!isDk2nu) {
     fChain->SetMakeClass(1);
     fChain->SetBranchAddress("run", &run, &b_data_run);
     fChain->SetBranchAddress("evtno", &evtno, &b_data_evtno);
     fChain->SetBranchAddress("protonN", &protonN, &b_data_protonN);
     fChain->SetBranchAddress("beamHWidth", &beamHWidth, &b_data_beamHWidth);
     fChain->SetBranchAddress("beamVWidth", &beamVWidth, &b_data_beamVWidth);
     fChain->SetBranchAddress("beamX", &beamX, &b_data_beamX);
     fChain->SetBranchAddress("beamY", &beamY, &b_data_beamY);
     fChain->SetBranchAddress("protonX", &protonX, &b_data_protonX);
   fChain->SetBranchAddress("protonY", &protonY, &b_data_protonY);
   fChain->SetBranchAddress("protonZ", &protonZ, &b_data_protonZ);
   fChain->SetBranchAddress("protonPx", &protonPx, &b_data_protonPx);
   fChain->SetBranchAddress("protonPy", &protonPy, &b_data_protonPy);
   fChain->SetBranchAddress("protonPz", &protonPz, &b_data_protonPz);
   fChain->SetBranchAddress("nuTarZ", &nuTarZ, &b_data_nuTarZ);
   fChain->SetBranchAddress("hornCurrent", &hornCurrent, &b_data_hornCurrent);
   fChain->SetBranchAddress("Ndxdz", &Ndxdz, &b_data_Ndxdz);
   fChain->SetBranchAddress("Ndydz", &Ndydz, &b_data_Ndydz);
   fChain->SetBranchAddress("Npz", &Npz, &b_data_Npz);
   fChain->SetBranchAddress("Nenergy", &Nenergy, &b_data_Nenergy);
   fChain->SetBranchAddress("NdxdzNear[5]", NdxdzNear, &b_data_NdxdzNear);
   fChain->SetBranchAddress("NdydzNear[5]", NdydzNear, &b_data_NdydzNear);
   fChain->SetBranchAddress("NenergyN[5]", NenergyN, &b_data_NenergyN);
   fChain->SetBranchAddress("NWtNear[5]", NWtNear, &b_data_NWtNear);
   fChain->SetBranchAddress("NdxdzFar[3]", NdxdzFar, &b_data_NdxdzFar);
   fChain->SetBranchAddress("NdydzFar[3]", NdydzFar, &b_data_NdydzFar);
   fChain->SetBranchAddress("NenergyF[3]", NenergyF, &b_data_NenergyF);
   fChain->SetBranchAddress("NWtFar[3]", NWtFar, &b_data_NWtFar);
   fChain->SetBranchAddress("Norig", &Norig, &b_data_Norig);
   fChain->SetBranchAddress("Ndecay", &Ndecay, &b_data_Ndecay);
   fChain->SetBranchAddress("Ntype", &Ntype, &b_data_Ntype);
   fChain->SetBranchAddress("Vx", &Vx, &b_data_Vx);
   fChain->SetBranchAddress("Vy", &Vy, &b_data_Vy);
   fChain->SetBranchAddress("Vz", &Vz, &b_data_Vz);
   fChain->SetBranchAddress("pdPx", &pdPx, &b_data_pdPx);
   fChain->SetBranchAddress("pdPy", &pdPy, &b_data_pdPy);
   fChain->SetBranchAddress("pdPz", &pdPz, &b_data_pdPz);
   fChain->SetBranchAddress("ppdxdz", &ppdxdz, &b_data_ppdxdz);
   fChain->SetBranchAddress("ppdydz", &ppdydz, &b_data_ppdydz);
   fChain->SetBranchAddress("pppz", &pppz, &b_data_pppz);
   fChain->SetBranchAddress("ppenergy", &ppenergy, &b_data_ppenergy);
   fChain->SetBranchAddress("ppmedium", &ppmedium, &b_data_ppmedium);
   fChain->SetBranchAddress("ptype", &ptype, &b_data_ptype);
   fChain->SetBranchAddress("ptrkid", &ptrkid, &b_data_ptrkid);
   fChain->SetBranchAddress("ppvx", &ppvx, &b_data_ppvx);
   fChain->SetBranchAddress("ppvy", &ppvy, &b_data_ppvy);
   fChain->SetBranchAddress("ppvz", &ppvz, &b_data_ppvz);
   fChain->SetBranchAddress("muparpx", &muparpx, &b_data_muparpx);
   fChain->SetBranchAddress("muparpy", &muparpy, &b_data_muparpy);
   fChain->SetBranchAddress("muparpz", &muparpz, &b_data_muparpz);
   fChain->SetBranchAddress("mupare", &mupare, &b_data_mupare);
   fChain->SetBranchAddress("Necm", &Necm, &b_data_Necm);
   fChain->SetBranchAddress("Nimpwt", &Nimpwt, &b_data_Nimpwt);
   fChain->SetBranchAddress("xpoint", &xpoint, &b_data_xpoint);
   fChain->SetBranchAddress("ypoint", &ypoint, &b_data_ypoint);
   fChain->SetBranchAddress("zpoint", &zpoint, &b_data_zpoint);
   fChain->SetBranchAddress("tvx", &tvx, &b_data_tvx);
   fChain->SetBranchAddress("tvy", &tvy, &b_data_tvy);
   fChain->SetBranchAddress("tvz", &tvz, &b_data_tvz);
   fChain->SetBranchAddress("tpx", &tpx, &b_data_tpx);
   fChain->SetBranchAddress("tpy", &tpy, &b_data_tpy);
   fChain->SetBranchAddress("tpz", &tpz, &b_data_tpz);
   fChain->SetBranchAddress("tptype", &tptype, &b_data_tptype);
   fChain->SetBranchAddress("tgen", &tgen, &b_data_tgen);
   }

   else {
     dk2nu = new bsim::Dk2Nu;
     fChain->SetBranchAddress("dk2nu",&dk2nu);
     fChain->GetEntry(0);
   }
   Notify();
}

Bool_t eventRates::Notify()
{
   // The Notify() function is called when a new file is opened. This
   // can be either for a new TTree in a TChain or when when a new TTree
   // is started when using PROOF. It is normally not necessary to make changes
   // to the generated code, but the routine can be extended by the
   // user if needed. The return value is currently not used.

   return kTRUE;
}

void eventRates::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}
Int_t eventRates::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}
#endif // #ifdef eventRates_cxx