NueAnaReader.cxx

Go to the documentation of this file.

#include <iostream>
#include "TChain.h"
#include "TFile.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "NueAnaReader.h"
#include "NueAna/NuePOT.h"
#include "NueAna/MCNNBestMatch.h"
#include <fstream>

using std::cout;
using std::endl;

const Float_t NueAnaReader::POT_file_carrot_FD=2.91*1e20;
const Float_t NueAnaReader::POT_file_daikon_FD=6.5*1e20;
const Float_t NueAnaReader::POT_file_carrot_ND=1.0024*1e16;
const Float_t NueAnaReader::POT_file_daikon_ND=1.936*1e16;

// ****************** CONSTRUCTOR & DESTRUCTOR ****************

NueAnaReader::NueAnaReader(){

  //Make chains
  chain = new TChain("ana_nue");
  chain_pot = new TChain("pottree");

  nuerecord = new NueRecord();
  nuepot = new NuePOT();

  //defaults
  num_best_matches=50;
  h_weight=new TH1F();
  //weightfile = new TFile("./files/MakeHistosForSpectralReweighting.root","READ");

}

NueAnaReader::NueAnaReader(const char* name){

  //Make chains
  chain = new TChain("ana_nue");
  chain_pot = new TChain("pottree");

  filenum=chain->Add(name);
  chain_pot->Add(name);

  cout  << "Added " << filenum << " to NueAnaReader" << endl;

  nuerecord = new NueRecord();
  nuepot = new NuePOT();

  Setup();

  //defaults
  num_best_matches=50;  
  h_weight=new TH1F();

}

NueAnaReader::~NueAnaReader(){

  delete chain;
  delete chain_pot;
  delete nuerecord;
  delete nuepot;

}

// ******** SETUP (PRIVATE METHOD)
void NueAnaReader::Setup(){

  //Branches
  //--------
  if(chain->GetBranchStatus("NueRecord")){
    chain->SetBranchAddress("NueRecord",&nuerecord);
  }
  if(chain_pot->GetBranchStatus("NuePOT")){
    chain_pot->SetBranchAddress("NuePOT",&nuepot);
    chain_pot->GetEntry(0);
  }
     
}
  
// ****** Public Methods *************
Int_t NueAnaReader::Add(const char* name){

  Int_t filenum = chain->Add(name);
  cout << "Added " << chain_pot->Add(name) << " to the pottree chain" << endl;
  Setup();
  return filenum;

}

//----------------------------------------
Int_t NueAnaReader::GetEntry(Int_t srCtr_in){

  return chain->GetEntry(srCtr_in);

}

//----------------------------------------
Int_t NueAnaReader::GetEntries(){

  return chain->GetEntries();

}
//----------------------------------------
Bool_t NueAnaReader::NuePresel(){

  Bool_t passpresel=false;
  if(nuerecord!=NULL){
    if(nuerecord->anainfo.inFiducialVolume == 1){// && nuerecord->srevent.vtxZ<29){
       if(nuerecord->srevent.phNueGeV > 1 && nuerecord->srevent.phNueGeV < 8){
         if(nuerecord->srtrack.planes < 25){
           if(nuerecord->srtrack.trklikePlanes < 16){
             passpresel=true;
           }
         }
       }
     }
  } else {
    cout  << "Warning ! Emtpy NueRecord in NuePresel" << endl;
  }

  return passpresel;

}

//--------------------------------------
Float_t NueAnaReader::CalculatePOT(){

  Int_t srCtr_pot=0;
  Float_t pot_out=0;
  NuePOT *nuepot = new NuePOT();
  chain_pot->SetBranchAddress("NuePOT",&nuepot);

  while(chain_pot->GetEntry(srCtr_pot)){
    ++srCtr_pot;
    pot_out += nuepot->pot;
  }

  delete nuepot;

  return pot_out;

}

//--------------------------------------
Float_t NueAnaReader::CalculatePOT(Int_t filenum, Int_t filetype){

  Float_t pot_out=0;

  if(filetype==0){//<--carrot FD
    pot_out=filenum*POT_file_carrot_FD*1e-12;
    cout << "Using carrot FD POT" << endl;
  }
  if(filetype==1){//<--daikon FD
    pot_out=filenum*POT_file_daikon_FD*1e-12;
    cout << "Using daikon FD POT " << endl;
  }
  if(filetype==2){//<--carrot ND
    pot_out=filenum*POT_file_carrot_ND*1e-12;
    cout << "Using carrot ND POT" << endl;
  }
  if(filetype==3){//<--daikon ND
    pot_out=filenum*POT_file_daikon_ND*1e-12;
    cout << "Using daikon ND POT " << endl;
  }
  
  return pot_out;

}

//--------------------------------------------
MCNNBestMatch* NueAnaReader::GetBestMatch(Int_t ii){

  MCNNBestMatch *castedbmatch = 0;
  TClonesArray& bmatcharray = *(nuerecord->mcnnv.bmatch);
  castedbmatch = dynamic_cast<MCNNBestMatch *>(bmatcharray[ii]);
  return castedbmatch;


}

//---------------------------------------------
void NueAnaReader::SetMaxNumBestMatches(Int_t nbm_in){
  
  num_best_matches=nbm_in;

}

