doReweight_dk2nu_original.C
Go to the documentation of this file.
1 
2 #include "dk2nu/tree/dkmeta.h"
3 #include "dk2nu/tree/dk2nu.h"
4 #include "MakeReweight.h"
5 
6 #include <string>
7 #include <vector>
8 #include <iostream>
9 #include <stdlib.h>
10 #include <iomanip>
11 
12 #include "TH1D.h"
13 
14 //Some constants::
15 
16 const int NbinsE = 120;
17 const double emin = 0.;
18 const double emax = 120.;
19 const int Nnuhel = 4;
20 const char* nuhel[Nnuhel] = {"numu","numubar","nue","nuebar"};
21 
22 using namespace NeutrinoFluxReweight;
23 
24 class nu_g4numi;
25 class Numi2Pdg;
26 
27 int idx_hel(int pdgdcode);
28 
29 /*!
30  * Run the reweighting for a single file (inputFile) for
31  * particular MIPP covariance matrix given in input.xml file.
32  */
33 void doReweight_dk2nu(const char* inputFile, const char* outputFile, const char* optionsFile){
34 
35  const char* thisDir = getenv("PPFX_DIR");
36  const char* OutputDir=thisDir;
37 
38  std::cout<< "Instance of MakeReweight()" <<std::endl;
39  MakeReweight* makerew = MakeReweight::getInstance();
40  makerew->SetOptions(optionsFile);
41 
42  std::cout<<"Making an output file to store histograms"<<std::endl;
43  TFile* fOut = new TFile(outputFile,"recreate");
44  std::cout<<"File name: "<<fOut->GetName()<<std::endl;
45 
46  fOut->mkdir("nom");
47  for(int ii=0;ii<Nnuhel;ii++){
48  fOut->mkdir(Form("%s_thintarget",nuhel[ii]));
49  fOut->mkdir(Form("%s_mippnumi",nuhel[ii]));
50  fOut->mkdir(Form("%s_attenuation",nuhel[ii]));
51  fOut->mkdir(Form("%s_others",nuhel[ii]));
52  fOut->mkdir(Form("%s_total",nuhel[ii]));
53  }
54  std::cout<<"Done making the output file"<<std::endl;
55 
56  int Nuniverses = makerew->GetNumberOfUniversesUsed();
57 
58  TH1D* hnom[Nnuhel];
59  TH1D* hcv[Nnuhel];
60  TH1D* hthin[Nnuhel][Nuniverses];
61  TH1D* hmipp[Nnuhel][Nuniverses];
62  TH1D* hatt[Nnuhel][Nuniverses];
63  TH1D* hothers[Nnuhel][Nuniverses];
64  TH1D* htotal[Nnuhel][Nuniverses];
65 
66  for(int ii=0;ii<Nnuhel;ii++){
67  hnom[ii] = new TH1D(Form("hnom_%s",nuhel[ii]),"",NbinsE,emin,emax);
68  hcv[ii] = new TH1D(Form("hcv_%s",nuhel[ii]),"",NbinsE,emin,emax);
69  for(int jj=0;jj< Nuniverses; jj++){
70  hthin[ii][jj] = new TH1D(Form("hthin_%s_%d", nuhel[ii],jj),"",NbinsE,emin,emax);
71  hmipp[ii][jj] = new TH1D(Form("hmipp_%s_%d", nuhel[ii],jj),"",NbinsE,emin,emax);
72  hatt[ii][jj] = new TH1D(Form("hatt_%s_%d", nuhel[ii],jj),"",NbinsE,emin,emax);
73  hothers[ii][jj] = new TH1D(Form("hothers_%s_%d",nuhel[ii],jj),"",NbinsE,emin,emax);
74  htotal[ii][jj] = new TH1D(Form("htotal_%s_%d", nuhel[ii],jj),"",NbinsE,emin,emax);
75  }
76  }
77 
78  //Loading ntuples:
79  TChain* chain_evts = new TChain("dk2nuTree");
80  TChain* chain_meta = new TChain("dkmetaTree");
81  bsim::Dk2Nu* dk2nu = new bsim::Dk2Nu;
82  bsim::DkMeta* dkmeta = new bsim::DkMeta;
83 
84  std::cout<<" Adding ntuple at: "<<inputFile<<std::endl;
85 
86  chain_evts->Add(inputFile);
87  chain_evts->SetBranchAddress("dk2nu",&dk2nu);
88  int nentries = chain_evts->GetEntries();
89 
90  chain_meta->Add(inputFile);
91  chain_meta->SetBranchAddress("dkmeta",&dkmeta);
92  chain_meta->GetEntry(0); //all entries are the same
93 
94  std::vector<double> vwgt_mipp_pi;
95  std::vector<double> vwgt_mipp_K;
96  std::vector<double> vwgt_abs;
97  std::vector<double> vwgt_att;
98  std::vector<double> vwgt_ttpCpi;
99  std::vector<double> vwgt_ttpCk;
100  std::vector<double> vwgt_ttnCpi;
101  std::vector<double> vwgt_ttpCnu;
102  std::vector<double> vwgt_ttnua;
103  std::vector<double> vwgt_ttmesinc;
104  std::vector<double> vwgt_oth;
105 
106  std::cout<<"N of entries: "<<nentries<<std::endl;
107 
108  for(int ii=0;ii<nentries;ii++){
109  if(ii%1000==0)std::cout<<ii/1000<<" k evts"<<std::endl;
110  vwgt_mipp_pi.clear();
111  vwgt_mipp_K.clear();
112  vwgt_abs.clear();
113  vwgt_att.clear();
114  vwgt_ttpCpi.clear();
115  vwgt_ttpCk.clear();
116  vwgt_ttnCpi.clear();
117  vwgt_ttpCnu.clear();
118  vwgt_ttmesinc.clear();
119  vwgt_ttnua.clear();
120  vwgt_oth.clear();
121 
122  chain_evts->GetEntry(ii);
123  makerew->calculateWeights(dk2nu,dkmeta);
124 
125  double fluxWGT = ( (dk2nu->nuray)[1].wgt )*(dk2nu->decay.nimpwt)/3.1416;
126  int nuidx = idx_hel(dk2nu->decay.ntype);
127  double nuenergy = (dk2nu->nuray)[1].E;
128 
129  if(nuidx<0){
130  std::cout<<"=> Wrong neutrino file"<<std::endl;
131  }
132  hnom[nuidx]->Fill(nuenergy,fluxWGT);
133  hcv[nuidx]->Fill(nuenergy,fluxWGT*makerew->GetCVWeight());
134 
135  vwgt_mipp_pi = makerew->GetWeights("MIPPNumiPionYields");
136  vwgt_mipp_K = makerew->GetWeights("MIPPNumiKaonYields");
137  vwgt_abs = makerew->GetWeights("TotalAbsorption");
138  vwgt_att = makerew->GetWeights("TargetAttenuation");
139  vwgt_ttpCpi = makerew->GetWeights("ThinTargetpCPion");
140  vwgt_ttpCk = makerew->GetWeights("ThinTargetpCKaon");
141  vwgt_ttnCpi = makerew->GetWeights("ThinTargetnCPion");
142  vwgt_ttpCnu = makerew->GetWeights("ThinTargetpCNucleon");
143  vwgt_ttmesinc= makerew->GetWeights("ThinTargetMesonIncident");
144  vwgt_ttnua = makerew->GetWeights("ThinTargetnucleonA");
145  vwgt_oth = makerew->GetWeights("Other");
146 
147 
148  for(int jj=0;jj<Nuniverses;jj++){
149  double wgt_thin = vwgt_ttpCpi[jj]*vwgt_ttpCk[jj]*vwgt_ttnCpi[jj]*vwgt_ttpCnu[jj]*vwgt_ttmesinc[jj]*vwgt_ttnua[jj];
150  double wgt_mipp = vwgt_mipp_pi[jj]*vwgt_mipp_K[jj];
151  double wgt_att = vwgt_att[jj]*vwgt_abs[jj];
152  hthin[nuidx][jj]->Fill(nuenergy,fluxWGT*wgt_thin);
153  hmipp[nuidx][jj]->Fill(nuenergy,fluxWGT*wgt_mipp);
154  hatt[nuidx][jj]->Fill(nuenergy,fluxWGT*wgt_att);
155  hothers[nuidx][jj]->Fill(nuenergy,fluxWGT*vwgt_oth[jj]);
156  htotal[nuidx][jj]->Fill(nuenergy,fluxWGT*wgt_thin*wgt_mipp*wgt_att*vwgt_oth[jj]);
157 
158  }
159  }
160 
161  std::cout<<"storing general histos"<<std::endl;
162  fOut->cd();
163  for(int ii=0;ii<Nnuhel;ii++){
164  fOut->cd("nom");
165  hnom[ii]->Write();
166  hcv[ii]->Write();
167  for(int jj=0;jj< Nuniverses; jj++){
168  fOut->cd(Form("%s_thintarget",nuhel[ii])); hthin[ii][jj]->Write();
169  fOut->cd(Form("%s_mippnumi",nuhel[ii])); hmipp[ii][jj]->Write();
170  fOut->cd(Form("%s_attenuation",nuhel[ii])); hatt[ii][jj]->Write();
171  fOut->cd(Form("%s_others",nuhel[ii])); hothers[ii][jj]->Write();
172  fOut->cd(Form("%s_total",nuhel[ii])); htotal[ii][jj]->Write();
173  }
174  }
175 
176  //Releasing memory:
177  makerew->resetInstance();
178 
179  std::cout<<"End of run()"<<std::endl;
180 
181 }
182 
183 int idx_hel(int pdgcode){
184  int idx = -1;
185  if(pdgcode == 14)idx = 0;
186  if(pdgcode == -14)idx = 1;
187  if(pdgcode == 12)idx = 2;
188  if(pdgcode == -12)idx = 3;
189  return idx;
190 }
191 
192 //////////////
193 
194 #ifndef __CINT__
195 int main(int argc, const char* argv[]){
196 
197  doReweight_dk2nu(argv[1],argv[2],argv[3]);
198  return 0;
199 }
200 #endif
nuhel
const char * nuhel[Nnuhel]
Definition: doReweight_dk2nu_original.C:20
idx_hel
int idx_hel(int pdgdcode)
Definition: doReweight_dk2nu_original.C:183
NeutrinoFluxReweight::MakeReweight::GetWeights
std::vector< double > GetWeights(std::string nameReweighter)
get the vector of the weights for a given reweighter
Definition: MakeReweight.cpp:137
NeutrinoFluxReweight::MakeReweight
A class to make the reweight event by event.
Definition: MakeReweight.h:32
NeutrinoFluxReweight::MakeReweight::GetCVWeight
double GetCVWeight()
get the cv weights
Definition: MakeReweight.cpp:148
nu_g4numi
Definition: nu_g4numi.h:9
NeutrinoFluxReweight::MakeReweight::getInstance
static MakeReweight * getInstance()
Definition: MakeReweight.cpp:159
doReweight_dk2nu
void doReweight_dk2nu(const char *inputFile, const char *outputFile, const char *optionsFile)
Definition: doReweight_dk2nu_original.C:33
emax
const double emax
Definition: doReweight_dk2nu_original.C:18
cout
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
Definition: MessageFacilityAsStandAlonePackage.txt:82
main
int main(int argc, const char *argv[])
Definition: doReweight_dk2nu_original.C:195
E
const Double_t E[nknots]
Definition: testXsec.C:25
cet::getenv
std::string getenv(std::string const &name)
Definition: getenv.cc:15
NeutrinoFluxReweight::MakeReweight::GetNumberOfUniversesUsed
int GetNumberOfUniversesUsed()
number of universes used in this run
Definition: MakeReweight.cpp:152
NeutrinoFluxReweight::MakeReweight::calculateWeights
void calculateWeights(nu_g4numi *nu, const char *tgtcfg, const char *horncfg)
Definition: MakeReweight.cpp:94
emin
const double emin
Definition: doReweight_dk2nu_original.C:17
Nnuhel
const int Nnuhel
Definition: doReweight_dk2nu_original.C:19
NbinsE
const int NbinsE
Definition: doReweight_dk2nu_original.C:16
NeutrinoFluxReweight::MakeReweight::SetOptions
void SetOptions(std::string fileIn)
Definition: MakeReweight.cpp:11
Generated by   doxygen 1.8.11