make_nus17_pot.C
Go to the documentation of this file.
1 // Calculate analysis exposure, selected events as a function of run
2 
3 #include <bitset>
4 #include <cstdlib>
5 #include <iostream>
6 #include "NuXAna/macros/Nus17/Nus17Exposures.h"
7 #include "CAFAna/Cuts/Cuts.h"
9 #include "CAFAna/Cuts/NueCutsFirstAna.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Core/Loaders.h"
12 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Vars/HistAxes.h"
17 #include "CAFAna/Core/EventList.h"
20 using namespace ana;
21 
22 #include "TFile.h"
23 #include "TH1D.h"
24 
25 //fiducial mass accounting from docdb-13556
26 
28  {
29  return spill->spilltimesec;
30  });
31 const SpillVar kRun( [](const caf::SRSpillProxy* spill)
32  {
33  return spill->run;
34  });
35 const SpillVar fPOT( [](const caf::SRSpillProxy* spill)
36  {
37  return spill->spillpot;
38  });
39 
40 
41 const SpillVar kLiveMass(
42  [](const caf::SRSpillProxy* spill)
43  {
44  std::string binary = std::bitset<14>(spill->dibmask).to_string();
45  double live = 12e-6;
46  if (spill->run < 17946) live = 24e-6;
47  double sum = 0;
48  int dbcount = 1;
49  int mydb = 0;
50  for (int i=13; i>-1; --i){
51  int test = binary[i] - '0';
52  if (dbcount >= spill->dibfirst &&
53  dbcount <= spill->diblast &&
54  test == 1) mydb++;
55  if (dbcount >= spill->dibfirst &&
56  dbcount <= spill->diblast &&
57  test == 0){
58  if (mydb >= 4) sum+= live*kNus17FiducialMass[mydb];
59  mydb = 0;
60  }
61  dbcount++;
62  }
63  if (mydb >= 4) sum += live*kNus17FiducialMass[mydb];
64  return sum;
65  });
66 const SpillVar kLive( [](const caf::SRSpillProxy* spill)
67  {
68  if (spill->run < 17946)
69  return 24e-6;
70  else
71  return 12e-6;
72  });
73 const SpillVar fPOTMass(
74  [](const caf::SRSpillProxy* spill)
75  {
76  std::string binary = std::bitset<14>(spill->dibmask).to_string();
77  double pot = 0;
78  int dbcount = 1;
79  int mydb = 0;
80  for (int i=13; i>-1; --i){
81  int test = binary[i] - '0';
82  if (dbcount >= spill->dibfirst &&
83  dbcount <= spill->diblast &&
84  test == 1) mydb++;
85  if (dbcount >= spill->dibfirst &&
86  dbcount <= spill->diblast &&
87  test == 0){
88  if (mydb >= 4) pot+= spill->spillpot*kNus17FiducialMass[mydb];
89  mydb = 0;
90  }
91  dbcount++;
92  }
93  if (mydb >= 4) pot += spill->spillpot*kNus17FiducialMass[mydb];
94  return pot;
95  });
96 
97 const Var kSliceTime(
98  [](const caf::SRProxy* sr)
99  {
100  return sr->hdr.unixtime;
101  });
102 
103 
105 {
106  TH1::SetDefaultSumw2();
107 
108  //binned in unix time encompassing first analysis period
109  TH1D* hLiveMass = new TH1D("hLiveMass","livetime exposure;run;exposure (s #times kT)",12473, 12939.5, 25412.5);
110  TH1D* hLive = new TH1D("hLive","livetime;run;livetime (s)",12473, 12939.5, 25412.5);
111  TH1D* hPOTMass = new TH1D("hPOTMass","exposure;run;exposure (POT #times kT)",12473, 12939.5, 25412.5);
112  TH1D* hPOT = new TH1D("hPOT","POT;run;POT",12473, 12939.5, 25412.5);
113  TH1D* hLiveMassT = new TH1D("hLiveMassT","livetime exposure;run;exposure (s #times kT)",9900,1391000000,1490000000);
114  TH1D* hLiveT = new TH1D("hLiveT","livetime;run;livetime (s)",9900,1391000000,1490000000);
115  TH1D* hPOTMassT = new TH1D("hPOTMassT","exposure;run;exposure (POT #times kT)",9900,1391000000,1490000000);
116  TH1D* hPOTT = new TH1D("hPOTT","POT;run;POT",9900,1391000000,1490000000);
117 
118  const std::vector<std::string> fNumiUnblind(MakeNumiUnblindList());
119 
120  std::vector<std::string> fPeriod;
121  if(period == 0) fPeriod = fNumiUnblind;
122  else if(period == 1) fPeriod = fnumi_period1;
123  else if(period == 2) fPeriod = fnumi_period2;
124  else if(period == 3) fPeriod = fnumi_period3;
125  else if(period == 5) fPeriod = fnumi_period5;
126  else {
127  std::cout << "Invalid period. Select a number from 0-5, excl. 4" << std::endl;
128  std::cout << "0 selects all periods." << std::endl;
129  abort();
130  }
131  //load data files
132  SpectrumLoader loaderFDData(fPeriod, kCosmic);
133 
134  loaderFDData.SetSpillCut(kStandardSpillCuts);
135  loaderFDData.AddSpillHistogram(hLiveMass, kRun, kStandardSpillCuts,
136  kLiveMass);
137  loaderFDData.AddSpillHistogram(hLive, kRun, kStandardSpillCuts,
138  kLive);
139  loaderFDData.AddSpillHistogram(hPOTMass, kRun, kStandardSpillCuts,
140  fPOTMass);
141  loaderFDData.AddSpillHistogram(hPOT, kRun, kStandardSpillCuts,
142  fPOT);
143  loaderFDData.AddSpillHistogram(hLiveMassT, kSpillTime, kStandardSpillCuts,
144  kLiveMass);
145  loaderFDData.AddSpillHistogram(hLiveT, kSpillTime, kStandardSpillCuts,
146  kLive);
147  loaderFDData.AddSpillHistogram(hPOTMassT, kSpillTime, kStandardSpillCuts,
148  fPOTMass);
149  loaderFDData.AddSpillHistogram(hPOTT, kSpillTime, kStandardSpillCuts,
150  fPOT);
151  loaderFDData.Go();
152 
153 
154  std::cout << "********* Period " << std::to_string(period)
155  << " *********" << std::endl;
156 
157  std::cout << "Live Integral: " << hLive->Integral() << std::endl;
158  std::cout << "LiveT Integral: " << hLiveT->Integral() << std::endl;
159  std::cout << "LiveMass Integral: " << hLiveMass->Integral() << std::endl;
160  std::cout << "LiveMassT Integral: " << hLiveMassT->Integral() << std::endl;
161  std::cout << "POT Integral: " << hPOT->Integral() << std::endl;
162  std::cout << "POTT Integral: " << hPOTT->Integral() << std::endl;
163  std::cout << "POTMass Integral: " << hPOTMass->Integral() << std::endl;
164  std::cout << "POTMassT Integral: " << hPOTMassT->Integral() << std::endl;
165 
166 
167  std::string outname = "nus17_pot_period"+std::to_string(period)+".root";
168 
169  TFile* outFile = new TFile(outname.c_str(), "RECREATE");
170  hLiveMass->Write();
171  hLive->Write();
172  hPOTMass->Write();
173  hPOT->Write();
174  hLiveMassT->Write();
175  hLiveT->Write();
176  hPOTMassT->Write();
177  hPOTT->Write();
178  outFile->Close();
179 
180 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< std::string > fnumi_period3
Definition: NusLoadProd3.h:77
const SpillVar kLiveMass([](const caf::SRSpillProxy *spill){std::string binary=std::bitset< 14 >(spill->dibmask).to_string();double live=12e-6;if(spill->run< 17946) live=24e-6;double sum=0;int dbcount=1;int mydb=0;for(int i=13;i >-1;--i){int test=binary[i]- '0';if(dbcount >=spill->dibfirst && dbcount<=spill->diblast && test==1) mydb++;if(dbcount >=spill->dibfirst && dbcount<=spill->diblast && test==0){if(mydb >=4) sum+=live *kNus17FiducialMass[mydb];mydb=0;}dbcount++;}if(mydb >=4) sum+=live *kNus17FiducialMass[mydb];return sum;})
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
caf::Proxy< float > spillpot
Definition: SRProxy.h:1407
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< short unsigned int > dibfirst
Definition: SRProxy.h:1362
const SpillVar kLive([](const caf::SRSpillProxy *spill){if(spill->run< 17946) return 24e-6;else return 12e-6;})
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:34
std::vector< std::string > MakeNumiUnblindList()
Definition: NusLoadProd3.h:84
void SetSpillCut(const SpillCut &cut)
std::vector< std::string > fnumi_period1
Definition: NusLoadProd3.h:71
const SpillVar kSpillTime([](const caf::SRSpillProxy *spill){return spill->spilltimesec;})
TFile * outFile
Definition: PlotXSec.C:135
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
const SpillVar fPOTMass([](const caf::SRSpillProxy *spill){std::string binary=std::bitset< 14 >(spill->dibmask).to_string();double pot=0;int dbcount=1;int mydb=0;for(int i=13;i >-1;--i){int test=binary[i]- '0';if(dbcount >=spill->dibfirst && dbcount<=spill->diblast && test==1) mydb++;if(dbcount >=spill->dibfirst && dbcount<=spill->diblast && test==0){if(mydb >=4) pot+=spill->spillpot *kNus17FiducialMass[mydb];mydb=0;}dbcount++;}if(mydb >=4) pot+=spill->spillpot *kNus17FiducialMass[mydb];return pot;})
caf::Proxy< long unsigned int > spilltimesec
Definition: SRProxy.h:1409
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void make_nus17_pot(int period)
std::vector< std::string > fnumi_period2
Definition: NusLoadProd3.h:74
Proxy for caf::SRSpill.
Definition: SRProxy.h:1346
OStream cout
Definition: OStream.cxx:6
caf::Proxy< unsigned int > run
Definition: SRProxy.h:1406
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::vector< std::string > fnumi_period5
Definition: NusLoadProd3.h:80
caf::Proxy< short unsigned int > dibmask
Definition: SRProxy.h:1364
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< float > unixtime
Definition: SRProxy.h:255
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Var kRun
Definition: Vars.cxx:20
Template for Var and SpillVar.
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
const SpillVar fPOT([](const caf::SRSpillProxy *spill){return spill->spillpot;})
enum BeamMode string