//----------------------------------------------
Float_t NueAnaReader::GetfracCC(){

  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Int_t match_nue=0;
  Int_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1){
      ++match_nue;
    }
    ++matches;
    ++bm;
  }
  
  Float_t fracCC=-1;

  if(matches>0){
    fracCC=match_nue*1./matches;
  }
  
  return fracCC;

}


//----------------------------------------------
Float_t NueAnaReader::GetfracCCy(Float_t y_cut){

  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Int_t match_nue=0;
  Int_t matches=0;
  Int_t bm=0;

  while(matches<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);

    //Int_t filetype=mcnnbm->run/100000;
    //Int_t max_run=0;
    //if(filetype==210){
    //  max_run=21011200;
    // } else {
    //  max_run=21111100;
    // }
    
    //if(mcnnbm->run <= max_run){
      
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y<y_cut){
        ++match_nue;
      }
      ++matches;
      //}

    ++bm;
  }
  
  Float_t fracCCy=-1;

  if(matches>0){
    fracCCy=match_nue*1./matches;
  }
  
  return fracCCy;

}


//---------------------------------------------------
Float_t NueAnaReader::Getymean(Float_t y_cut){
  
  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t ely=0;
  Int_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y < y_cut){
      ely += mcnnbm->true_y;
      ++matches;
    }
    ++bm;
  }
  
  Float_t result=-1;

  if(matches>0){
    result=ely*1./matches;
  }
  
  return result;

}

//------------------------------------------------------
Float_t NueAnaReader::GetMeanFracQMatched(Float_t y_cut){
  
  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t mfqm=0;
  Int_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y < y_cut){
      mfqm += mcnnbm->fracQmatched;
      ++matches;
    }
    ++bm;
  }
  
  Float_t result=-1;

  if(matches>0){
    result=mfqm*1./matches;
  }
  
  return result;

}

//------------------------------------------------------
Float_t NueAnaReader::GetdlnLmean(){
  
  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t mfqm=0;
  Int_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1){
      mfqm += mcnnbm->dlnL;
      ++matches;
    }
    ++bm;
  }
  
  Float_t result=-1;

  if(matches>0){
    result=mfqm*1./matches;
  }
  
  return result;

}

//----------------------------------------------
void NueAnaReader::SetLibraryWeightingFile(char* weightfile_in, char* histname_in){

  TFile *weightfile = new TFile(weightfile_in,"READ");
  h_weight = (TH1F*)weightfile->Get(histname_in);

}

//----------------------------------------------
Float_t NueAnaReader::GetWeightedfracCC(){

  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t match_nue=0;
  Float_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);

    Float_t weight=1.0;
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1) weight=h_weight->GetBinContent(h_weight->FindBin(mcnnbm->true_enu));
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1){
      match_nue+=weight;;
    }
    matches+=weight;
    ++bm;
  }
  
  Float_t fracCC=-1;

  if(matches>0){
    fracCC=match_nue*1./matches;
  }
  
  return fracCC;

}

//----------------------------------------------
Float_t NueAnaReader::GetWeightedfracCCy(Float_t y_cut){

  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t match_nue=0;
  Float_t matches=0;
  Int_t bm=0;

  while(matches<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);

    //Int_t filetype=mcnnbm->run/100000;
    //Int_t max_run=0;
    //if(filetype==210){
    //  max_run=21011200;
    // } else {
    //  max_run=21111100;
    // }
    
    //if(mcnnbm->run <= max_run){

    //Setting a non-1 weight for nueCC best matches
    Float_t weight=1.0;
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1) weight=h_weight->GetBinContent(h_weight->FindBin(mcnnbm->true_enu));
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y<y_cut){
      match_nue+=weight;
    }
    matches+=weight;
      

      //}

    ++bm;
  }
  
  Float_t fracCCy=-1;

  if(matches>0){
    fracCCy=match_nue*1./matches;
  }
  
  return fracCCy;

}

//------------------------------------------------------
Float_t NueAnaReader::GetWeightedMeanFracQMatched(Float_t y_cut){
  
  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t mfqm=0;
  Float_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);

    Float_t weight=1.0;
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1) weight=h_weight->GetBinContent(h_weight->FindBin(mcnnbm->true_enu));
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y < y_cut){
      mfqm += mcnnbm->fracQmatched*weight;
      matches+=weight;
    }
    ++bm;
  }
  
  Float_t result=-1;

  if(matches>0){
    result=mfqm*1./matches;
  }
  
  return result;

}

//---------------------------------------------------
Float_t NueAnaReader::GetWeightedymean(Float_t y_cut){
  
  Int_t nresults=nuerecord->mcnnv.bestmatches;
  Float_t ely=0;
  Float_t matches=0;
  Int_t bm=0;

  while(bm<num_best_matches && bm<nresults){
    
    MCNNBestMatch *mcnnbm = GetBestMatch(bm);

    Float_t weight=1.0;
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1) weight=h_weight->GetBinContent(h_weight->FindBin(mcnnbm->true_enu));
    
    if(TMath::Abs(mcnnbm->nuFlavor)==12 && mcnnbm->interactionType==1 && mcnnbm->true_y < y_cut){
      ely += mcnnbm->true_y*weight;
      matches+=weight;
    }
    ++bm;
  }
  
  Float_t result=-1;

  if(matches>0){
    result=ely*1./matches;
  }
  
  return result;

}

---