estimate_cosmics_triggeronly.C

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////                                                                                                      
// estimate_cosmics_triggeronly.C                                         //                                                                                                      
//                                                                        //
// follows get_cosmic_sample.C                                            //                                                                                                     
// uses cosmic trigger sample for spectrum shape (more stats)             //
//                                                                        //
// 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_triggeronly(std::string horn   = "fhc")
{

  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);
  }


  double pot_sample      = kAna2017Period3POT;
  double livetime_sample = kAna2017Period3Livetime;
  // // in case of using period5
  // pot_sample = kAna2017Period5POT;
  // livetime_sample = kAna2017Period5Livetime;


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

  std::vector<TH1*> hcosmic;
  std::cout << "Getting histograms from files" << std::endl;
  for(int quantId = 0; quantId < 5; quantId++){
    hcosmic.push_back((TH1*)fcosmic->Get(Form("numuEOpt_q%d", quantId)));
  }
  

  std::cout << "Getting pot, livetime, count" << std::endl;  
  TH1* hlivetime_cosmic  = (TH1*)fcosmic->Get("livetime");
  float livetime_cosmic  = hlivetime_cosmic->Integral();
  // check counts agree in total
  // and get counts per quantile

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

  // note livetime scaling:
  // livetime sideband is not equal to livetime from beamspill
  float scale = livetime / livetime_cosmic;
  std::cout << "\ncosmic trigger livetime = " << livetime_cosmic << ", ana livetime = " << livetime << " -> scale = " << scale << std::endl;


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

  // 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++){
    hcosmic[quantId]->Scale( events[quantId] / (hcosmic[quantId]->Integral()) );
    if(quantId==0) std::cout << "Total integral = " << hcosmic[quantId]->Integral() << std::endl;
    else std::cout << "quantile " << quantId << " = " << hcosmic[quantId]->Integral() << std::endl;
  }
  
  std::cout << "\nWriting to output file" << std::endl;
  std::string name_out = "cosmicstriggeronly_" + horn + "__numu2018.root";
  TFile fout(name_out.c_str(),"recreate");
  for(int quantId=0; quantId<5; quantId++){  
    hcosmic[quantId]->Write(Form("cosmics_q%d", quantId));
  }
  
  fout.Close();

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

}