studyAncestorsL.cc

Go to the documentation of this file.

// 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>

#include <TStyle.h>

#include <string>
#include <iostream>
#include "TROOT.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TColor.h"

#include "dk2nu/tree/dk2nu.h"
#include "dk2nu/tree/dkmeta.h"  
 
class ancestor {
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;

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

   ancestor();

   virtual ~ancestor();
   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);

 private:

   
};

ancestor::ancestor()
{

  /*
  // 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;
  }
  */

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

  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00250.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00249.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00248.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00247.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00246.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00245.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00244.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00243.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00242.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00241.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_NoSnout_neutrino_00240.dk2nu.root");
  fFileVec.push_back("/pnfs/dune/scratch/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017_NoSnout/neutrino/flux/g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_neutrino_00239.dk2nu.root");


  std::string 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); 

}

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

Int_t ancestor::GetEntry(Long64_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Long64_t ancestor::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 ancestor::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;
   
   dk2nu = new bsim::Dk2Nu;
   fChain->SetBranchAddress("dk2nu",&dk2nu);
   fChain->GetEntry(0);

   Notify();
}

Bool_t ancestor::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 ancestor::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}
Int_t ancestor::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}

void ancestor::Loop()
{


  dk2nu = new bsim::Dk2Nu;
   fChain->SetBranchAddress("dk2nu",&dk2nu);

   double z = 2267; // z position where you want to plot hadron parent distributions

   TH2D *h2_all = new TH2D("h_all","",20,-100,100,20,-100,100);
   TH2D *h2_fardet = new TH2D("h_fardet","",20,-100,100,20,-100,100);
   TH2D *h2_fardet_peak = new TH2D("h_fardet_peak","",20,-100,100,20,-100,100);
   TH1D *hx_fardet = new TH1D("hx_fardet","",40,-100,100);
   TH1D *hy_fardet = new TH1D("hy_fardet","",40,-100,100);

   TH2D *h_vxy_all = new TH2D("h_vxy_all","",20,-100,100,20,-100,100);
   TH2D *h_vxy = new TH2D("h_vxy","",20,-100,100,20,-100,100);
   TH2D *h_vxy_peak = new TH2D("h_vxy_peak","",20,-100,100,20,-100,100);
   TH1D *h_vx = new TH1D("h_vx","",40,-100,100);     
   TH1D *h_vy = new TH1D("h_vy","",40,-100,100);     
   TH1D *h_vx_peak = new TH1D("h_vx_peak","",40,-100,100);     
   TH1D *h_vy_peak = new TH1D("h_vy_peak","",40,-100,100);     

   TH1D *h_fluxfardet = new TH1D("h_fluxfardet","",40,0,20);

   TCanvas *c = new TCanvas();

   int iread = 0;
   Long64_t nentries = fChain->GetEntries();
   std::cout << "Total number of Entries = " << nentries << std::endl;   

   double ancestors_tot = 0;
   double ancestors_10 = 0;
   double ancestors_20 = 0;
   double ancestors_30 = 0;
   double ancestors_40 = 0;
   double ancestors_50 = 0;
   double ancestors_60 = 0;
   double ancestors_70 = 0;
   double ancestors_77p5 = 0;
   double ancestors_80 = 0;
   double ancestors_90 = 0;
   double ancestors_100 = 0;

   for (Long64_t jentry=0; jentry<nentries;jentry++) 
   {
      fChain->GetEntry(jentry);  
      ++iread;
      
      if(iread % 10000 == 0)
      {
         std::cout << "Reading Entry " << iread << std::endl;
      }


      double Nimpwt = dk2nu->decay.nimpwt;
      int Ntype = dk2nu->decay.ntype;
      int Nancestors = dk2nu->ancestor.size();
      double fardetwt = dk2nu->nuray[2].wgt;

      h_fluxfardet->Fill(dk2nu->nuray[2].E,Nimpwt*fardetwt);

      if(dk2nu->decay.vz>z-100 && dk2nu->decay.vz<z+100) {
        h_vxy_all->Fill(dk2nu->decay.vx,dk2nu->decay.vy,Nimpwt);
        h_vxy->Fill(dk2nu->decay.vx,dk2nu->decay.vy,Nimpwt*fardetwt);
        h_vx->Fill(dk2nu->decay.vx,Nimpwt*fardetwt);
        h_vy->Fill(dk2nu->decay.vy,Nimpwt*fardetwt);

        if(dk2nu->nuray[2].E > 0.5 &&
           dk2nu->nuray[2].E < 4.0 && 
           dk2nu->decay.ntype==14) {
          h_vxy_peak->Fill(dk2nu->decay.vx,dk2nu->decay.vy,Nimpwt*fardetwt);
          h_vx_peak->Fill(dk2nu->decay.vx,Nimpwt*fardetwt);
          h_vy_peak->Fill(dk2nu->decay.vy,Nimpwt*fardetwt);
        }
      }

      for(int i = 0; i<Nancestors; i++) {
        // Skip neutrinos -- they don't care if the snout is there
        if(dk2nu->ancestor[i].pdg == 12 ||
           dk2nu->ancestor[i].pdg == -12 ||
           dk2nu->ancestor[i].pdg == 14 ||
           dk2nu->ancestor[i].pdg == -14)
          continue;


        if(dk2nu->ancestor[i].startz <= z &&
           (i+1 == Nancestors || dk2nu->ancestor[i+1].startz>z)) {
          double projx = dk2nu->ancestor[i+1].startx - ((dk2nu->ancestor[i+1].startz-z)*dk2nu->ancestor[i].stoppx/dk2nu->ancestor[i].stoppz);
          double projy = dk2nu->ancestor[i+1].starty - ((dk2nu->ancestor[i+1].startz-z)*dk2nu->ancestor[i].stoppy/dk2nu->ancestor[i].stoppz);
          h2_all->Fill(projx,projy,Nimpwt);
          h2_fardet->Fill(projx,projy,Nimpwt*fardetwt*1);
          hx_fardet->Fill(projx,Nimpwt*fardetwt*1);
          hy_fardet->Fill(projy,Nimpwt*fardetwt*1);

          if(dk2nu->nuray[2].E > 0.5 &&
             dk2nu->nuray[2].E < 4.0 && 
             dk2nu->decay.ntype==14) {
            h2_fardet_peak->Fill(projx,projy,Nimpwt*fardetwt*1);
            ancestors_tot += Nimpwt*fardetwt;
            double radius = sqrt(projx*projx+projy*projy);
            if(radius<10)
              ancestors_10 += Nimpwt*fardetwt;
            if(radius<20)
              ancestors_20 += Nimpwt*fardetwt;
            if(radius<30)
              ancestors_30 += Nimpwt*fardetwt;
            if(radius<40)
              ancestors_40 += Nimpwt*fardetwt;
            if(radius<50)
              ancestors_50 += Nimpwt*fardetwt;
            if(radius<60)
              ancestors_60 += Nimpwt*fardetwt;
            if(radius<70)
              ancestors_70 += Nimpwt*fardetwt;
            if(radius<77.5)
              ancestors_77p5 += Nimpwt*fardetwt;
            if(radius<80)
              ancestors_80 += Nimpwt*fardetwt;
            if(radius<90)
              ancestors_90 += Nimpwt*fardetwt;
            if(radius<100)
              ancestors_100 += Nimpwt*fardetwt;
          }

        }
        
      }
     
   }//end loop over entries

   std::cout<<"Total ancestors "<<ancestors_tot<<std::endl;
   std::cout<<"Total within 10 cm "<<ancestors_10/ancestors_tot<<std::endl;
   std::cout<<"Total within 20 cm "<<ancestors_20/ancestors_tot<<std::endl;
   std::cout<<"Total within 30 cm "<<ancestors_30/ancestors_tot<<std::endl;
   std::cout<<"Total within 40 cm "<<ancestors_40/ancestors_tot<<std::endl;
   std::cout<<"Total within 50 cm "<<ancestors_50/ancestors_tot<<std::endl;
   std::cout<<"Total within 60 cm "<<ancestors_60/ancestors_tot<<std::endl;
   std::cout<<"Total within 70 cm "<<ancestors_70/ancestors_tot<<std::endl;
   std::cout<<"Total within 77.5 cm "<<ancestors_77p5/ancestors_tot<<std::endl;
   std::cout<<"Total within 80 cm "<<ancestors_80/ancestors_tot<<std::endl;
   std::cout<<"Total within 90 cm "<<ancestors_90/ancestors_tot<<std::endl;
   std::cout<<"Total within 100 cm "<<ancestors_100/ancestors_tot<<std::endl;

   
   std::string str_z = std::to_string(z);
   h2_all->SetTitle("All Neutrinos");
   h2_fardet->SetTitle("Neutrinos at Far Detector");
   h2_fardet_peak->SetTitle("Peak #nu_#mu neutrinos at Far Detector");
   h_vxy_all->SetTitle("All Neutrino");
   h_vxy->SetTitle("Neutrinos at Far Detector");
   h_vx->SetTitle("Neutrinos at Far Detector");
   h_vy->SetTitle("Neutrinos at Far Detector");
   h_vxy_peak->SetTitle("Peak #nu_#mu neutrinos at Far Detector");
   h_vx_peak->SetTitle("Peak #nu_#mu neutrinos at Far Detector");
   h_vy_peak->SetTitle("Peak #nu_#mu neutrinos at Far Detector");
   h2_all->GetXaxis()->SetTitle(("Ancestor X position (cm) at z = "+str_z).c_str());
   h2_all->GetYaxis()->SetTitle(("Ancestor Y position (cm) at z = "+str_z).c_str());
   h2_fardet->GetXaxis()->SetTitle(("Ancestor X position (cm) at z = "+str_z).c_str());
   h2_fardet->GetYaxis()->SetTitle(("Ancestor Y position (cm) at z = "+str_z).c_str());
   h2_fardet_peak->GetXaxis()->SetTitle(("Ancestor X position (cm) at z = "+str_z).c_str());
   h2_fardet_peak->GetYaxis()->SetTitle(("Ancestor Y position (cm) at z = "+str_z).c_str());
   h_vxy->GetXaxis()->SetTitle(("Neutrino parent decay X position (cm) at z = "+str_z).c_str());
   h_vxy_all->GetXaxis()->SetTitle(("Neutrino parent decay X position (cm) at z = "+str_z).c_str());
   h_vx->GetXaxis()->SetTitle(("Neutrino parent decay X position (cm) at z = "+str_z).c_str());
   h_vy->GetXaxis()->SetTitle(("Neutrino parent decay Y position (cm) at z = "+str_z).c_str());
   h_vxy->GetYaxis()->SetTitle(("Neutrino parent decay Y postion (cm) at z = "+str_z).c_str());
   h_vxy_all->GetYaxis()->SetTitle(("Neutrino parent decay Y postion (cm) at z = "+str_z).c_str());
   h_vxy_peak->GetXaxis()->SetTitle(("Neutrino parent decay position (cm) at z = "+str_z).c_str());
   h_vx_peak->GetXaxis()->SetTitle(("Neutrino parent decay X position (cm) at z = "+str_z).c_str());
   h_vy_peak->GetXaxis()->SetTitle(("Neutrino parent decay Y position at (cm) z = "+str_z).c_str());
   h_vxy_peak->GetYaxis()->SetTitle(("Neutrino parent decay postion at (cm) z = "+str_z).c_str());

   gStyle->SetOptStat(0);
   gStyle->SetPalette(56);

   h2_all->SetContour(99);
   h2_all->Draw("COLZ");
   c->Print("Ancestors_Through_Window_AllNu.eps");
   c->Print("Ancestors_Through_Window_AllNu.png");
   h2_fardet->Draw("COLZ");
   c->Print("Ancestors_Through_Window_FarDetNu_2D.eps");
   //   c->Print("Ancestors_Through_Window_FarDetNu_2D.png");
   h2_fardet_peak->Draw("COLZ");
   c->Print("Ancestors_Through_Window_FarDetNuPeak_2D.eps");
   //c->Print("Ancestors_Through_Window_FarDetNuPeak_2D.png");
   hx_fardet->Draw();
   c->Print("Ancestors_Through_Window_FarDetNu_x.eps");
   c->Print("Ancestors_Through_Window_FarDetNu_x.png");
   hy_fardet->Draw();
   c->Print("Ancestors_Through_Window_FarDetNu_y.eps");
   c->Print("Ancestors_Through_Window_FarDetNu_y.png");
   h_fluxfardet->Draw();
   c->Print("Flux_fardet.eps");
   c->Print("Flux_fardet.png");


   h_vxy->Draw("COLZ");
   c->Print("Decay_Position_Near_Window_FarDetNu_2D.eps");
   //c->Print("Decay_Position_Near_Window_FarDetNu_2D.png");  

   h_vxy_all->Draw("COLZ");
   c->Print("Decay_Position_Near_Window_AllNu_2D.eps");
   c->Print("Decay_Position_Near_Window_AllNu_2D.png");  

   h_vxy_peak->Draw("COLZ");
   c->Print("Decay_Position_Near_Window_FarDetNuPeak_2D.eps");
   //c->Print("Decay_Position_Near_Window_FarDetNuPeak_2D.png");  

   h_vx->Draw("");
   c->Print("Decay_Position_Near_Window_FarDetNu_x.eps");
   c->Print("Decay_Position_Near_Window_FarDetNu_x.eps");
   h_vy->Draw("");
   c->Print("Decay_Position_Near_Window_FarDetNu_y.eps");
   c->Print("Decay_Position_Near_Window_FarDetNu_y.eps");

   h_vx_peak->Draw("");
   c->Print("Decay_Position_Near_Window_FarDetNuPeak_x.eps");
   c->Print("Decay_Position_Near_Window_FarDetNuPeak_x.eps");
   h_vy_peak->Draw("");
   c->Print("Decay_Position_Near_Window_FarDetNuPeak_y.eps");
   c->Print("Decay_Position_Near_Window_FarDetNuPeak_y.eps");

}




int main(int argc, char* argv[]) {

  ancestor t;
  std::cout<< "BLAH"<<std::endl;
  t.Loop();  

}