estimate_cosmics.C

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////                                                                                                      
// estimate_cosmics.C                                                     //
//                                                                        //
// follows get_cosmic_sample.C                                            //                                                                                                     
// uses cosmic trigger sample for spectrum shape (more stats)             //
// uses numi data for the count (real fd conditions from sidebands)       //
// if only using the cosmic trigger data, run *_triggeronly.C             // 
//                                                                        //
// Script by:      Diana Mendez (d.p.mendez@sussex.ac.uk) April 4 2018    //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include "3FlavorAna/Ana2018/numu/includes.h"

using namespace ana;

void estimate_cosmics(std::string horn  = "rhc")
{

  std::cout << "\nEstimating " << horn << " cosmics" << std::endl;
  std::cout << "==========================   " << std::endl;
  std::cout << "  *\n    *\n      *        FD" << std::endl;
  std::cout << "        *\n          *       " << std::endl;
  std::cout << "            *                " << std::endl;
  std::cout << "==========================\n " << std::endl;

  double pot, livetime;
  std::string period_numi = "full", period_cosmic = "full";

  if(horn=="fhc"){
    livetime = kAna2018FHCLivetime;
    pot      = kAna2018FHCPOT;
  }
  else if(horn=="rhc"){
    livetime = kAna2018RHCLivetime;
    pot      = kAna2018RHCPOT;
  }
  else{
    std::cout << "=== unrecognised horn " << horn << "===" << std::endl;
    exit (EXIT_FAILURE);
  }


  std::string file_name_numi   = "results/cosmicnumi_" + horn + "_" + period_cosmic + "__numu2018.root";
  std::string file_name_cosmic = "results/cosmiccosmic_" + horn + "_" + period_cosmic + "__numu2018.root";
  TFile *fnumi   = TFile::Open(file_name_numi.c_str());
  TFile *fcosmic = TFile::Open(file_name_cosmic.c_str());


  // we'll write original numi and trigger histograms in the output
  // so clone hcosmic for final estimate (instead of directly scaling the trigger)
  std::vector<TH1*> hnumi;
  std::vector<TH1*> hcosmic;
  std::vector<TH1*> hestimate;
  std::cout << "Getting histograms from files" << std::endl;
  for(int quantId = 0; quantId < 5; quantId++){
    hnumi.push_back((TH1*)fnumi->Get(Form("numuEOpt_q%d", quantId)));
    hcosmic.push_back((TH1*)fcosmic->Get(Form("numuEOpt_q%d", quantId)));
    TH1* htemp = (TH1*)hcosmic[quantId]->Clone(); hestimate.push_back(htemp);
  }
  

  // dont trust the count in the files (to do: check)
  // get the count directly from the integral
  std::cout << "Getting pot, livetime, count" << std::endl;  
  TH1* hlivetime_numi  = (TH1*)fnumi->Get("livetime"); float livetime_numi  = hlivetime_numi->Integral();
  TH1* hlivetime_cosmic  = (TH1*)fcosmic->Get("livetime"); float livetime_cosmic  = hlivetime_cosmic->Integral();

  // note livetime scaling:
  // livetime sideband is not equal to livetime from beamspill which is what analysis call livetime
  float scale = livetime / livetime_numi;
  std::cout << "\nsideband livetime = " << livetime_numi << ", ana livetime = " << livetime << " -> scale = " << scale << std::endl;


  std::cout << "Getting counts per quantile" << std::endl;
  std::vector<float> count_numi, count_cosmic;
  for(int quantId = 0; quantId < 5; quantId++){
    float counttemp;
    counttemp = hnumi[quantId]->Integral(); count_numi.push_back(counttemp);
    counttemp = hcosmic[quantId]->Integral(); count_cosmic.push_back(counttemp);
  }

  std::vector<float> events;
  std::vector<float> errors;
  for(int quantId=0; quantId<5; quantId++){
    events.push_back( count_numi[quantId] * scale ); 
    errors.push_back( events[quantId] * sqrt(count_numi[quantId]) / count_numi[quantId] );
    if(quantId==0) std::cout << "\nTotal estimate = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_numi[quantId] << " raw events" << std::endl;
    else std::cout << "quantile " << quantId <<  " = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_numi[quantId] << " raw events" << std::endl;
  }

  // we get the shape from the trigger and the count from numi
  // now scale the trigger to number of events
  std::cout << "\nChecking cosmic integrals are right" << std::endl;
  std::cout << "(should agree with event estimates)" << std::endl;
  for(int quantId=0; quantId<5; quantId++){
    hestimate[quantId]->Scale( events[quantId] / (hestimate[quantId]->Integral()) );
    if(quantId==0) std::cout << "Total integral = " << hestimate[quantId]->Integral() << std::endl;
    else std::cout << "quantile " << quantId << " = " << hestimate[quantId]->Integral() << std::endl;
  }
  
  std::cout << "\nWriting to output file" << std::endl;
  std::string name_out = "cosmics_" + horn + "__numu2018.root";
  TFile fout(name_out.c_str(),"recreate");
  for(int quantId=0; quantId<5; quantId++){
    hnumi[quantId]->Write(Form("numi_q%d", quantId));
    hcosmic[quantId]->Write(Form("trigger_q%d", quantId));
    hestimate[quantId]->Write(Form("cosmics_q%d", quantId));
  }
  
  fout.Close();

  std::cout << "\nDone!" << std::endl;

}