readNtuple.C

Go to the documentation of this file.

#define readNtuple_cxx
#include "readNtuple.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "TLorentzVector.h"
#include <iostream>
#include "TH3F.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TTree.h"
#include "string.h"
#include "TString.h"

using namespace std;
using namespace TMVA;

void readNtuple::Loop(const char* TypeNameNu, const char* TypeNameInt)
{

  bool ccnumu=false;
  bool ccnumubar=false;
  bool ccnue=false;
  bool ccnuebar=false;
  bool nc=false;

  TString nutype = (TString)TypeNameNu;
  if( nutype=="ccnumu" ){
    ccnumu=true;
    cout<<"Analyze the CC numu, ";
  } 
  else if( nutype=="ccnumubar" ){
    ccnumubar=true;
    cout<<"Analyze the CC numubar, ";
  }
  else if( nutype=="ccnue" ){
    ccnue=true;
    cout<<"Analyze the CC nue, ";
  }
  else if( nutype=="ccnuebar") {
    ccnuebar=true;
    cout<<"Analyze the CC nuebar, ";
  }
  else if( nutype=="nc" ){
     nc=true;
     cout<<"Analyze the NC, ";
  }
  else{
    cout<<"Please enter the right type of nu interactions: ccnumu, ccnumubar, ccnue, ccnuebar, nc"<<endl;
    exit(0);
  }

  bool inclusive=false;
  bool qeonly=false;
  bool disonly=false;
  bool resonly=false;
  bool meconly=false;
  bool cohonly=false;
  
  TString protype=(TString)TypeNameInt;
  if( protype=="all" ){
     inclusive=true;
     cout<<"inclusive events."<<endl;
  }
  else if( protype=="qe"){
     qeonly=true;
     cout<<"QE only events."<<endl;
  }
  else if( protype="res"){
     resonly=true;
     cout<<"RES only events."<<endl;
  }
  else if( protype=="dis"){
     disonly=true;
     cout<<"DIS only events."<<endl;
  }
  else if( protype="mec"){
     meconly=true;
     cout<<"MEC only events."<<endl;
  }
  else if( protype="coh"){
     cohonly=true;
     cout<<"COH only events."<<endl;
  }
  else{
    cout<<"please enerty the right interaction process: all, qe, res, dis, mec, coh"<<endl;
    exit(0);
  }
  
  if (fChain == 0) return;

  Int_t nentries = fChain->GetEntries();
  cout<<"Total events are "<<nentries<<endl;
  
  //read-in the flux spectrum
  TFile *flux=new TFile("doc/nu_flux_v08.root");
  TString fluxhist="";
  if( ccnumu ) fluxhist="numu_flux";
  else if( ccnumubar ) fluxhist="numubar_flux";
  else if( ccnue ) fluxhist="nue_flux";
  else if( ccnuebar ) fluxhist="nuebar_flux";
  else fluxhist="numu_flux";

  TH1F *hflux=(TH1F*)flux->Get(fluxhist);//nu flux in unit of per E10 POT per m2
  double integrate_flux=hflux->Integral();

  //# of nucleons in fully active detector
  //the Z is cut off at 12.75m
  double Nnucleons=11.5279;//unit E31

  //vertex cuts
  double minVz=0.;//cm
  double maxVz=1275.0;//cm
  double minVy=-150.0;//cm
  double maxVy=150.0; //cm
  double minVx=-150.0;//cm
  double maxVx=150.0;//cm

  double corr_Nnucleons= (maxVz-minVz)*(maxVy-minVy)*(maxVx-minVx)/(1275.0*400.0*400.0);
  Nnucleons = Nnucleons*corr_Nnucleons;

  cout<<"# of nucleons are "<<Nnucleons<<endl;

  //pot of MC
  double mcpot=0.00001;//E21

  int Nevts=0;
  int Nfinal=0;
 
  for( int ii=0; ii<nentries; ii++)
  //for( int ii=0; ii<10000; ii++)
    {
      fChain->GetEntry(ii);
      ++Nevts;
      if( Nevts%100000==0 ) cout<<Nevts<<" processed ..."<<endl;

      if( ccnumu ){
        if( nu_ccnc != 0 ) continue;
        if( nu_PDG != 14 ) continue;
        if( nu_origPDG != 14 ) continue;
      } 
      else if( ccnumubar ){
        if( nu_ccnc != 0 ) continue;
        if( nu_PDG != -14 ) continue;
        if( nu_origPDG != -14 ) continue;
      }
      else if( ccnue ){
        if( nu_ccnc != 0 ) continue;
        if( nu_PDG != 12 ) continue;
        if( nu_origPDG != 12 ) continue;
      }
      else if( ccnuebar ){
        if( nu_ccnc != 0 ) continue;
        if( nu_PDG != -12 ) continue;
        if( nu_origPDG != -12 ) continue;
      }
      else if( nc ){
        if( nu_ccnc != 1 ) continue;
        if( nu_PDG != nu_origPDG ) continue;
      }

      if( qeonly ){
        if( nu_mode != 0 ) continue;
      }
      else if( resonly ){
        if( nu_mode != 1 ) continue;
      }
      else if( disonly ){
        if( nu_mode != 2 ) continue;
      }
      else if( meconly ){
        if( nu_mode != 10 ) continue;
      }
      else if( cohonly ){
        if( nu_mode != 3 ) continue;
      }

      if( nu_Vx>maxVx || nu_Vx<minVx ) continue; 
      if( nu_Vy>maxVy || nu_Vy<minVy ) continue; 
      if( nu_Vz>maxVz || nu_Vz<minVz ) continue; 

      int Xbin=hflux->FindBin(nu_trueE);
      double flux_Ebin=0.0;
      if( hflux->GetBinContent(Xbin)>0. ) flux_Ebin=1.0/hflux->GetBinContent(Xbin);
      double total_weight=flux_Ebin/(mcpot*Nnucleons);
      double diff_weight=1.0/(integrate_flux*mcpot*Nnucleons); 

      myhists->Xsec_Enu->Fill(nu_trueE,total_weight);
      myhists->diffXsec_lepE->Fill(lep_E,diff_weight);
      TLorentzVector lepton;
      lepton.SetPxPyPzE(lep_Px,lep_Py,lep_Pz,lep_E);
      myhists->diffXsec_lepcostheta->Fill(cos(lepton.Eta()),diff_weight); 

      ++Nfinal;
    }//loop all events

    //dividing the bin width for diff. Xsec
    for( int ibin=1; ibin<myhists->diffXsec_lepE->GetNbinsX()+1; ++ibin ){
      double binwidth=myhists->diffXsec_lepE->GetBinWidth(ibin);
      double bincon=myhists->diffXsec_lepE->GetBinContent(ibin)/binwidth;
      double binerr=0.;
      if( bincon>0.0 ) binerr=(myhists->diffXsec_lepE->GetBinError(ibin)/myhists->diffXsec_lepE->GetBinContent(ibin))*bincon;
      myhists->diffXsec_lepE->SetBinContent(ibin,bincon);
      myhists->diffXsec_lepE->SetBinError(ibin,binerr);
    }
  
  cout<<"# selected events are "<<Nfinal<<endl;
  
}