Belle II Software  release-06-01-15
NtuplePhase1_v6Module.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 // Own include
10 #include <beast/analysis/modules/NtuplePhase1_v6Module.h>
11 #include <beast/analysis/modules/BEASTTree_v5.h>
12 
13 // framework - DataStore
14 #include <framework/datastore/DataStore.h>
15 #include <framework/datastore/StoreObjPtr.h>
16 
17 // framework aux
18 #include <framework/logging/Logger.h>
19 #include <framework/core/RandomNumbers.h>
20 
21 // DataStore classes
22 #include <framework/io/RootIOUtilities.h>
23 #include <framework/dataobjects/EventMetaData.h>
24 
25 #include <TMath.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TString.h>
29 
30 using namespace std;
31 
32 namespace Belle2 {
38  using namespace RootIOUtilities;
39 
40  //-----------------------------------------------------------------
41  // Register module
42  //-----------------------------------------------------------------
43 
44  REG_MODULE(NtuplePhase1_v6)
45 
46  //-----------------------------------------------------------------
47  // Implementation
48  //-----------------------------------------------------------------
49 
50  NtuplePhase1_v6Module::NtuplePhase1_v6Module() : Module()
51  {
52  // set module description (e.g. insert text)
53  setDescription("Read SKB PVs, simulated measurements of BEAST sensors, and write scaled simulated Ntuple in BEAST phase 1 data format");
54 
55 
56  // Add parameters
57  addParam("inputFileNames", m_inputFileNames,
58  "List of files with SKB PVs ");
59 
60  addParam("outputFileName", m_outputFileName, "Output file name");
61 
62  //addParam("input_ts", m_input_ts, "Input time stamp start and stop");
63 
64  addParam("input_Time_eqv", m_input_Time_eqv, "time-eqv");
65 
66  addParam("input_I_HER", m_input_I_HER, "HER current");
67  addParam("input_I_LER", m_input_I_LER, "LER current");
68 
69  addParam("input_P_HER", m_input_P_HER, "HER pressure");
70  addParam("input_P_LER", m_input_P_LER, "LER pressure");
71 
72  addParam("input_sigma_x_HER", m_input_sigma_x_HER, "HER beam size");
73  addParam("input_sigma_x_LER", m_input_sigma_x_LER, "LER beam size");
74  addParam("input_sigma_y_HER", m_input_sigma_y_HER, "HER beam size");
75  addParam("input_sigma_y_LER", m_input_sigma_y_LER, "LER beam size");
76 
77  addParam("input_bunchNb_HER", m_input_bunchNb_HER, "HER bunch number");
78  addParam("input_bunchNb_LER", m_input_bunchNb_LER, "LER bunch number");
79 
80  addParam("input_data_bunchNb_HER", m_input_data_bunchNb_HER, "HER bunch number");
81  addParam("input_data_bunchNb_LER", m_input_data_bunchNb_LER, "LER bunch number");
82  addParam("input_data_SingleBeam", m_input_data_SingleBeam, "LER/HER/Both");
83 
84  addParam("input_LT_SAD_RLR", m_input_LT_SAD_RLR, "SAD Low Ring Loss Rate");
85  addParam("input_HT_SAD_RLR", m_input_HT_SAD_RLR, "SAD High Ring Loss Rate");
86  addParam("input_LC_SAD_RLR", m_input_LC_SAD_RLR, "SAD Low Ring Loss Rate");
87  addParam("input_HC_SAD_RLR", m_input_HC_SAD_RLR, "SAD High Ring Loss Rate");
88  addParam("input_LB_SAD_RLR", m_input_LB_SAD_RLR, "SAD Low Ring Loss Rate");
89  addParam("input_HB_SAD_RLR", m_input_HB_SAD_RLR, "SAD High Ring Loss Rate");
90  addParam("input_LC_SAD_RLR_av", m_input_LC_SAD_RLR_av, "SAD Low Ring Loss Rate");
91  addParam("input_HC_SAD_RLR_av", m_input_HC_SAD_RLR_av, "SAD High Ring Loss Rate");
92  addParam("input_LB_SAD_RLR_av", m_input_LB_SAD_RLR_av, "SAD Low Ring Loss Rate");
93  addParam("input_HB_SAD_RLR_av", m_input_HB_SAD_RLR_av, "SAD High Ring Loss Rate");
94 
95  addParam("input_BGSol", m_input_BGSol, "BG solution 0 or 1");
96  addParam("input_ToSol", m_input_ToSol, "To solution 0 or 1");
97 
98  addParam("input_Z", m_input_Z, "Z number");
99 
100  addParam("input_GasCorrection", m_input_GasCorrection, "GasCorrection");
101 
102  addParam("input_part", m_input_part, "Part");
103 
104  addParam("inputHistoFileNames", m_inputHistoFileNames,
105  "List of root files with histograms");
106 
107  addParam("inputRateHistoNames", m_inputRateHistoNames,
108  "List of rate histograms");
109  addParam("inputDoseHistoNames", m_inputDoseHistoNames,
110  "List of dose histograms");
111  addParam("inputRateHistoNamesVrs", m_inputRateHistoNamesVrs,
112  "List of rate histograms");
113  addParam("inputDoseHistoNamesVrs", m_inputDoseHistoNamesVrs,
114  "List of dose histograms");
115 
116 
117  addParam("input_PIN_width", m_input_PIN_width, "PIN width");
118  addParam("input_HE3_EfCor", m_input_HE3_EfCor, "HE3 inefficiency correction");
119 
120  // initialize other private data members
121  m_file = NULL;
122  m_tree = NULL;
123  m_treeBEAST = NULL;
124  m_beast.clear();
125 
126  m_numEntries = 0;
127  m_entryCounter = 0;
128 
129  m_DayBin = 0;
130  }
131 
132  NtuplePhase1_v6Module::~NtuplePhase1_v6Module()
133  {
134  }
135 
136  void NtuplePhase1_v6Module::initialize()
137  {
138  // read TFile with histograms
139 
140  // expand possible wildcards
141  m_inputFileNames = expandWordExpansions(m_inputFileNames);
142  if (m_inputFileNames.empty()) {
143  B2FATAL("No valid files specified!");
144  }
145 
146  // check files
147  TDirectory* dir = gDirectory;
148  for (const string& fileName : m_inputFileNames) {
149  TFile* f = TFile::Open(fileName.c_str(), "READ");
150  if (!f or !f->IsOpen()) {
151  B2FATAL("Couldn't open input file " + fileName);
152  }
153  delete f;
154  }
155  dir->cd();
156 
157  // get event TTree
158  //m_tree = new TChain(c_treeNames[DataStore::c_Event].c_str());
159  m_tree = new TChain("tout");
160  for (const string& fileName : m_inputFileNames) {
161  m_tree->AddFile(fileName.c_str());
162  }
163  m_numEvents = m_tree->GetEntries();
164  if (m_numEvents == 0) B2ERROR(c_treeNames[DataStore::c_Event] << " has no entires");
165  m_eventCount = 0;
166 
167  m_tree->SetBranchAddress("ts", &(m_beast.ts));
168  m_tree->SetBranchAddress("event", &(m_beast.event));
169  m_tree->SetBranchAddress("run", &(m_beast.run));
170  m_tree->SetBranchAddress("subrun", &(m_beast.subrun));
171  m_tree->SetBranchAddress("SKB_HER_injectionFlag", &(m_beast.SKB_HER_injectionFlag));
172  m_tree->SetBranchAddress("SKB_LER_injectionFlag", &(m_beast.SKB_LER_injectionFlag));
173  m_tree->SetBranchAddress("SKB_HER_injectionFlag_safe", &(m_beast.SKB_HER_injectionFlag_safe));
174  m_tree->SetBranchAddress("SKB_LER_injectionFlag_safe", &(m_beast.SKB_LER_injectionFlag_safe));
175  m_tree->SetBranchAddress("SKB_HER_abortFlag", &(m_beast.SKB_HER_abortFlag));
176  m_tree->SetBranchAddress("SKB_LER_abortFlag", &(m_beast.SKB_LER_abortFlag));
177  m_tree->SetBranchAddress("SKB_HER_abortFlag_safe", &(m_beast.SKB_HER_abortFlag_safe));
178  m_tree->SetBranchAddress("SKB_LER_abortFlag_safe", &(m_beast.SKB_LER_abortFlag_safe));
179  m_tree->SetBranchAddress("SKB_Status", &(m_beast.SKB_Status));
180  m_tree->SetBranchAddress("SKB_HER_injectionRate", &(m_beast.SKB_HER_injectionRate));
181  m_tree->SetBranchAddress("SKB_LER_injectionRate", &(m_beast.SKB_LER_injectionRate));
182  m_tree->SetBranchAddress("SKB_HER_lifetime", &(m_beast.SKB_HER_lifetime));
183  m_tree->SetBranchAddress("SKB_LER_lifetime", &(m_beast.SKB_LER_lifetime));
184  m_tree->SetBranchAddress("SKB_LER_current", &(m_beast.SKB_LER_current));
185  m_tree->SetBranchAddress("SKB_HER_current", &(m_beast.SKB_HER_current));
186  m_tree->SetBranchAddress("SKB_LER_injectionEfficiency", &(m_beast.SKB_LER_injectionEfficiency));
187  m_tree->SetBranchAddress("SKB_HER_injectionEfficiency", &(m_beast.SKB_HER_injectionEfficiency));
188  m_tree->SetBranchAddress("SKB_beamLoss_ionChambers_mean", &(m_beast.SKB_beamLoss_ionChambers_mean));
189  m_tree->SetBranchAddress("SKB_beamLoss_PINdiodes_mean", &(m_beast.SKB_beamLoss_PINdiodes_mean));
190  m_tree->SetBranchAddress("SKB_beamLoss_nearCollimators", &(m_beast.SKB_beamLoss_nearCollimators));
191  m_tree->SetBranchAddress("SKB_beamLoss_aroundMasks", &(m_beast.SKB_beamLoss_aroundMasks));
192  m_tree->SetBranchAddress("SKB_LER_injectionCharge", &(m_beast.SKB_LER_injectionCharge));
193  m_tree->SetBranchAddress("SKB_HER_injectionCharge", &(m_beast.SKB_HER_injectionCharge));
194  m_tree->SetBranchAddress("SKB_LER_injectionRepetitionRate", &(m_beast.SKB_LER_injectionRepetitionRate));
195  m_tree->SetBranchAddress("SKB_HER_injectionRepetitionRate", &(m_beast.SKB_HER_injectionRepetitionRate));
196  m_tree->SetBranchAddress("SKB_LER_injectionNumberOfBunches", &(m_beast.SKB_LER_injectionNumberOfBunches));
197  m_tree->SetBranchAddress("SKB_HER_injectionNumberOfBunches", &(m_beast.SKB_HER_injectionNumberOfBunches));
198  m_tree->SetBranchAddress("SKB_LER_pressures", &(m_beast.SKB_LER_pressures));
199  m_tree->SetBranchAddress("SKB_HER_pressures", &(m_beast.SKB_HER_pressures));
200  m_tree->SetBranchAddress("SKB_LER_pressure_average", &(m_beast.SKB_LER_pressure_average));
201  m_tree->SetBranchAddress("SKB_HER_pressure_average", &(m_beast.SKB_HER_pressure_average));
202  m_tree->SetBranchAddress("SKB_LER_pressures_corrected", &(m_beast.SKB_LER_pressures_corrected));
203  m_tree->SetBranchAddress("SKB_HER_pressures_corrected", &(m_beast.SKB_HER_pressures_corrected));
204  m_tree->SetBranchAddress("SKB_LER_pressure_average_corrected", &(m_beast.SKB_LER_pressure_average_corrected));
205  m_tree->SetBranchAddress("SKB_HER_pressure_average_corrected", &(m_beast.SKB_HER_pressure_average_corrected));
206  m_tree->SetBranchAddress("SKB_HER_collimatorPositions_mm", &(m_beast.SKB_HER_collimatorPositions_mm));
207  m_tree->SetBranchAddress("SKB_HER_collimatorPositions_DMM", &(m_beast.SKB_HER_collimatorPositions_DMM));
208  m_tree->SetBranchAddress("SKB_HER_collimatorPositions_inX", &(m_beast.SKB_HER_collimatorPositions_inX));
209  m_tree->SetBranchAddress("SKB_HER_collimatorPositions_inY", &(m_beast.SKB_HER_collimatorPositions_inY));
210  m_tree->SetBranchAddress("SKB_HER_collimatorPositions_fromBeam", &(m_beast.SKB_HER_collimatorPositions_fromBeam));
211  m_tree->SetBranchAddress("SKB_LER_collimatorPositions_mm", &(m_beast.SKB_LER_collimatorPositions_mm));
212  m_tree->SetBranchAddress("SKB_LER_collimatorPositions_X", &(m_beast.SKB_LER_collimatorPositions_X));
213  m_tree->SetBranchAddress("SKB_LER_collimatorPositions_Y", &(m_beast.SKB_LER_collimatorPositions_Y));
214  m_tree->SetBranchAddress("SKB_LER_collimatorPositions_fromBeam", &(m_beast.SKB_LER_collimatorPositions_fromBeam));
215  m_tree->SetBranchAddress("SKB_HER_beamSize_xray_X", &(m_beast.SKB_HER_beamSize_xray_X));
216  m_tree->SetBranchAddress("SKB_HER_beamSize_xray_Y", &(m_beast.SKB_HER_beamSize_xray_Y));
217  m_tree->SetBranchAddress("SKB_HER_correctedBeamSize_xray_Y", &(m_beast.SKB_HER_correctedBeamSize_xray_Y));
218  m_tree->SetBranchAddress("SKB_LER_beamSize_xray_X", &(m_beast.SKB_LER_beamSize_xray_X));
219  m_tree->SetBranchAddress("SKB_LER_beamSize_xray_Y", &(m_beast.SKB_LER_beamSize_xray_Y));
220  m_tree->SetBranchAddress("SKB_LER_correctedBeamSize_xray_Y", &(m_beast.SKB_LER_correctedBeamSize_xray_Y));
221  m_tree->SetBranchAddress("SKB_LER_beamSize_SR_X", &(m_beast.SKB_LER_beamSize_SR_X));
222  m_tree->SetBranchAddress("SKB_LER_beamSize_SR_Y", &(m_beast.SKB_LER_beamSize_SR_Y));
223  m_tree->SetBranchAddress("SKB_HER_beamSize_SR_X", &(m_beast.SKB_HER_beamSize_SR_X));
224  m_tree->SetBranchAddress("SKB_HER_beamSize_SR_Y", &(m_beast.SKB_HER_beamSize_SR_Y));
225  m_tree->SetBranchAddress("SKB_HER_integratedCurrent", &(m_beast.SKB_HER_integratedCurrent));
226  m_tree->SetBranchAddress("SKB_LER_integratedCurrent", &(m_beast.SKB_LER_integratedCurrent));
227  m_tree->SetBranchAddress("SKB_LER_partialPressures_D06", &(m_beast.SKB_LER_partialPressures_D06));
228  m_tree->SetBranchAddress("SKB_LER_partialPressures_D02", &(m_beast.SKB_LER_partialPressures_D02));
229  m_tree->SetBranchAddress("SKB_LER_pressures_local", &(m_beast.SKB_LER_pressures_local));
230  m_tree->SetBranchAddress("SKB_HER_pressures_local", &(m_beast.SKB_HER_pressures_local));
231  m_tree->SetBranchAddress("SKB_LER_pressures_local_corrected", &(m_beast.SKB_LER_pressures_local_corrected));
232  m_tree->SetBranchAddress("SKB_HER_pressures_local_corrected", &(m_beast.SKB_HER_pressures_local_corrected));
233  m_tree->SetBranchAddress("SKB_LER_Zeff_D02", &(m_beast.SKB_LER_Zeff_D02));
234  m_tree->SetBranchAddress("SKB_LER_Zeff_D06", &(m_beast.SKB_LER_Zeff_D06));
235  m_tree->SetBranchAddress("CSI_sumE", &(m_beast.CSI_data_sumE));
236  m_tree->SetBranchAddress("BGO_dose", &(m_beast.BGO_data_dose));
237  m_tree->SetBranchAddress("PIN_dose", &(m_beast.PIN_data_dose));
238  m_tree->SetBranchAddress("DIA_dose", &(m_beast.DIA_data_dose));
239 
240  m_tree->SetBranchAddress("HE3_rate", &(m_beast.HE3_data_rate));
241  m_tree->SetBranchAddress("CSI_hitRate", &(m_beast.CSI_data_rate));
242  m_tree->SetBranchAddress("CSI_binnedE", &(m_beast.CSI_data_Ebin));
243  m_tree->SetBranchAddress("SCI_rate", &(m_beast.QCSS_data_rate));
244  m_tree->SetBranchAddress("CLW_N_MIPs_online", &(m_beast.CLAWS_data_rate));
245 
246  if (m_numEvents > 0) {
247  m_tree->GetEntry(0);
248  m_DayBin = (int)((m_beast.ts - 1454943600) / 60. / 60. / 24.);
249  }
250 
251  // expand possible wildcards
252  m_inputHistoFileNames = expandWordExpansions(m_inputHistoFileNames);
253  if (m_inputFileNames.empty()) {
254  B2FATAL("No valid files specified!");
255  }
256 
257  fctRate_HB = new TF1("fctRate_HB", "[0] * x*x * log([1] / TMath::Power(x,1./3.) + [2])", 1.0, 19.0);
258  fctRate_LB = new TF1("fctRate_LB", "[0] * x*x * log([1] / TMath::Power(x,1./3.) + [2])", 1.0, 19.0);
259  fctRate_HC = new TF1("fctRate_HC", "[0] * x*x / TMath::Power( ([1] / TMath::Power(x,1./3.) + [2]), 2.)", 1.0, 19.0);
260  fctRate_LC = new TF1("fctRate_LC", "[0] * x*x / TMath::Power( ([1] / TMath::Power(x,1./3.) + [2]), 2.)", 1.0, 19.0);
261  fctRate_HB->SetParameters(0.183373, 0.117173, 1.23431);
262  fctRate_LB->SetParameters(0.900838, 0.0455552, 1.10098);
263  fctRate_HC->SetParameters(1.80992, -0.000115401, 8.4047);
264  fctRate_LC->SetParameters(0.210872, -4.50637e-06, 1.64209);
265  m_input_Z_scaling[0] = fctRate_HC->Eval(m_input_Z[0]) / fctRate_HC->Eval(7);
266  m_input_Z_scaling[1] = fctRate_LC->Eval(m_input_Z[1]) / fctRate_LC->Eval(7);
267  m_input_Z_scaling[2] = fctRate_HB->Eval(m_input_Z[2]) / fctRate_HB->Eval(7);
268  m_input_Z_scaling[3] = fctRate_LB->Eval(m_input_Z[3]) / fctRate_LB->Eval(7);
269 
270  if (m_input_Z[0] == 0) m_input_Z_scaling[0] = 0;
271  if (m_input_Z[1] == 0) m_input_Z_scaling[1] = 0;
272  if (m_input_Z[2] == 0) m_input_Z_scaling[2] = 0;
273  if (m_input_Z[3] == 0) m_input_Z_scaling[3] = 0;
274 
275  double volume;
276  double rho;
277  double mass = 0.;
278  const double RadConv = 6.24e7; // 1 mrad = 6.24e7 MeV/kg
279 
280  // check files
281  TDirectory* dirh = gDirectory;
282  TFile* fh[6];
283  int iter = 0;
284  for (const string& fileNameTmp : m_inputHistoFileNames) {
285  const TString& fileName = (TString)fileNameTmp;
286  fh[iter] = TFile::Open(fileName, "READ");
287  if (!fh[iter] or !fh[iter]->IsOpen()) {
288  B2FATAL("Couldn't open input file " + fileName);
289  }
290  if (fileName.Contains("Touschek") || fileName.Contains("Coulomb") || fileName.Contains("Brems")) {
291  for (const string& HistoRateNameTmp : m_inputRateHistoNames) {
292  const TString& HistoRateName = (TString)HistoRateNameTmp;
293 
294  TH1F* h1D;
295  if (HistoRateName.Contains("csi")) h1D = (TH1F*)fh[iter]->Get(TString::Format("csi_rate_%d", m_DayBin));
296  else h1D = (TH1F*)fh[iter]->Get(HistoRateName);
297 
298  for (int i = 0; i < h1D->GetNbinsX(); i++) {
299  double counts = h1D->GetBinContent(i + 1);
300  double rate = counts / m_input_Time_eqv;
301 
302  if (fileName.Contains("Coulomb")) {
303  if (fileName.Contains("HER")) rate *= m_input_Z_scaling[0];
304  if (fileName.Contains("LER")) rate *= m_input_Z_scaling[1];
305  }
306  if (fileName.Contains("Brems")) {
307  if (fileName.Contains("HER")) rate *= m_input_Z_scaling[2];
308  if (fileName.Contains("LER")) rate *= m_input_Z_scaling[3];
309  }
310 
311  if (HistoRateName.Contains("Def")) rate *= m_input_HE3_EfCor[i];
312 
313  if (fileName.Contains("HER")) {
314  if (HistoRateName.Contains("qcss") && fileName.Contains("Touschek")) m_input_HT_QCSS_rate.push_back(rate); //Hz
315  if (HistoRateName.Contains("claws") && fileName.Contains("Touschek")) m_input_HT_CLAWS_rate.push_back(rate); //Hz
316  if (HistoRateName.Contains("csi") && fileName.Contains("Touschek")) m_input_HT_CSI_rate.push_back(rate); //Hz
317  if (HistoRateName.Contains("Def") && fileName.Contains("Touschek")) m_input_HT_HE3_rate.push_back(rate); //Hz
318  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Touschek")) m_input_HT_TPC_rate.push_back(rate); //Hz
319  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Touschek")) m_input_HT_TPC_angular_rate.push_back(rate); //Hz
320  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Touschek")) m_input_HT_TPC_angular_dose.push_back(rate); //Hz
321  if (HistoRateName.Contains("qcss") && fileName.Contains("Brems")) m_input_HB_QCSS_rate_av.push_back(rate); //Hz
322  if (HistoRateName.Contains("claws") && fileName.Contains("Brems")) m_input_HB_CLAWS_rate_av.push_back(rate); //Hz
323  if (HistoRateName.Contains("csi") && fileName.Contains("Brems")) m_input_HB_CSI_rate_av.push_back(rate); //Hz
324  if (HistoRateName.Contains("Def") && fileName.Contains("Brems")) m_input_HB_HE3_rate_av.push_back(rate); //Hz
325  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Brems")) m_input_HB_TPC_rate_av.push_back(rate); //Hz
326  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Brems")) m_input_HB_TPC_angular_rate_av.push_back(rate); //Hz
327  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Brems")) m_input_HB_TPC_angular_dose_av.push_back(rate); //Hz
328  if (HistoRateName.Contains("qcss") && fileName.Contains("Coulomb")) m_input_HC_QCSS_rate_av.push_back(rate); //Hz
329  if (HistoRateName.Contains("claws") && fileName.Contains("Coulomb")) m_input_HC_CLAWS_rate_av.push_back(rate); //Hz
330  if (HistoRateName.Contains("csi") && fileName.Contains("Coulomb")) m_input_HC_CSI_rate_av.push_back(rate); //Hz
331  if (HistoRateName.Contains("Def") && fileName.Contains("Coulomb")) m_input_HC_HE3_rate_av.push_back(rate); //Hz
332  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Coulomb")) m_input_HC_TPC_rate_av.push_back(rate); //Hz
333  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Coulomb")) m_input_HC_TPC_angular_rate_av.push_back(rate); //Hz
334  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Coulomb")) m_input_HC_TPC_angular_dose_av.push_back(rate); //Hz
335  }
336  if (fileName.Contains("LER")) {
337  if (HistoRateName.Contains("qcss") && fileName.Contains("Touschek")) m_input_LT_QCSS_rate.push_back(rate); //Hz
338  if (HistoRateName.Contains("claws") && fileName.Contains("Touschek")) m_input_LT_CLAWS_rate.push_back(rate); //Hz
339  if (HistoRateName.Contains("csi") && fileName.Contains("Touschek")) m_input_LT_CSI_rate.push_back(rate); //Hz
340  if (HistoRateName.Contains("Def") && fileName.Contains("Touschek")) m_input_LT_HE3_rate.push_back(rate); //Hz
341  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Touschek")) m_input_LT_TPC_rate.push_back(rate); //Hz
342  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Touschek")) m_input_LT_TPC_angular_rate.push_back(rate); //Hz
343  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Touschek")) m_input_LT_TPC_angular_dose.push_back(rate); //Hz
344  if (HistoRateName.Contains("qcss") && fileName.Contains("Brems")) m_input_LB_QCSS_rate_av.push_back(rate); //Hz
345  if (HistoRateName.Contains("claws") && fileName.Contains("Brems")) m_input_LB_CLAWS_rate_av.push_back(rate); //Hz
346  if (HistoRateName.Contains("csi") && fileName.Contains("Brems")) m_input_LB_CSI_rate_av.push_back(rate); //Hz
347  if (HistoRateName.Contains("Def") && fileName.Contains("Brems")) m_input_LB_HE3_rate_av.push_back(rate); //Hz
348  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Brems")) m_input_LB_TPC_rate_av.push_back(rate); //Hz
349  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Brems")) m_input_LB_TPC_angular_rate_av.push_back(rate); //Hz
350  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Brems")) m_input_LB_TPC_angular_dose_av.push_back(rate); //Hz
351  if (HistoRateName.Contains("qcss") && fileName.Contains("Coulomb")) m_input_LC_QCSS_rate_av.push_back(rate); //Hz
352  if (HistoRateName.Contains("claws") && fileName.Contains("Coulomb")) m_input_LC_CLAWS_rate_av.push_back(rate); //Hz
353  if (HistoRateName.Contains("csi") && fileName.Contains("Coulomb")) m_input_LC_CSI_rate_av.push_back(rate); //Hz
354  if (HistoRateName.Contains("Def") && fileName.Contains("Coulomb")) m_input_LC_HE3_rate_av.push_back(rate); //Hz
355  if (HistoRateName.Contains("tpc_rate") && fileName.Contains("Coulomb")) m_input_LC_TPC_rate_av.push_back(rate); //Hz
356  if (HistoRateName.Contains("tpc_angular_rate") && fileName.Contains("Coulomb")) m_input_LC_TPC_angular_rate_av.push_back(rate); //Hz
357  if (HistoRateName.Contains("tpc_angular_dose") && fileName.Contains("Coulomb")) m_input_LC_TPC_angular_dose_av.push_back(rate); //Hz
358  }
359  }
360  delete h1D;
361  }
362  for (const string& HistoDoseNameTmp : m_inputDoseHistoNames) {
363  const TString& HistoDoseName = (TString)HistoDoseNameTmp;
364 
365  int imax = 0;
366  if (HistoDoseName.Contains("csi")) imax = 18;
367  if (HistoDoseName.Contains("bgo")) imax = 8;
368  if (HistoDoseName.Contains("pin")) {
369  imax = 64;
370  volume = 0.265 * 0.265 * m_input_PIN_width; //cm^3
371  rho = 2.32; //g/cm^3
372  mass = rho * volume * 1e-3; //g to kg
373  }
374  if (HistoDoseName.Contains("dia")) {
375  imax = 4;
376  volume = 0.4 * 0.4 * 0.05; //cm^3
377  rho = 3.53; //g/cm^3
378  mass = rho * volume * 1e-3; //g to kg
379  }
380  if (HistoDoseName.Contains("dosi")) {
381  imax = 4;
382  volume = 0.265 * 0.265 * 0.01; //cm^3
383  rho = 2.32; //g/cm^3
384  mass = rho * volume * 1e-3; //g to kg
385  }
386  if (HistoDoseName.Contains("tpc")) {
387  imax = 2;
388  volume = 10.8537 * 2.0 * 1.68; //cm^3
389  rho = 0.00066908; //g/cm^3
390  mass = rho * volume * 1e-3; //g to kg
391  }
392  for (int i = 0; i < imax; i++) {
393 
394  TH1F* he;
395  if (HistoDoseName.Contains("csi")) {
396  he = (TH1F*)fh[iter]->Get(TString::Format("csi_dedep_%d_%d", i, m_DayBin));
397  } else {
398  he = (TH1F*)fh[iter]->Get(TString::Format("%s_%d", HistoDoseName.Data(), i));
399  }
400 
401  //double step = ((double)he->GetXaxis()->GetXmax() - (double)he->GetXaxis()->GetXmin()) / ((double)he->GetNbinsX());
402  double esum = 0;
403  for (int j = 0; j < he->GetNbinsX(); j++) {
404  double co = he->GetBinContent(j + 1);
405  double va = he->GetXaxis()->GetBinCenter(j + 1);
406  double esumbin = va * co;
407 
408  if (fileName.Contains("Coulomb")) {
409  if (fileName.Contains("HER")) esumbin *= m_input_Z_scaling[0];
410  if (fileName.Contains("LER")) esumbin *= m_input_Z_scaling[1];
411  }
412  if (fileName.Contains("Brems")) {
413  if (fileName.Contains("HER")) esumbin *= m_input_Z_scaling[2];
414  if (fileName.Contains("LER")) esumbin *= m_input_Z_scaling[3];
415  }
416 
417  esum += esumbin;// / step;
418  if (HistoDoseName.Contains("csi_energy")) {
419  if (fileName.Contains("HER")) {
420  if (fileName.Contains("Touschek")) m_input_HT_CSI_dose_binE.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
421  if (fileName.Contains("Coulomb")) m_input_HC_CSI_dose_binE_av.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
422  if (fileName.Contains("Brems")) m_input_HB_CSI_dose_binE_av.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
423  }
424  if (fileName.Contains("LER")) {
425  if (fileName.Contains("Touschek")) m_input_LT_CSI_dose_binE.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
426  if (fileName.Contains("Coulomb")) m_input_LC_CSI_dose_binE_av.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
427  if (fileName.Contains("Brems")) m_input_LB_CSI_dose_binE_av.push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
428  }
429  }
430  }
431  if (fileName.Contains("HER")) {
432  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
433  && fileName.Contains("Touschek")) m_input_HT_CSI_dose.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
434  if (HistoDoseName.Contains("bgo")
435  && fileName.Contains("Touschek")) m_input_HT_BGO_dose.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
436  if (HistoDoseName.Contains("pin")
437  && fileName.Contains("Touschek")) m_input_HT_PIN_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
438  if (HistoDoseName.Contains("dosi")
439  && fileName.Contains("Touschek")) m_input_HT_DOSI.push_back(esum / m_input_Time_eqv / mass / RadConv);
440  if (HistoDoseName.Contains("tpc_dose")
441  && fileName.Contains("Touschek")) m_input_HT_TPC_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
442  /*if (HistoDoseName.Contains("tpc_angular_dose")
443  && fileName.Contains("Touschek")) m_input_HT_TPC_angular_dose.push_back(esum / m_input_Time_eqv / mass / RadConv *
444  1e-3); //keV to MeV*/
445  if (HistoDoseName.Contains("dia")
446  && fileName.Contains("Touschek")) m_input_HT_DIA_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
447  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
448  && fileName.Contains("Brems")) m_input_HB_CSI_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
449  if (HistoDoseName.Contains("bgo")
450  && fileName.Contains("Brems")) m_input_HB_BGO_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
451  if (HistoDoseName.Contains("pin")
452  && fileName.Contains("Brems")) m_input_HB_PIN_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
453  if (HistoDoseName.Contains("dosi")
454  && fileName.Contains("Brems")) m_input_HB_DOSI_av.push_back(esum / m_input_Time_eqv / mass / RadConv);
455  if (HistoDoseName.Contains("tpc_dose")
456  && fileName.Contains("Brems")) m_input_HB_TPC_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
457  /*if (HistoDoseName.Contains("tpc_angular_dose")
458  && fileName.Contains("Brems")) m_input_HB_TPC_angular_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv *
459  1e-3); //keV to MeV*/
460  if (HistoDoseName.Contains("dia")
461  && fileName.Contains("Brems")) m_input_HB_DIA_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
462  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
463  && fileName.Contains("Coulomb")) m_input_HC_CSI_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
464  if (HistoDoseName.Contains("bgo")
465  && fileName.Contains("Coulomb")) m_input_HC_BGO_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
466  if (HistoDoseName.Contains("pin")
467  && fileName.Contains("Coulomb")) m_input_HC_PIN_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
468  if (HistoDoseName.Contains("dosi")
469  && fileName.Contains("Coulomb")) m_input_HC_DOSI_av.push_back(esum / m_input_Time_eqv / mass / RadConv);
470  if (HistoDoseName.Contains("tpc_dose")
471  && fileName.Contains("Coulomb")) m_input_HC_TPC_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
472  /*if (HistoDoseName.Contains("tpc_angular_dose")
473  && fileName.Contains("Coulomb")) m_input_HC_TPC_angular_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv *
474  1e-3); //keV to MeV*/
475  if (HistoDoseName.Contains("dia")
476  && fileName.Contains("Coulomb")) m_input_HC_DIA_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
477  }
478  if (fileName.Contains("LER")) {
479  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
480  && fileName.Contains("Touschek")) m_input_LT_CSI_dose.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
481  if (HistoDoseName.Contains("bgo")
482  && fileName.Contains("Touschek")) m_input_LT_BGO_dose.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
483  if (HistoDoseName.Contains("pin")
484  && fileName.Contains("Touschek")) m_input_LT_PIN_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
485  if (HistoDoseName.Contains("dosi")
486  && fileName.Contains("Touschek")) m_input_LT_DOSI.push_back(esum / m_input_Time_eqv / mass / RadConv);
487  if (HistoDoseName.Contains("tpc_dose")
488  && fileName.Contains("Touschek")) m_input_LT_TPC_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
489  /*if (HistoDoseName.Contains("tpc_angular_dose")
490  && fileName.Contains("Touschek")) m_input_LT_TPC_angular_dose.push_back(esum / m_input_Time_eqv / mass / RadConv *
491  1e-3); //keV to MeV*/
492  if (HistoDoseName.Contains("dia")
493  && fileName.Contains("Touschek")) m_input_LT_DIA_dose.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
494  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
495  && fileName.Contains("Brems")) m_input_LB_CSI_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
496  if (HistoDoseName.Contains("bgo")
497  && fileName.Contains("Brems")) m_input_LB_BGO_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
498  if (HistoDoseName.Contains("pin")
499  && fileName.Contains("Brems")) m_input_LB_PIN_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
500  if (HistoDoseName.Contains("dosi")
501  && fileName.Contains("Brems")) m_input_LB_DOSI_av.push_back(esum / m_input_Time_eqv / mass / RadConv);
502  if (HistoDoseName.Contains("tpc_dose")
503  && fileName.Contains("Brems")) m_input_LB_TPC_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
504  /*if (HistoDoseName.Contains("tpc_angular_dose")
505  && fileName.Contains("Brems")) m_input_LB_TPC_angular_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv *
506  1e-3); //keV to MeV*/
507  if (HistoDoseName.Contains("dia")
508  && fileName.Contains("Brems")) m_input_LB_DIA_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
509  if (HistoDoseName.Contains("csi") && HistoDoseName.Contains("edep")
510  && fileName.Contains("Coulomb")) m_input_LC_CSI_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
511  if (HistoDoseName.Contains("bgo")
512  && fileName.Contains("Coulomb")) m_input_LC_BGO_dose_av.push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
513  if (HistoDoseName.Contains("pin")
514  && fileName.Contains("Coulomb")) m_input_LC_PIN_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
515  if (HistoDoseName.Contains("dosi")
516  && fileName.Contains("Coulomb")) m_input_LC_DOSI_av.push_back(esum / m_input_Time_eqv / mass / RadConv);
517  if (HistoDoseName.Contains("tpc_dose")
518  && fileName.Contains("Coulomb")) m_input_LC_TPC_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
519  /*if (HistoDoseName.Contains("tpc_angular_dose")
520  && fileName.Contains("Coulomb")) m_input_LC_TPC_angular_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv *
521  1e-3); //keV to MeV*/
522  if (HistoDoseName.Contains("dia")
523  && fileName.Contains("Coulomb")) m_input_LC_DIA_dose_av.push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
524  }
525  delete he;
526  }
527  }
528  }
529  if (fileName.Contains("Coulomb") || fileName.Contains("Brems")) {
530  for (const string& HistoRateNameTmp : m_inputRateHistoNamesVrs) {
531  const TString& HistoRateName = (TString)HistoRateNameTmp;
532 
533  TH2F* h2D;
534  if (HistoRateName.Contains("csi")) h2D = (TH2F*)fh[iter]->Get(TString::Format("csi_rs_drate_%d", m_DayBin));
535  else h2D = (TH2F*)fh[iter]->Get(HistoRateName);
536 
537  for (int k = 0; k < h2D->GetNbinsY(); k++) {
538  for (int i = 0; i < h2D->GetNbinsX(); i++) {
539  double counts = h2D->GetBinContent(i + 1, k + 1);
540  double rate = counts / m_input_Time_eqv;
541 
542  if (fileName.Contains("Coulomb")) {
543  if (fileName.Contains("HER")) rate *= m_input_Z_scaling[0];
544  if (fileName.Contains("LER")) rate *= m_input_Z_scaling[1];
545  }
546  if (fileName.Contains("Brems")) {
547  if (fileName.Contains("HER")) rate *= m_input_Z_scaling[2];
548  if (fileName.Contains("LER")) rate *= m_input_Z_scaling[3];
549  }
550 
551  if (HistoRateName.Contains("Def")) rate *= m_input_HE3_EfCor[i];
552 
553  if (fileName.Contains("Coulomb_HER")) {
554  if (HistoRateName.Contains("qcss")) m_input_HC_QCSS_rate[k].push_back(rate); //Hz
555  if (HistoRateName.Contains("claws")) m_input_HC_CLAWS_rate[k].push_back(rate); //Hz
556  if (HistoRateName.Contains("csi")) m_input_HC_CSI_rate[k].push_back(rate); //Hz
557  if (HistoRateName.Contains("Def")) m_input_HC_HE3_rate[k].push_back(rate); //Hz
558  if (HistoRateName.Contains("tpc_rate")) m_input_HC_TPC_rate[k].push_back(rate); //Hz
559  if (HistoRateName.Contains("tpc_angular_rate")) m_input_HC_TPC_angular_rate[k].push_back(rate); //Hz
560  if (HistoRateName.Contains("tpc_angular_dose")) m_input_HC_TPC_angular_dose[k].push_back(rate); //Hz
561  }
562  if (fileName.Contains("Coulomb_LER")) {
563  if (HistoRateName.Contains("qcss")) m_input_LC_QCSS_rate[k].push_back(rate); //Hz
564  if (HistoRateName.Contains("claws")) m_input_LC_CLAWS_rate[k].push_back(rate); //Hz
565  if (HistoRateName.Contains("csi")) m_input_LC_CSI_rate[k].push_back(rate); //Hz
566  if (HistoRateName.Contains("Def")) m_input_LC_HE3_rate[k].push_back(rate); //Hz
567  if (HistoRateName.Contains("tpc_rate")) m_input_LC_TPC_rate[k].push_back(rate); //Hz
568  if (HistoRateName.Contains("tpc_angular_rate")) m_input_LC_TPC_angular_rate[k].push_back(rate); //Hz
569  if (HistoRateName.Contains("tpc_angular_dose")) m_input_LC_TPC_angular_dose[k].push_back(rate); //Hz
570  }
571  if (fileName.Contains("Brems_HER")) {
572  if (HistoRateName.Contains("qcss")) m_input_HB_QCSS_rate[k].push_back(rate); //Hz
573  if (HistoRateName.Contains("claws")) m_input_HB_CLAWS_rate[k].push_back(rate); //Hz
574  if (HistoRateName.Contains("csi")) m_input_HB_CSI_rate[k].push_back(rate); //Hz
575  if (HistoRateName.Contains("Def")) m_input_HB_HE3_rate[k].push_back(rate); //Hz
576  if (HistoRateName.Contains("tpc_rate")) m_input_HB_TPC_rate[k].push_back(rate); //Hz
577  if (HistoRateName.Contains("tpc_angular_rate")) m_input_HB_TPC_angular_rate[k].push_back(rate); //Hz
578  if (HistoRateName.Contains("tpc_angular_dose")) m_input_HB_TPC_angular_dose[k].push_back(rate); //Hz
579  }
580  if (fileName.Contains("Brems_LER")) {
581  if (HistoRateName.Contains("qcss")) m_input_LB_QCSS_rate[k].push_back(rate); //Hz
582  if (HistoRateName.Contains("claws")) m_input_LB_CLAWS_rate[k].push_back(rate); //Hz
583  if (HistoRateName.Contains("csi")) m_input_LB_CSI_rate[k].push_back(rate); //Hz
584  if (HistoRateName.Contains("Def")) m_input_LB_HE3_rate[k].push_back(rate); //Hz
585  if (HistoRateName.Contains("tpc_rate")) m_input_LB_TPC_rate[k].push_back(rate); //Hz
586  if (HistoRateName.Contains("tpc_angular_rate")) m_input_LB_TPC_angular_rate[k].push_back(rate); //Hz
587  if (HistoRateName.Contains("tpc_angular_dose")) m_input_LB_TPC_angular_dose[k].push_back(rate); //Hz
588  }
589  }
590  }
591  delete h2D;
592  }
593 
594  for (const string& HistoDoseNameTmp : m_inputDoseHistoNamesVrs) {
595  const TString& HistoDoseName = (TString)HistoDoseNameTmp;
596  int imax = 0;
597  if (HistoDoseName.Contains("csi")) imax = 18;
598  if (HistoDoseName.Contains("bgo")) imax = 8;
599  if (HistoDoseName.Contains("pin")) {
600  imax = 64;
601  volume = 0.265 * 0.265 * m_input_PIN_width; //cm^3
602  rho = 2.32; //g/cm^3
603  mass = rho * volume * 1e-3; //g to kg
604  }
605  if (HistoDoseName.Contains("dia")) {
606  imax = 4;
607  volume = 0.4 * 0.4 * 0.05; //cm^3
608  rho = 3.53; //g/cm^3
609  mass = rho * volume * 1e-3; //g to kg
610  }
611  if (HistoDoseName.Contains("dosi")) {
612  imax = 4;
613  volume = 0.265 * 0.265 * 0.01; //cm^3
614  rho = 2.32; //g/cm^3
615  mass = rho * volume * 1e-3; //g to kg
616  }
617  if (HistoDoseName.Contains("tpc")) {
618  imax = 2;
619  volume = 10.8537 * 2.0 * 1.68; //cm^3
620  rho = 0.00066908; //g/cm^3
621  mass = rho * volume * 1e-3; //g to kg
622  }
623  for (int i = 0; i < imax; i++) {
624 
625  TH2F* he;
626  if (HistoDoseName.Contains("csi")) {
627  he = (TH2F*)fh[iter]->Get(TString::Format("csi_rs_dedep_%d_%d", i, m_DayBin));
628  } else {
629  he = (TH2F*)fh[iter]->Get(TString::Format("%s_%d", HistoDoseName.Data(), i));
630  }
631 
632  //double step = ((double)he->GetXaxis()->GetXmax() - (double)he->GetXaxis()->GetXmin()) / ((double)he->GetNbinsX());
633  for (int k = 0; k < he->GetNbinsY(); k++) {
634  double esum = 0;
635  for (int j = 0; j < he->GetNbinsX(); j++) {
636  double co = he->GetBinContent(j + 1, k + 1);
637  double va = he->GetXaxis()->GetBinCenter(j + 1);
638  double esumbin = va * co;
639 
640  if (fileName.Contains("Coulomb")) {
641  if (fileName.Contains("HER")) esumbin *= m_input_Z_scaling[0];
642  if (fileName.Contains("LER")) esumbin *= m_input_Z_scaling[1];
643  }
644  if (fileName.Contains("Brems")) {
645  if (fileName.Contains("HER")) esumbin *= m_input_Z_scaling[2];
646  if (fileName.Contains("LER")) esumbin *= m_input_Z_scaling[3];
647  }
648 
649  esum += esumbin;// / step;
650  if (HistoDoseName.Contains("csi_energy")) {
651  if (fileName.Contains("HER")) {
652  if (fileName.Contains("Coulomb")) m_input_HC_CSI_dose_binE[k].push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
653  if (fileName.Contains("Brems")) m_input_HB_CSI_dose_binE[k].push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
654  }
655  if (fileName.Contains("LER")) {
656  if (fileName.Contains("Coulomb")) m_input_LC_CSI_dose_binE[k].push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
657  if (fileName.Contains("Brems")) m_input_LB_CSI_dose_binE[k].push_back(esumbin / m_input_Time_eqv * 1e-3); //MeV to GeV
658  }
659  }
660  }
661 
662  if (fileName.Contains("Coulomb_HER")) {
663  if (HistoDoseName.Contains("csi")
664  && HistoDoseName.Contains("edep")) m_input_HC_CSI_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
665  if (HistoDoseName.Contains("bgo")) m_input_HC_BGO_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
666  if (HistoDoseName.Contains("pin")) m_input_HC_PIN_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
667  if (HistoDoseName.Contains("dosi")) m_input_HC_DOSI[k].push_back(esum / m_input_Time_eqv / mass / RadConv);
668  if (HistoDoseName.Contains("tpc_dose")) m_input_HC_TPC_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
669  1e-3); //keV to MeV
670  /*if (HistoDoseName.Contains("tpc_angular_dose")) m_input_HC_TPC_angular_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
671  1e-3); //keV to MeV*/
672  if (HistoDoseName.Contains("dia")) m_input_HC_DIA_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
673  }
674  if (fileName.Contains("Coulomb_LER")) {
675  if (HistoDoseName.Contains("csi")
676  && HistoDoseName.Contains("edep")) m_input_LC_CSI_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
677  if (HistoDoseName.Contains("bgo")) m_input_LC_BGO_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
678  if (HistoDoseName.Contains("pin")) m_input_LC_PIN_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
679  if (HistoDoseName.Contains("dosi")) m_input_LC_DOSI[k].push_back(esum / m_input_Time_eqv / mass / RadConv);
680  if (HistoDoseName.Contains("tpc_dose")) m_input_LC_TPC_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
681  1e-3); //keV to MeV
682  /*if (HistoDoseName.Contains("tpc_angular_dose")) m_input_LC_TPC_angular_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
683  1e-3); //keV to MeV*/
684  if (HistoDoseName.Contains("dia")) m_input_LC_DIA_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
685  }
686  if (fileName.Contains("Brems_HER")) {
687  if (HistoDoseName.Contains("csi")
688  && HistoDoseName.Contains("edep")) m_input_HB_CSI_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
689  if (HistoDoseName.Contains("bgo")) m_input_HB_BGO_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
690  if (HistoDoseName.Contains("pin")) m_input_HB_PIN_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
691  if (HistoDoseName.Contains("dosi")) m_input_HB_DOSI[k].push_back(esum / m_input_Time_eqv / mass / RadConv);
692  if (HistoDoseName.Contains("tpc_dose")) m_input_HB_TPC_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
693  1e-3); //keV to MeV
694  /*if (HistoDoseName.Contains("tpc_angular_dose")) m_input_HB_TPC_angular_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
695  1e-3); //keV to MeV*/
696  if (HistoDoseName.Contains("dia")) m_input_HB_DIA_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
697  }
698  if (fileName.Contains("Brems_LER")) {
699  if (HistoDoseName.Contains("csi")
700  && HistoDoseName.Contains("edep")) m_input_LB_CSI_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
701  if (HistoDoseName.Contains("bgo")) m_input_LB_BGO_dose[k].push_back(esum / m_input_Time_eqv * 1e-3); //MeV to GeV
702  if (HistoDoseName.Contains("pin")) m_input_LB_PIN_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
703  if (HistoDoseName.Contains("dosi")) m_input_LB_DOSI[k].push_back(esum / m_input_Time_eqv / mass / RadConv);
704  if (HistoDoseName.Contains("tpc_dose")) m_input_LB_TPC_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
705  1e-3); //keV to MeV
706  /*if (HistoDoseName.Contains("tpc_angular_dose")) m_input_LB_TPC_angular_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv *
707  1e-3); //keV to MeV*/
708  if (HistoDoseName.Contains("dia")) m_input_LB_DIA_dose[k].push_back(esum / m_input_Time_eqv / mass / RadConv * 1e-3); //keV to MeV
709  }
710  }
711  delete he;
712  }
713  }
714 
715  }
716  iter++;
717  }
718  dirh->cd();
719  dir->cd();
720  m_numEntries = m_tree->GetEntries();
721  cout << "m_numEntries " << m_numEntries << endl;
722  m_entryCounter = 0;
723  m_exp = 0;
724  // data store objects registration
725 
726  StoreObjPtr<EventMetaData> evtMetaData;
727  evtMetaData.registerInDataStore();
728 
729  m_file = new TFile(m_outputFileName.c_str(), "RECREATE");
730  m_treeTruth = new TTree("truth", "Truth table (simulation)");
731 
732  m_treeBEAST = new TTree("tout", "BEAST data tree (simulation)");
733 
734  m_treeBEAST->Branch("ts", &(m_beast.ts));
735  m_treeBEAST->Branch("event", &(m_beast.event));
736  m_treeBEAST->Branch("run", &(m_beast.run));
737  m_treeBEAST->Branch("subrun", &(m_beast.subrun));
738 
739  m_treeBEAST->Branch("SKB_HER_injectionFlag", &(m_beast.SKB_HER_injectionFlag));
740  m_treeBEAST->Branch("SKB_LER_injectionFlag", &(m_beast.SKB_LER_injectionFlag));
741  m_treeBEAST->Branch("SKB_HER_injectionFlag_safe", &(m_beast.SKB_HER_injectionFlag_safe));
742  m_treeBEAST->Branch("SKB_LER_injectionFlag_safe", &(m_beast.SKB_LER_injectionFlag_safe));
743  m_treeBEAST->Branch("SKB_HER_abortFlag", &(m_beast.SKB_HER_abortFlag));
744  m_treeBEAST->Branch("SKB_LER_abortFlag", &(m_beast.SKB_LER_abortFlag));
745  m_treeBEAST->Branch("SKB_HER_abortFlag_safe", &(m_beast.SKB_HER_abortFlag_safe));
746  m_treeBEAST->Branch("SKB_LER_abortFlag_safe", &(m_beast.SKB_LER_abortFlag_safe));
747  m_treeBEAST->Branch("SKB_Status", &(m_beast.SKB_Status));
748  m_treeBEAST->Branch("SKB_HER_injectionRate", &(m_beast.SKB_HER_injectionRate));
749  m_treeBEAST->Branch("SKB_LER_injectionRate", &(m_beast.SKB_LER_injectionRate));
750  m_treeBEAST->Branch("SKB_HER_lifetime", &(m_beast.SKB_HER_lifetime));
751  m_treeBEAST->Branch("SKB_LER_lifetime", &(m_beast.SKB_LER_lifetime));
752  m_treeBEAST->Branch("SKB_LER_current", &(m_beast.SKB_LER_current));
753  m_treeBEAST->Branch("SKB_HER_current", &(m_beast.SKB_HER_current));
754  m_treeBEAST->Branch("SKB_LER_injectionEfficiency", &(m_beast.SKB_LER_injectionEfficiency));
755  m_treeBEAST->Branch("SKB_HER_injectionEfficiency", &(m_beast.SKB_HER_injectionEfficiency));
756  m_treeBEAST->Branch("SKB_beamLoss_ionChambers_mean", &(m_beast.SKB_beamLoss_ionChambers_mean));
757  m_treeBEAST->Branch("SKB_beamLoss_PINdiodes_mean", &(m_beast.SKB_beamLoss_PINdiodes_mean));
758  m_treeBEAST->Branch("SKB_beamLoss_nearCollimators", &(m_beast.SKB_beamLoss_nearCollimators));
759  m_treeBEAST->Branch("SKB_beamLoss_aroundMasks", &(m_beast.SKB_beamLoss_aroundMasks));
760  m_treeBEAST->Branch("SKB_LER_injectionCharge", &(m_beast.SKB_LER_injectionCharge));
761  m_treeBEAST->Branch("SKB_HER_injectionCharge", &(m_beast.SKB_HER_injectionCharge));
762  m_treeBEAST->Branch("SKB_LER_injectionRepetitionRate", &(m_beast.SKB_LER_injectionRepetitionRate));
763  m_treeBEAST->Branch("SKB_HER_injectionRepetitionRate", &(m_beast.SKB_HER_injectionRepetitionRate));
764  m_treeBEAST->Branch("SKB_LER_injectionNumberOfBunches", &(m_beast.SKB_LER_injectionNumberOfBunches));
765  m_treeBEAST->Branch("SKB_HER_injectionNumberOfBunches", &(m_beast.SKB_HER_injectionNumberOfBunches));
766  m_treeBEAST->Branch("SKB_LER_pressures", &(m_beast.SKB_LER_pressures));
767  m_treeBEAST->Branch("SKB_HER_pressures", &(m_beast.SKB_HER_pressures));
768  m_treeBEAST->Branch("SKB_LER_pressure_average", &(m_beast.SKB_LER_pressure_average));
769  m_treeBEAST->Branch("SKB_HER_pressure_average", &(m_beast.SKB_HER_pressure_average));
770  m_treeBEAST->Branch("SKB_LER_pressures_corrected", &(m_beast.SKB_LER_pressures_corrected));
771  m_treeBEAST->Branch("SKB_HER_pressures_corrected", &(m_beast.SKB_HER_pressures_corrected));
772  m_treeBEAST->Branch("SKB_LER_pressure_average_corrected", &(m_beast.SKB_LER_pressure_average_corrected));
773  m_treeBEAST->Branch("SKB_HER_pressure_average_corrected", &(m_beast.SKB_HER_pressure_average_corrected));
774  m_treeBEAST->Branch("SKB_HER_collimatorPositions_mm", &(m_beast.SKB_HER_collimatorPositions_mm));
775  m_treeBEAST->Branch("SKB_HER_collimatorPositions_DMM", &(m_beast.SKB_HER_collimatorPositions_DMM));
776  m_treeBEAST->Branch("SKB_HER_collimatorPositions_inX", &(m_beast.SKB_HER_collimatorPositions_inX));
777  m_treeBEAST->Branch("SKB_HER_collimatorPositions_inY", &(m_beast.SKB_HER_collimatorPositions_inY));
778  m_treeBEAST->Branch("SKB_HER_collimatorPositions_fromBeam", &(m_beast.SKB_HER_collimatorPositions_fromBeam));
779  m_treeBEAST->Branch("SKB_LER_collimatorPositions_mm", &(m_beast.SKB_LER_collimatorPositions_mm));
780  m_treeBEAST->Branch("SKB_LER_collimatorPositions_X", &(m_beast.SKB_LER_collimatorPositions_X));
781  m_treeBEAST->Branch("SKB_LER_collimatorPositions_Y", &(m_beast.SKB_LER_collimatorPositions_Y));
782  m_treeBEAST->Branch("SKB_LER_collimatorPositions_fromBeam", &(m_beast.SKB_LER_collimatorPositions_fromBeam));
783  m_treeBEAST->Branch("SKB_HER_beamSize_xray_X", &(m_beast.SKB_HER_beamSize_xray_X));
784  m_treeBEAST->Branch("SKB_HER_beamSize_xray_Y", &(m_beast.SKB_HER_beamSize_xray_Y));
785  m_treeBEAST->Branch("SKB_HER_correctedBeamSize_xray_Y", &(m_beast.SKB_HER_correctedBeamSize_xray_Y));
786  m_treeBEAST->Branch("SKB_LER_beamSize_xray_X", &(m_beast.SKB_LER_beamSize_xray_X));
787  m_treeBEAST->Branch("SKB_LER_beamSize_xray_Y", &(m_beast.SKB_LER_beamSize_xray_Y));
788  m_treeBEAST->Branch("SKB_LER_correctedBeamSize_xray_Y", &(m_beast.SKB_LER_correctedBeamSize_xray_Y));
789  m_treeBEAST->Branch("SKB_LER_beamSize_SR_X", &(m_beast.SKB_LER_beamSize_SR_X));
790  m_treeBEAST->Branch("SKB_LER_beamSize_SR_Y", &(m_beast.SKB_LER_beamSize_SR_Y));
791  m_treeBEAST->Branch("SKB_HER_beamSize_SR_X", &(m_beast.SKB_HER_beamSize_SR_X));
792  m_treeBEAST->Branch("SKB_HER_beamSize_SR_Y", &(m_beast.SKB_HER_beamSize_SR_Y));
793  m_treeBEAST->Branch("SKB_HER_integratedCurrent", &(m_beast.SKB_HER_integratedCurrent));
794  m_treeBEAST->Branch("SKB_LER_integratedCurrent", &(m_beast.SKB_LER_integratedCurrent));
795  m_treeBEAST->Branch("SKB_LER_partialPressures_D06", &(m_beast.SKB_LER_partialPressures_D06));
796  m_treeBEAST->Branch("SKB_LER_partialPressures_D02", &(m_beast.SKB_LER_partialPressures_D02));
797  m_treeBEAST->Branch("SKB_LER_pressures_local", &(m_beast.SKB_LER_pressures_local));
798  m_treeBEAST->Branch("SKB_HER_pressures_local", &(m_beast.SKB_HER_pressures_local));
799  m_treeBEAST->Branch("SKB_LER_pressures_local_corrected", &(m_beast.SKB_LER_pressures_local_corrected));
800  m_treeBEAST->Branch("SKB_HER_pressures_local_corrected", &(m_beast.SKB_HER_pressures_local_corrected));
801  m_treeBEAST->Branch("SKB_LER_Zeff_D02", &(m_beast.SKB_LER_Zeff_D02));
802  m_treeBEAST->Branch("SKB_LER_Zeff_D06", &(m_beast.SKB_LER_Zeff_D06));
803 
804  m_treeBEAST->Branch("CSI_data_sumE", &(m_beast.CSI_data_sumE));
805  m_treeBEAST->Branch("BGO_data_dose", &(m_beast.BGO_data_dose));
806  m_treeBEAST->Branch("PIN_data_dose", &(m_beast.PIN_data_dose));
807  m_treeBEAST->Branch("DIA_data_dose", &(m_beast.DIA_data_dose));
808 
809  m_treeBEAST->Branch("HE3_data_rate", &(m_beast.HE3_data_rate));
810  m_treeBEAST->Branch("CSI_data_rate", &(m_beast.CSI_data_rate));
811  m_treeBEAST->Branch("QCSS_data_rate", &(m_beast.QCSS_data_rate));
812  m_treeBEAST->Branch("CLAWS_data_rate", &(m_beast.CLAWS_data_rate));
813 
814  m_treeBEAST->Branch("CSI_data_Ebin", &(m_beast.CSI_data_Ebin));
815 
816  m_treeBEAST->Branch("PIN_dose", &(m_beast.PIN_dose));
817  m_treeBEAST->Branch("BGO_energy", &(m_beast.BGO_energy));
818  m_treeBEAST->Branch("HE3_rate", &(m_beast.HE3_rate));
819  //m_treeBEAST->Branch("TPC_rate", &(m_beast.TPC_rate), "TPC_rate[2][5]/F");
820  //m_treeBEAST->Branch("TPC_dose", &(m_beast.TPC_dose), "TPC_dose[2][5]/F");
821  //m_treeBEAST->Branch("TPC_angular_rate", &(m_beast.TPC_angular_rate), "TPC_angular_rate[2][9][18]/F");
822  //m_treeBEAST->Branch("TPC_angular_dose", &(m_beast.TPC_angular_dose), "TPC_angular_dose[2][9][18]/F");
823  m_treeBEAST->Branch("CSI_sumE", &(m_beast.CSI_sumE));
824  m_treeBEAST->Branch("CSI_Ebin", &(m_beast.CSI_Ebin));
825  m_treeBEAST->Branch("CSI_hitRate", &(m_beast.CSI_hitRate));
826  m_treeBEAST->Branch("DIA_dose", &(m_beast.DIA_dose));
827  m_treeBEAST->Branch("CLAWS_rate", &(m_beast.CLAWS_rate));
828  m_treeBEAST->Branch("QCSS_rate", &(m_beast.QCSS_rate));
829 
830  m_treeBEAST->Branch("DOSI_av", &(m_beast.DOSI_av));
831  m_treeBEAST->Branch("DOSI", &(m_beast.DOSI));
832  m_treeBEAST->Branch("PIN_dose_av", &(m_beast.PIN_dose_av));
833  m_treeBEAST->Branch("BGO_energy_av", &(m_beast.BGO_energy_av));
834  m_treeBEAST->Branch("HE3_rate_av", &(m_beast.HE3_rate_av));
835  m_treeBEAST->Branch("TPC_rate_av", &(m_beast.TPC_rate_av), "TPC_rate_av[2][5]/F");
836  //m_treeBEAST->Branch("TPC_dose_av", &(m_beast.TPC_dose_av), "TPC_dose_av[2][5]/F");
837  m_treeBEAST->Branch("TPC_dose_av", &(m_beast.TPC_dose_av));
838  m_treeBEAST->Branch("TPC_angular_rate_av", &(m_beast.TPC_angular_rate_av), "TPC_angular_rate_av[2][9][18]/F");
839  m_treeBEAST->Branch("TPC_angular_dose_av", &(m_beast.TPC_angular_dose_av), "TPC_angular_dose_av[2][9][18]/F");
840  m_treeBEAST->Branch("CSI_sumE_av", &(m_beast.CSI_sumE_av));
841  m_treeBEAST->Branch("CSI_Ebin_av", &(m_beast.CSI_Ebin_av));
842  m_treeBEAST->Branch("CSI_hitRate_av", &(m_beast.CSI_hitRate_av));
843  m_treeBEAST->Branch("DIA_dose_av", &(m_beast.DIA_dose_av));
844  m_treeBEAST->Branch("CLAWS_rate_av", &(m_beast.CLAWS_rate_av));
845  m_treeBEAST->Branch("QCSS_rate_av", &(m_beast.QCSS_rate_av));
846 
847  m_treeBEAST->Branch("SAD_HER_lifetime", &(m_beast.SAD_HER_lifetime));
848  m_treeBEAST->Branch("SAD_LER_lifetime", &(m_beast.SAD_LER_lifetime));
849  m_treeBEAST->Branch("SAD_HER_lifetime_av", &(m_beast.SAD_HER_lifetime_av));
850  m_treeBEAST->Branch("SAD_LER_lifetime_av", &(m_beast.SAD_LER_lifetime_av));
851  m_treeBEAST->Branch("SAD_HER_RLR", &(m_beast.SAD_HER_RLR));
852  m_treeBEAST->Branch("SAD_LER_RLR", &(m_beast.SAD_LER_RLR));
853  m_treeBEAST->Branch("SKB_HER_RLR", &(m_beast.SKB_HER_RLR));
854  m_treeBEAST->Branch("SKB_LER_RLR", &(m_beast.SKB_LER_RLR));
855  m_treeBEAST->Branch("SAD_HER_RLR_av", &(m_beast.SAD_HER_RLR_av));
856  m_treeBEAST->Branch("SAD_LER_RLR_av", &(m_beast.SAD_LER_RLR_av));
857 
858  m_treeBEAST->Branch("mc_reweight_LERT", &(m_beast.mc_reweight_LERT));
859  m_treeBEAST->Branch("mc_reweight_LERC", &(m_beast.mc_reweight_LERC));
860  m_treeBEAST->Branch("mc_reweight_LERB", &(m_beast.mc_reweight_LERB));
861  m_treeBEAST->Branch("mc_reweight_HERT", &(m_beast.mc_reweight_HERT));
862  m_treeBEAST->Branch("mc_reweight_HERC", &(m_beast.mc_reweight_HERC));
863  m_treeBEAST->Branch("mc_reweight_HERB", &(m_beast.mc_reweight_HERB));
864 
865 
866  /*
867  m_treeBEAST->Branch("TPC_neutrons_N", &(m_beast.TPC_neutrons_N));
868  m_treeBEAST->Branch("TPC_neutrons_tracks_E", &(m_beast.TPC_neutrons_tracks_E));
869  m_treeBEAST->Branch("TPC_neutrons_phi", &(m_beast.TPC_neutrons_phi));
870  m_treeBEAST->Branch("TPC_neutrons_theta", &(m_beast.TPC_neutrons_theta));
871  m_treeBEAST->Branch("TPC_neutrons_length", &(m_beast.TPC_neutrons_length));
872  m_treeBEAST->Branch("TPC_alphas_top_N", &(m_beast.TPC_alphas_top_N));
873  m_treeBEAST->Branch("TPC_alphas_bot_N", &(m_beast.TPC_alphas_bot_N));
874  m_treeBEAST->Branch("TPC_alphas_top_tracks_dEdx", &(m_beast.TPC_alphas_top_tracks_dEdx));
875  m_treeBEAST->Branch("TPC_alphas_bot_tracks_dEdx", &(m_beast.TPC_alphas_bot_tracks_dEdx));
876  m_treeBEAST->Branch("TPC_xrays_N", &(m_beast.TPC_xrays_N));
877  m_treeBEAST->Branch("TPC_xrays_sumE", &(m_beast.TPC_xrays_sumE));
878  */
879  m_treeTruth->Branch("SAD_I_HER", &(m_input_I_HER));
880  m_treeTruth->Branch("SAD_I_LER", &(m_input_I_LER));
881  m_treeTruth->Branch("SAD_P_HER", &(m_input_P_HER));
882  m_treeTruth->Branch("SAD_P_LER", &(m_input_P_LER));
883  m_treeTruth->Branch("SAD_sigma_x_HER", &(m_input_sigma_x_HER));
884  m_treeTruth->Branch("SAD_sigma_x_LER", &(m_input_sigma_x_LER));
885  m_treeTruth->Branch("SAD_sigma_y_HER", &(m_input_sigma_y_HER));
886  m_treeTruth->Branch("SAD_sigma_y_LER", &(m_input_sigma_y_LER));
887  m_treeTruth->Branch("SAD_bunchNb_HER", &(m_input_bunchNb_HER));
888  m_treeTruth->Branch("SAD_bunchNb_LER", &(m_input_bunchNb_LER));
889 
890  for (int i = 0; i < 12; i ++) {
891  m_treeTruth->Branch(TString::Format("MC_LC_DIA_dose_%d", i), &(m_input_LC_DIA_dose[i]));
892  m_treeTruth->Branch(TString::Format("MC_HC_DIA_dose_%d", i), &(m_input_HC_DIA_dose[i]));
893  m_treeTruth->Branch(TString::Format("MC_LB_DIA_dose_%d", i), &(m_input_LB_DIA_dose[i]));
894  m_treeTruth->Branch(TString::Format("MC_HB_DIA_dose_%d", i), &(m_input_HB_DIA_dose[i]));
895 
896  m_treeTruth->Branch(TString::Format("MC_LC_PIN_dose_%d", i), &(m_input_LC_PIN_dose[i]));
897  m_treeTruth->Branch(TString::Format("MC_HC_PIN_dose_%d", i), &(m_input_HC_PIN_dose[i]));
898  m_treeTruth->Branch(TString::Format("MC_LB_PIN_dose_%d", i), &(m_input_LB_PIN_dose[i]));
899  m_treeTruth->Branch(TString::Format("MC_HB_PIN_dose_%d", i), &(m_input_HB_PIN_dose[i]));
900 
901  m_treeTruth->Branch(TString::Format("MC_LC_BGO_dose_%d", i), &(m_input_LC_BGO_dose[i]));
902  m_treeTruth->Branch(TString::Format("MC_HC_BGO_dose_%d", i), &(m_input_HC_BGO_dose[i]));
903  m_treeTruth->Branch(TString::Format("MC_LB_BGO_dose_%d", i), &(m_input_LB_BGO_dose[i]));
904  m_treeTruth->Branch(TString::Format("MC_HB_BGO_dose_%d", i), &(m_input_HB_BGO_dose[i]));
905 
906  m_treeTruth->Branch(TString::Format("MC_LC_HE3_rate_%d", i), &(m_input_LC_HE3_rate[i]));
907  m_treeTruth->Branch(TString::Format("MC_HC_HE3_rate_%d", i), &(m_input_HC_HE3_rate[i]));
908  m_treeTruth->Branch(TString::Format("MC_LB_HE3_rate_%d", i), &(m_input_LB_HE3_rate[i]));
909  m_treeTruth->Branch(TString::Format("MC_HB_HE3_rate_%d", i), &(m_input_HB_HE3_rate[i]));
910 
911  //m_treeTruth->Branch(TString::Format("MC_LC_TPC_rate_%d", i), &(m_input_LC_TPC_rate[i]));
912  //m_treeTruth->Branch(TString::Format("MC_HC_TPC_rate_%d", i), &(m_input_HC_TPC_rate[i]));
913  //m_treeTruth->Branch(TString::Format("MC_LB_TPC_rate_%d", i), &(m_input_LB_TPC_rate[i]));
914  //m_treeTruth->Branch(TString::Format("MC_HB_TPC_rate_%d", i), &(m_input_HB_TPC_rate[i]));
915 
916  //m_treeTruth->Branch(TString::Format("MC_LC_TPC_dose_%d", i), &(m_input_LC_TPC_dose[i]));
917  //m_treeTruth->Branch(TString::Format("MC_HC_TPC_dose_%d", i), &(m_input_HC_TPC_dose[i]));
918  //m_treeTruth->Branch(TString::Format("MC_LB_TPC_dose_%d", i), &(m_input_LB_TPC_dose[i]));
919  //m_treeTruth->Branch(TString::Format("MC_HB_TPC_dose_%d", i), &(m_input_HB_TPC_dose[i]));
920 
921  //m_treeTruth->Branch(TString::Format("MC_LC_TPC_angular_rate_%d", i), &(m_input_LC_TPC_angular_rate[i]));
922  //m_treeTruth->Branch(TString::Format("MC_HC_TPC_angular_rate_%d", i), &(m_input_HC_TPC_angular_rate[i]));
923  //m_treeTruth->Branch(TString::Format("MC_LB_TPC_angular_rate_%d", i), &(m_input_LB_TPC_angular_rate[i]));
924  //m_treeTruth->Branch(TString::Format("MC_HB_TPC_angular_rate_%d", i), &(m_input_HB_TPC_angular_rate[i]));
925 
926  //m_treeTruth->Branch(TString::Format("MC_LC_TPC_angular_dose_%d", i), &(m_input_LC_TPC_angular_dose[i]));
927  //m_treeTruth->Branch(TString::Format("MC_HC_TPC_angular_dose_%d", i), &(m_input_HC_TPC_angular_dose[i]));
928  //m_treeTruth->Branch(TString::Format("MC_LB_TPC_angular_dose_%d", i), &(m_input_LB_TPC_angular_dose[i]));
929  //m_treeTruth->Branch(TString::Format("MC_HB_TPC_angular_dose_%d", i), &(m_input_HB_TPC_angular_dose[i]));
930 
931  m_treeTruth->Branch(TString::Format("MC_LC_CSI_rate_%d", i), &(m_input_LC_CSI_rate[i]));
932  m_treeTruth->Branch(TString::Format("MC_HC_CSI_rate_%d", i), &(m_input_HC_CSI_rate[i]));
933  m_treeTruth->Branch(TString::Format("MC_LB_CSI_rate_%d", i), &(m_input_LB_CSI_rate[i]));
934  m_treeTruth->Branch(TString::Format("MC_HB_CSI_rate_%d", i), &(m_input_HB_CSI_rate[i]));
935 
936  m_treeTruth->Branch(TString::Format("MC_LC_CSI_dose_%d", i), &(m_input_LC_CSI_dose[i]));
937  m_treeTruth->Branch(TString::Format("MC_HC_CSI_dose_%d", i), &(m_input_HC_CSI_dose[i]));
938  m_treeTruth->Branch(TString::Format("MC_LB_CSI_dose_%d", i), &(m_input_LB_CSI_dose[i]));
939  m_treeTruth->Branch(TString::Format("MC_HB_CSI_dose_%d", i), &(m_input_HB_CSI_dose[i]));
940 
941  m_treeTruth->Branch(TString::Format("MC_LC_CLAWS_rate_%d", i), &(m_input_LC_CLAWS_rate[i]));
942  m_treeTruth->Branch(TString::Format("MC_HC_CLAWS_rate_%d", i), &(m_input_HC_CLAWS_rate[i]));
943  m_treeTruth->Branch(TString::Format("MC_LB_CLAWS_rate_%d", i), &(m_input_LB_CLAWS_rate[i]));
944  m_treeTruth->Branch(TString::Format("MC_HB_CLAWS_rate_%d", i), &(m_input_HB_CLAWS_rate[i]));
945 
946  m_treeTruth->Branch(TString::Format("MC_LC_QCSS_rate_%d", i), &(m_input_LC_QCSS_rate[i]));
947  m_treeTruth->Branch(TString::Format("MC_HC_QCSS_rate_%d", i), &(m_input_HC_QCSS_rate[i]));
948  m_treeTruth->Branch(TString::Format("MC_LB_QCSS_rate_%d", i), &(m_input_LB_QCSS_rate[i]));
949  m_treeTruth->Branch(TString::Format("MC_HB_QCSS_rate_%d", i), &(m_input_HB_QCSS_rate[i]));
950  }
951 
952  m_treeTruth->Branch("MC_LT_DIA_dose", &(m_input_LT_DIA_dose));
953  m_treeTruth->Branch("MC_HT_DIA_dose", &(m_input_HT_DIA_dose));
954  m_treeTruth->Branch("MC_LC_DIA_dose_av", &(m_input_LC_DIA_dose_av));
955  m_treeTruth->Branch("MC_HC_DIA_dose_av", &(m_input_HC_DIA_dose_av));
956  m_treeTruth->Branch("MC_LB_DIA_dose_av", &(m_input_LB_DIA_dose_av));
957  m_treeTruth->Branch("MC_HB_DIA_dose_av", &(m_input_HB_DIA_dose_av));
958 
959  m_treeTruth->Branch("MC_LT_PIN_dose", &(m_input_LT_PIN_dose));
960  m_treeTruth->Branch("MC_HT_PIN_dose", &(m_input_HT_PIN_dose));
961  m_treeTruth->Branch("MC_LC_PIN_dose_av", &(m_input_LC_PIN_dose_av));
962  m_treeTruth->Branch("MC_HC_PIN_dose_av", &(m_input_HC_PIN_dose_av));
963  m_treeTruth->Branch("MC_LB_PIN_dose_av", &(m_input_LB_PIN_dose_av));
964  m_treeTruth->Branch("MC_HB_PIN_dose_av", &(m_input_HB_PIN_dose_av));
965 
966  m_treeTruth->Branch("MC_LT_BGO_dose", &(m_input_LT_BGO_dose));
967  m_treeTruth->Branch("MC_HT_BGO_dose", &(m_input_HT_BGO_dose));
968  m_treeTruth->Branch("MC_LC_BGO_dose_av", &(m_input_LC_BGO_dose_av));
969  m_treeTruth->Branch("MC_HC_BGO_dose_av", &(m_input_HC_BGO_dose_av));
970  m_treeTruth->Branch("MC_LB_BGO_dose_av", &(m_input_LB_BGO_dose_av));
971  m_treeTruth->Branch("MC_HB_BGO_dose_av", &(m_input_HB_BGO_dose_av));
972 
973  m_treeTruth->Branch("MC_LT_HE3_rate", &(m_input_LT_HE3_rate));
974  m_treeTruth->Branch("MC_HT_HE3_rate", &(m_input_HT_HE3_rate));
975  m_treeTruth->Branch("MC_LC_HE3_rate_av", &(m_input_LC_HE3_rate_av));
976  m_treeTruth->Branch("MC_HC_HE3_rate_av", &(m_input_HC_HE3_rate_av));
977  m_treeTruth->Branch("MC_LB_HE3_rate_av", &(m_input_LB_HE3_rate_av));
978  m_treeTruth->Branch("MC_HB_HE3_rate_av", &(m_input_HB_HE3_rate_av));
979 
980  //m_treeTruth->Branch("MC_LT_TPC_rate", &(m_input_LT_TPC_rate));
981  //m_treeTruth->Branch("MC_HT_TPC_rate", &(m_input_HT_TPC_rate));
982  m_treeTruth->Branch("MC_LC_TPC_rate_av", &(m_input_LC_TPC_rate_av));
983  m_treeTruth->Branch("MC_HC_TPC_rate_av", &(m_input_HC_TPC_rate_av));
984  m_treeTruth->Branch("MC_LB_TPC_rate_av", &(m_input_LB_TPC_rate_av));
985  m_treeTruth->Branch("MC_HB_TPC_rate_av", &(m_input_HB_TPC_rate_av));
986 
987  //m_treeTruth->Branch("MC_LT_TPC_dose", &(m_input_LT_TPC_dose));
988  //m_treeTruth->Branch("MC_HT_TPC_dose", &(m_input_HT_TPC_dose));
989  m_treeTruth->Branch("MC_LC_TPC_dose_av", &(m_input_LC_TPC_dose_av));
990  m_treeTruth->Branch("MC_HC_TPC_dose_av", &(m_input_HC_TPC_dose_av));
991  m_treeTruth->Branch("MC_LB_TPC_dose_av", &(m_input_LB_TPC_dose_av));
992  m_treeTruth->Branch("MC_HB_TPC_dose_av", &(m_input_HB_TPC_dose_av));
993 
994  //m_treeTruth->Branch("MC_LT_TPC_angular_rate", &(m_input_LT_TPC_angular_rate));
995  //m_treeTruth->Branch("MC_HT_TPC_angular_rate", &(m_input_HT_TPC_angular_rate));
996  m_treeTruth->Branch("MC_LC_TPC_angular_rate_av", &(m_input_LC_TPC_angular_rate_av));
997  m_treeTruth->Branch("MC_HC_TPC_angular_rate_av", &(m_input_HC_TPC_angular_rate_av));
998  m_treeTruth->Branch("MC_LB_TPC_angular_rate_av", &(m_input_LB_TPC_angular_rate_av));
999  m_treeTruth->Branch("MC_HB_TPC_angular_rate_av", &(m_input_HB_TPC_angular_rate_av));
1000 
1001  //m_treeTruth->Branch("MC_LT_TPC_angular_dose", &(m_input_LT_TPC_angular_dose));
1002  //m_treeTruth->Branch("MC_HT_TPC_angular_dose", &(m_input_HT_TPC_angular_dose));
1003  m_treeTruth->Branch("MC_LC_TPC_angular_dose_av", &(m_input_LC_TPC_angular_dose_av));
1004  m_treeTruth->Branch("MC_HC_TPC_angular_dose_av", &(m_input_HC_TPC_angular_dose_av));
1005  m_treeTruth->Branch("MC_LB_TPC_angular_dose_av", &(m_input_LB_TPC_angular_dose_av));
1006  m_treeTruth->Branch("MC_HB_TPC_angular_dose_av", &(m_input_HB_TPC_angular_dose_av));
1007 
1008  m_treeTruth->Branch("MC_LT_CSI_rate", &(m_input_LT_CSI_rate));
1009  m_treeTruth->Branch("MC_HT_CSI_rate", &(m_input_HT_CSI_rate));
1010  m_treeTruth->Branch("MC_LC_CSI_rate_av", &(m_input_LC_CSI_rate_av));
1011  m_treeTruth->Branch("MC_HC_CSI_rate_av", &(m_input_HC_CSI_rate_av));
1012  m_treeTruth->Branch("MC_LB_CSI_rate_av", &(m_input_LB_CSI_rate_av));
1013  m_treeTruth->Branch("MC_HB_CSI_rate_av", &(m_input_HB_CSI_rate_av));
1014 
1015  m_treeTruth->Branch("MC_LT_CSI_dose", &(m_input_LT_CSI_dose));
1016  m_treeTruth->Branch("MC_HT_CSI_dose", &(m_input_HT_CSI_dose));
1017  m_treeTruth->Branch("MC_LC_CSI_dose_av", &(m_input_LC_CSI_dose_av));
1018  m_treeTruth->Branch("MC_HC_CSI_dose_av", &(m_input_HC_CSI_dose_av));
1019  m_treeTruth->Branch("MC_LB_CSI_dose_av", &(m_input_LB_CSI_dose_av));
1020  m_treeTruth->Branch("MC_HB_CSI_dose_av", &(m_input_HB_CSI_dose_av));
1021 
1022  m_treeTruth->Branch("MC_LT_CLAWS_rate", &(m_input_LT_CLAWS_rate));
1023  m_treeTruth->Branch("MC_HT_CLAWS_rate", &(m_input_HT_CLAWS_rate));
1024  m_treeTruth->Branch("MC_LC_CLAWS_rate_av", &(m_input_LC_CLAWS_rate_av));
1025  m_treeTruth->Branch("MC_HC_CLAWS_rate_av", &(m_input_HC_CLAWS_rate_av));
1026  m_treeTruth->Branch("MC_LB_CLAWS_rate_av", &(m_input_LB_CLAWS_rate_av));
1027  m_treeTruth->Branch("MC_HB_CLAWS_rate_av", &(m_input_HB_CLAWS_rate_av));
1028 
1029  m_treeTruth->Branch("MC_LT_QCSS_rate", &(m_input_LT_QCSS_rate));
1030  m_treeTruth->Branch("MC_HT_QCSS_rate", &(m_input_HT_QCSS_rate));
1031  m_treeTruth->Branch("MC_LC_QCSS_rate_av", &(m_input_LC_QCSS_rate_av));
1032  m_treeTruth->Branch("MC_HC_QCSS_rate_av", &(m_input_HC_QCSS_rate_av));
1033  m_treeTruth->Branch("MC_LB_QCSS_rate_av", &(m_input_LB_QCSS_rate_av));
1034  m_treeTruth->Branch("MC_HB_QCSS_rate_av", &(m_input_HB_QCSS_rate_av));
1035 
1036  m_treeTruth->Fill();
1037  }
1038 
1039 
1040  void NtuplePhase1_v6Module::beginRun()
1041  {
1042  }
1043 
1044 
1045  void NtuplePhase1_v6Module::event()
1046  {
1047  m_beast.clear();
1048  // create data store objects
1049 
1050  StoreObjPtr<EventMetaData> evtMetaData;
1051  evtMetaData.create();
1052 
1053  if (m_eventCount == m_numEvents) {
1054  evtMetaData->setEndOfData(); // event processing
1055  return;
1056  }
1057 
1058  m_tree->GetEntry(m_eventCount);
1059 
1060  double Zeff_LER = 2.7;
1061  if (m_beast.SKB_LER_Zeff_D02 != 0 && m_beast.SKB_LER_Zeff_D02->size() > 0) Zeff_LER = m_beast.SKB_LER_Zeff_D02->at(0);
1062  //cout << "Zeff_DO2 " << Zeff_LER << endl;
1063  double Zeff_LC = 1;
1064  double Zeff_LB = 1;
1065  //if (Zeff_LER == 0) {
1066  //Zeff_LER = 2.7;
1067  //Zeff_LC = 1.0;
1068  //Zeff_LB = 1.0;
1069  //} else
1070  if (Zeff_LER > 0 && Zeff_LER < 40) {
1071  Zeff_LC = fctRate_LC->Eval(Zeff_LER) / fctRate_LC->Eval(7) / m_input_Z_scaling[1];
1072  Zeff_LB = fctRate_LB->Eval(Zeff_LER) / fctRate_LB->Eval(7) / m_input_Z_scaling[3];
1073  }
1074  //cout << "Zeff_LC " << Zeff_LC << " Zeff_LB " << Zeff_LB << " Zeff_LER " << Zeff_LER << endl;
1075  if (m_input_Z[1] == 2.7 && m_input_Z[3] == 2.7) {
1076  Zeff_LC = 1;
1077  Zeff_LB = 1;
1078  }
1079  double Zeff_HC = m_input_Z_scaling[0];
1080  double Zeff_HB = m_input_Z_scaling[2];
1081 
1082  double I_HER = 0;
1083  if (m_beast.SKB_HER_current != 0 && m_beast.SKB_HER_current->size() > 0) I_HER = m_beast.SKB_HER_current->at(0);
1084  if (m_input_I_HER[1] > 0) I_HER += gRandom->Gaus(0, m_input_I_HER[1]);
1085  double I_LER = 0;
1086  if (m_beast.SKB_LER_current != 0 && m_beast.SKB_LER_current->size() > 0) I_LER = m_beast.SKB_LER_current->at(0);
1087  if (m_input_I_LER[1] > 0) I_LER += gRandom->Gaus(0, m_input_I_LER[1]);
1088  /*double P_HER = 0;
1089  if (m_beast.SKB_HER_pressure_average != 0
1090  && m_beast.SKB_HER_pressure_average->size() > 0) P_HER = m_beast.SKB_HER_pressure_average->at(
1091  0) * 0.00750062 * 1e9 * m_input_GasCorrection[0];*/
1092  //if (m_input_P_HER[1] > 0) P_HER += gRandom->Gaus(0, m_input_P_HER[1]);
1093  //double P_LER = 0;
1094  /*if (m_beast.SKB_LER_pressure_average != 0
1095  && m_beast.SKB_LER_pressure_average->size() > 0) P_LER = m_beast.SKB_LER_pressure_average->at(
1096  0) * 0.00750062 * 1e9 * m_input_GasCorrection[1];*/
1097  //if (m_input_P_LER[1] > 0) P_LER += gRandom->Gaus(0, m_input_P_LER[1]);
1098  double P_corrected_HER = 0;
1099  if (m_beast.SKB_HER_pressure_average_corrected != 0
1100  && m_beast.SKB_HER_pressure_average_corrected->size() > 0) P_corrected_HER = m_beast.SKB_HER_pressure_average_corrected->at(
1101  0) * 0.00750062 * 1e9;
1102  if (m_input_P_HER[1] > 0) P_corrected_HER += gRandom->Gaus(0, m_input_P_HER[1]);
1103  double P_corrected_LER = 0;
1104  if (m_beast.SKB_LER_pressure_average_corrected != 0
1105  && m_beast.SKB_LER_pressure_average_corrected->size() > 0) P_corrected_LER = m_beast.SKB_LER_pressure_average_corrected->at(
1106  0) * 0.00750062 * 1e9;
1107  if (m_input_P_LER[1] > 0) P_corrected_LER += gRandom->Gaus(0, m_input_P_LER[1]);
1108 
1109  double sigma_y_HER = 0;
1110  double sigma_x_HER = 0;
1111  if (m_beast.SKB_HER_correctedBeamSize_xray_Y != 0
1112  && m_beast.SKB_HER_beamSize_xray_Y != 0
1113  && m_beast.SKB_HER_correctedBeamSize_xray_Y->size() > 0) {
1114  sigma_y_HER = m_beast.SKB_HER_correctedBeamSize_xray_Y->at(0);
1115 
1116  if (m_beast.SKB_HER_beamSize_SR_X != 0 && m_beast.SKB_HER_beamSize_SR_Y
1117  && m_beast.SKB_HER_beamSize_SR_X->size() > 0
1118  && m_beast.SKB_HER_beamSize_SR_Y->size() > 0)
1119  sigma_x_HER = m_beast.SKB_HER_correctedBeamSize_xray_Y->at(0) * m_beast.SKB_HER_beamSize_SR_X->at(
1120  0) / m_beast.SKB_HER_beamSize_SR_Y->at(0);
1121  }
1122  /*
1123  if (m_beast.SKB_HER_beamSize_xray_Y != 0
1124  && m_beast.SKB_HER_beamSize_xray_Y != 0
1125  && m_beast.SKB_HER_beamSize_xray_Y->size() > 0) sigma_y_HER = m_beast.SKB_HER_beamSize_xray_Y->at(0);
1126  */
1127  if (m_input_sigma_x_HER[1] > 0) sigma_x_HER += gRandom->Gaus(0, m_input_sigma_x_HER[1]);
1128  if (m_input_sigma_y_HER[1] > 0) sigma_y_HER += gRandom->Gaus(0, m_input_sigma_y_HER[1]);
1129  double sigma_y_LER = 0;
1130  double sigma_x_LER = 0;
1131  if (m_beast.SKB_LER_correctedBeamSize_xray_Y != 0
1132  && m_beast.SKB_LER_beamSize_xray_Y != 0
1133  && m_beast.SKB_LER_correctedBeamSize_xray_Y->size() > 0) {
1134  sigma_y_LER = m_beast.SKB_LER_correctedBeamSize_xray_Y->at(0);
1135 
1136  if (m_beast.SKB_LER_beamSize_SR_X != 0 && m_beast.SKB_LER_beamSize_SR_Y
1137  && m_beast.SKB_LER_beamSize_SR_X->size() > 0
1138  && m_beast.SKB_LER_beamSize_SR_Y->size() > 0)
1139  sigma_x_LER = m_beast.SKB_LER_correctedBeamSize_xray_Y->at(0) * m_beast.SKB_LER_beamSize_SR_X->at(
1140  0) / m_beast.SKB_LER_beamSize_SR_Y->at(0);
1141  }
1142  /*
1143  if (m_beast.SKB_LER_beamSize_xray_Y != 0
1144  && m_beast.SKB_LER_beamSize_xray_Y != 0
1145  && m_beast.SKB_LER_beamSize_xray_Y->size() > 0) sigma_y_LER = m_beast.SKB_LER_beamSize_xray_Y->at(0);
1146  */
1147  if (m_input_sigma_x_LER[1] > 0) sigma_x_LER += gRandom->Gaus(0, m_input_sigma_x_LER[1]);
1148  if (m_input_sigma_y_LER[1] > 0) sigma_y_LER += gRandom->Gaus(0, m_input_sigma_y_LER[1]);
1149  double bunch_nb_HER = 0;
1150  if (m_beast.SKB_HER_injectionNumberOfBunches != 0
1151  && m_beast.SKB_HER_injectionNumberOfBunches->size() > 0) bunch_nb_HER = m_beast.SKB_HER_injectionNumberOfBunches->at(0);
1152  if (bunch_nb_HER == 0) bunch_nb_HER = m_input_data_bunchNb_HER;
1153  if (m_input_bunchNb_HER[1] > 0) bunch_nb_HER += gRandom->Gaus(0, m_input_bunchNb_HER[1]);
1154  double bunch_nb_LER = 0;
1155  if (m_beast.SKB_LER_injectionNumberOfBunches != 0
1156  && m_beast.SKB_LER_injectionNumberOfBunches->size() > 0) bunch_nb_LER = m_beast.SKB_LER_injectionNumberOfBunches->at(0);
1157  if (bunch_nb_LER == 0) bunch_nb_LER = m_input_data_bunchNb_LER;
1158  if (m_input_bunchNb_LER[1] > 0) bunch_nb_LER += gRandom->Gaus(0, m_input_bunchNb_LER[1]);
1159 
1160  /*
1161  if (I_HER < 0 || I_HER > 2100.) I_HER = 0;
1162  if (I_LER < 0 || I_LER > 2100.) I_LER = 0;
1163  if (P_HER < 0 || P_HER > 760.) P_HER = 0;
1164  if (P_LER < 0 || P_LER > 760.) P_LER = 0;
1165  */
1166  if (I_HER < 0) I_HER = 0;
1167  if (I_LER < 0) I_LER = 0;
1168  //if (P_HER < 0) P_HER = 0;
1169  //if (P_LER < 0) P_LER = 0;
1170  if (P_corrected_HER < 0) P_corrected_HER = 0;
1171  if (P_corrected_LER < 0) P_corrected_LER = 0;
1172 
1173  if (m_input_data_SingleBeam == "LER") {
1174  I_HER = 0;
1175  //P_HER = 0;
1176  P_corrected_HER = 0;
1177  } else if (m_input_data_SingleBeam == "HER") {
1178  I_LER = 0;
1179  //P_LER = 0;
1180  P_corrected_LER = 0;
1181  }
1182 
1183  double Ib_HER = 0;
1184  if (bunch_nb_HER > 0) Ib_HER = I_HER / bunch_nb_HER * 1e-3; // mA -> A
1185  double Ib_LER = 0;
1186  if (bunch_nb_LER > 0) Ib_LER = I_LER / bunch_nb_LER * 1e-3; // mA -> A
1187 
1188  double Nb_HER = 0;
1189  if (Ib_HER > 0) Nb_HER = Ib_HER * 3000. / TMath::C() / (1.6e-19);
1190  double Nb_LER = 0;
1191  if (Ib_LER > 0) Nb_LER = Ib_LER * 3000. / TMath::C() / (1.6e-19);
1192 
1193  //double RLR_HER = 0;
1194  if (m_beast.SKB_HER_lifetime != 0 && m_beast.SKB_HER_lifetime->size() > 0 && Nb_HER > 0) {
1195  double RLR_HER = Nb_HER / (m_beast.SKB_HER_lifetime->at(0) * 60.) * 1e-9 * bunch_nb_HER;
1196  m_beast.SKB_HER_RLR.push_back(RLR_HER);
1197  }
1198  //double RLR_LER = 0;
1199  if (m_beast.SKB_LER_lifetime != 0 && m_beast.SKB_LER_lifetime->size() > 0 && Nb_LER > 0) {
1200  double RLR_LER = Nb_LER / (m_beast.SKB_LER_lifetime->at(0) * 60.) * 1e-9 * bunch_nb_LER;
1201  m_beast.SKB_LER_RLR.push_back(RLR_LER);
1202  }
1203 
1204  //Calculate Beam Gas scaling factor:
1205  // solution 0: Beam Gas \propo I x P => (IP)^data / (IP)^simu
1206  // solution 1: Beam Gas \propo I x P => (I_bP)^data / (I_bP)^simu where I_b = current / #bunch
1207  double ScaleFacBGav_HER = 0;
1208  double ScaleFacBGav_LER = 0;
1209  /*
1210  if (I_LER > 0 && P_LER > 0) {
1211  if (m_input_BGSol == 0
1212  && bunch_nb_LER > 0) ScaleFacBGav_LER = I_LER * P_LER / (m_input_I_LER[0] * m_input_P_LER[0]) / bunch_nb_LER *
1213  m_input_bunchNb_LER[0];
1214  if (m_input_BGSol == 1 && bunch_nb_LER > 0) ScaleFacBGav_LER = I_LER * P_LER / (m_input_I_LER[0] * m_input_P_LER[0]);
1215  }
1216  if (I_HER > 0 && P_HER > 0) {
1217  if (m_input_BGSol == 0
1218  && bunch_nb_HER > 0) ScaleFacBGav_HER = I_HER * P_HER / (m_input_I_HER[0] * m_input_P_HER[0]) / bunch_nb_HER *
1219  m_input_bunchNb_HER[0];
1220  if (m_input_BGSol == 1 && bunch_nb_HER > 0) ScaleFacBGav_HER = I_HER * P_HER / (m_input_I_HER[0] * m_input_P_HER[0]);
1221  }
1222  */
1223  if (I_LER > 0 && P_corrected_LER > 0) {
1224  if (m_input_BGSol == 0
1225  && bunch_nb_LER > 0) ScaleFacBGav_LER = I_LER * P_corrected_LER / (m_input_I_LER[0] * m_input_P_LER[0]) / bunch_nb_LER *
1226  m_input_bunchNb_LER[0];
1227  if (m_input_BGSol == 1 && bunch_nb_LER > 0) ScaleFacBGav_LER = I_LER * P_corrected_LER / (m_input_I_LER[0] * m_input_P_LER[0]);
1228  }
1229  if (I_HER > 0 && P_corrected_HER > 0) {
1230  if (m_input_BGSol == 0
1231  && bunch_nb_HER > 0) ScaleFacBGav_HER = I_HER * P_corrected_HER / (m_input_I_HER[0] * m_input_P_HER[0]) / bunch_nb_HER *
1232  m_input_bunchNb_HER[0];
1233  if (m_input_BGSol == 1 && bunch_nb_HER > 0) ScaleFacBGav_HER = I_HER * P_corrected_HER / (m_input_I_HER[0] * m_input_P_HER[0]);
1234  }
1235 
1236  //Calculate Beam Gas scaling factor: Beam Gas \propo I x P => (IP)^data / (IP)^simu
1237  double ScaleFacBG_HER[12];
1238  double ScaleFacBG_LER[12];
1239 
1240  if (m_beast.SKB_HER_pressures != 0 && m_input_GasCorrection[0] != 1) {
1241  for (int i = 0; i < (int)m_beast.SKB_HER_pressures->size(); i++) {
1242  ScaleFacBG_HER[i] = 0;
1243  double iP_HER = 0;
1244  iP_HER = m_beast.SKB_HER_pressures->at(i) * 0.00750062 * 1e9 * m_input_GasCorrection[0];
1245  if (m_input_P_HER[1] > 0) iP_HER += gRandom->Gaus(0, m_input_P_HER[1]);
1246  if (iP_HER < 0 || iP_HER > 260.) iP_HER = 0;
1247  if (I_HER > 0 && iP_HER > 0) {
1248  if (m_input_BGSol == 0
1249  && bunch_nb_HER > 0) ScaleFacBG_HER[i] = I_HER * iP_HER / (m_input_I_HER[0] * m_input_P_HER[0]) / bunch_nb_HER *
1250  m_input_bunchNb_HER[0];
1251  if (m_input_BGSol == 1 && bunch_nb_HER > 0) ScaleFacBG_HER[i] = I_HER * iP_HER / (m_input_I_HER[0] * m_input_P_HER[0]);
1252  }
1253  }
1254  }
1255  if (m_beast.SKB_LER_pressures != 0 && m_input_GasCorrection[1] != 1) {
1256  for (int i = 0; i < (int)m_beast.SKB_LER_pressures->size(); i++) {
1257  ScaleFacBG_LER[i] = 0;
1258  double iP_LER = 0;
1259  iP_LER = m_beast.SKB_LER_pressures->at(i) * 0.00750062 * 1e9 * m_input_GasCorrection[1];
1260  if (m_input_P_LER[1] > 0) iP_LER += gRandom->Gaus(0, m_input_P_LER[1]);
1261  if (iP_LER < 0 || iP_LER > 260.) iP_LER = 0;
1262  if (I_LER > 0 && iP_LER > 0) {
1263  if (m_input_BGSol == 0
1264  && bunch_nb_LER > 0) ScaleFacBG_LER[i] = I_LER * iP_LER / (m_input_I_LER[0] * m_input_P_LER[0]) / bunch_nb_LER *
1265  m_input_bunchNb_LER[0];
1266  if (m_input_BGSol == 1 && bunch_nb_LER > 0) ScaleFacBG_LER[i] = I_LER * iP_LER / (m_input_I_LER[0] * m_input_P_LER[0]);
1267  }
1268  }
1269  }
1270 
1271  if (m_beast.SKB_HER_pressures_corrected != 0 && m_input_GasCorrection[0] == 1) {
1272  for (int i = 0; i < (int)m_beast.SKB_HER_pressures_corrected->size(); i++) {
1273  ScaleFacBG_HER[i] = 0;
1274  double iP_HER = 0;
1275  iP_HER = m_beast.SKB_HER_pressures_corrected->at(i) * 0.00750062 * 1e9 * m_input_GasCorrection[0];
1276  if (m_input_P_HER[1] > 0) iP_HER += gRandom->Gaus(0, m_input_P_HER[1]);
1277  if (iP_HER < 0 || iP_HER > 260.) iP_HER = 0;
1278  if (I_HER > 0 && iP_HER > 0) {
1279  if (m_input_BGSol == 0
1280  && bunch_nb_HER > 0) ScaleFacBG_HER[i] = I_HER * iP_HER / (m_input_I_HER[0] * m_input_P_HER[0]) / bunch_nb_HER *
1281  m_input_bunchNb_HER[0];
1282  if (m_input_BGSol == 1 && bunch_nb_HER > 0) ScaleFacBG_HER[i] = I_HER * iP_HER / (m_input_I_HER[0] * m_input_P_HER[0]);
1283  }
1284  }
1285  }
1286  if (m_beast.SKB_LER_pressures_corrected != 0 && m_input_GasCorrection[1] == 1) {
1287  for (int i = 0; i < (int)m_beast.SKB_LER_pressures_corrected->size(); i++) {
1288  ScaleFacBG_LER[i] = 0;
1289  double iP_LER = 0;
1290  iP_LER = m_beast.SKB_LER_pressures_corrected->at(i) * 0.00750062 * 1e9 * m_input_GasCorrection[1];
1291  if (m_input_P_LER[1] > 0) iP_LER += gRandom->Gaus(0, m_input_P_LER[1]);
1292  if (iP_LER < 0 || iP_LER > 260.) iP_LER = 0;
1293  if (I_LER > 0 && iP_LER > 0) {
1294  if (m_input_BGSol == 0
1295  && bunch_nb_LER > 0) ScaleFacBG_LER[i] = I_LER * iP_LER / (m_input_I_LER[0] * m_input_P_LER[0]) / bunch_nb_LER *
1296  m_input_bunchNb_LER[0];
1297  if (m_input_BGSol == 1 && bunch_nb_LER > 0) ScaleFacBG_LER[i] = I_LER * iP_LER / (m_input_I_LER[0] * m_input_P_LER[0]);
1298  }
1299  }
1300  }
1301  //Calculate Touschek scaling factor: Touschek \propo I^2 / (bunch_nb x sigma_y) => (I^2/(bunch_nb x sigma_y))^data / (I^2/(bunch_nb x sigma_y))^simu
1302  double ScaleFacTo_HER = 0;
1303  double ScaleFacTo_LER = 0;
1304  if (bunch_nb_LER > 0 && sigma_y_LER > 0 && I_LER > 0) {
1305  if (m_input_ToSol == 0) ScaleFacTo_LER = TMath::Power(I_LER / m_input_I_LER[0],
1306  2) / (bunch_nb_LER / m_input_bunchNb_LER[0]) / (sigma_y_LER / m_input_sigma_y_LER[0]);
1307  else if (m_input_ToSol == 1) ScaleFacTo_LER = TMath::Power(I_LER / m_input_I_LER[0],
1308  2) / (bunch_nb_LER / m_input_bunchNb_LER[0]) /
1309  (sqrt(sigma_y_LER * sigma_x_LER / m_input_sigma_y_LER[0] / m_input_sigma_x_LER[0]));
1310  }
1311  if (bunch_nb_HER > 0 && sigma_y_HER > 0 && I_HER > 0) {
1312  if (m_input_ToSol == 0) ScaleFacTo_HER = TMath::Power(I_HER / m_input_I_HER[0],
1313  2) / (bunch_nb_HER / m_input_bunchNb_HER[0]) / (sigma_y_HER / m_input_sigma_y_HER[0]);
1314  else if (m_input_ToSol == 1) ScaleFacTo_HER = TMath::Power(I_HER / m_input_I_HER[0],
1315  2) / (bunch_nb_HER / m_input_bunchNb_HER[0]) /
1316  (sqrt(sigma_y_HER * sigma_x_HER / m_input_sigma_y_HER[0] / m_input_sigma_x_HER[0]));
1317  }
1318 
1319  //Save reweight information for Touschek and BeamGas-Coulomb & -Brems
1320  m_beast.mc_reweight_LERT.push_back(ScaleFacTo_LER);
1321  m_beast.mc_reweight_HERT.push_back(ScaleFacTo_HER);
1322  for (int i = 0; i < 12; i ++) {
1323  m_beast.mc_reweight_HERC.push_back(Zeff_HC * ScaleFacBG_HER[i]);
1324  m_beast.mc_reweight_HERB.push_back(Zeff_HB * ScaleFacBG_HER[i]);
1325  if (Zeff_LC != 1)
1326  m_beast.mc_reweight_LERC.push_back(Zeff_LC * ScaleFacBG_LER[i]);
1327  else
1328  m_beast.mc_reweight_LERC.push_back(m_input_Z_scaling[1] * ScaleFacBG_LER[i]);
1329  if (Zeff_LB != 1)
1330  m_beast.mc_reweight_LERB.push_back(Zeff_LB * ScaleFacBG_LER[i]);
1331  else
1332  m_beast.mc_reweight_LERB.push_back(m_input_Z_scaling[3] * ScaleFacBG_LER[i]);
1333  }
1334  //Only Touschek LER + HER
1335  if (m_input_part == 0) {
1336  ScaleFacBGav_HER = 0;
1337  ScaleFacBGav_LER = 0;
1338  for (int i = 0; i < 12; i ++) {
1339  ScaleFacBG_HER[i] = 0;
1340  ScaleFacBG_LER[i] = 0;
1341  }
1342  }
1343  //Only Beam Gas
1344  if (m_input_part == 1) {
1345  ScaleFacTo_HER = 0;
1346  ScaleFacTo_LER = 0;
1347  }
1348  /*
1349  //Only Beam Gas - LER Brems
1350  if (m_input_part == 2) {
1351  ScaleFacTo_HER = 0;
1352  ScaleFacTo_LER = 0;
1353  Zeff_LC = 0;
1354  Zeff_HB = 0;
1355  Zeff_HC = 0;
1356  }
1357  //Only Beam Gas - LER Coulomb
1358  if (m_input_part == 3) {
1359  ScaleFacTo_HER = 0;
1360  ScaleFacTo_LER = 0;
1361  Zeff_LB = 0;
1362  Zeff_HB = 0;
1363  Zeff_HC = 0;
1364  }
1365  //Only Beam Gas - HER Brems
1366  if (m_input_part == 4) {
1367  ScaleFacTo_HER = 0;
1368  ScaleFacTo_LER = 0;
1369  Zeff_HC = 0;
1370  Zeff_LB = 0;
1371  Zeff_LC = 0;
1372  }
1373  //Only Beam Gas - HER Coulomb
1374  if (m_input_part == 5) {
1375  ScaleFacTo_HER = 0;
1376  ScaleFacTo_LER = 0;
1377  Zeff_HB = 0;
1378  Zeff_LB = 0;
1379  Zeff_LC = 0;
1380  }
1381  */
1382  //Scale LER SAD RLR
1383  for (int i = 0; i < (int)m_input_LT_SAD_RLR.size(); i++) {
1384  float LBG = m_input_LB_SAD_RLR_av[i] * Zeff_LB + m_input_LC_SAD_RLR_av[i] * Zeff_LC;
1385  float BG = LBG * ScaleFacBGav_LER;
1386  float To = ScaleFacTo_LER * m_input_LT_SAD_RLR[i];
1387  //if (TMath::IsNaN(To)) To = 0;
1388  //if (TMath::IsNaN(BG)) BG = 0;
1389  m_beast.SAD_LER_RLR_av.push_back(BG + To);
1390  m_beast.SAD_LER_lifetime_av.push_back(Nb_LER / (BG + To) * 1e-9 / 60. * bunch_nb_LER);
1391  BG = 0;
1392  LBG = 0;
1393  for (int j = 0; j < 12; j++) {
1394  LBG = 0;
1395  if (m_input_LB_SAD_RLR.size() > 0) {
1396  LBG = m_input_LB_SAD_RLR[j] * Zeff_LB + m_input_LC_SAD_RLR[j] * Zeff_LC;
1397  BG += LBG * ScaleFacBG_LER[j];
1398  }
1399  }
1400  //if (TMath::IsNaN(To)) To = 0;
1401  //if (TMath::IsNaN(BG)) BG = 0;
1402  m_beast.SAD_LER_RLR.push_back(BG + To);
1403  m_beast.SAD_LER_lifetime.push_back(Nb_LER / (BG + To) * 1e-9 / 60. * bunch_nb_LER);
1404  }
1405 
1406  //Scale HER SAD RLR
1407  for (int i = 0; i < (int)m_input_HT_SAD_RLR.size(); i++) {
1408  float HBG = m_input_HB_SAD_RLR_av[i] + m_input_HC_SAD_RLR_av[i];
1409  float BG = HBG * ScaleFacBGav_HER;
1410  float To = ScaleFacTo_HER * m_input_HT_SAD_RLR[i];
1411  //if (TMath::IsNaN(To)) To = 0;
1412  //if (TMath::IsNaN(BG)) BG = 0;
1413  m_beast.SAD_HER_RLR_av.push_back(BG + To);
1414  m_beast.SAD_HER_lifetime_av.push_back(Nb_HER / (BG + To) * 1e-9 / 60. * bunch_nb_HER);
1415  BG = 0;
1416  HBG = 0;
1417  for (int j = 0; j < 12; j++) {
1418  HBG = 0;
1419  if (m_input_HB_SAD_RLR.size() > 0) {
1420  HBG = m_input_HB_SAD_RLR[j] + m_input_HC_SAD_RLR[j];
1421  BG += HBG * ScaleFacBG_HER[j];
1422  }
1423  }
1424  //if (TMath::IsNaN(To)) To = 0;
1425  //if (TMath::IsNaN(BG)) BG = 0;
1426  m_beast.SAD_HER_RLR.push_back(BG + To);
1427  m_beast.SAD_HER_lifetime.push_back(Nb_HER / (BG + To) * 1e-9 / 60. * bunch_nb_HER);
1428  }
1429 
1430  //Scale DIA
1431  for (int i = 0; i < (int)m_input_LT_DIA_dose.size(); i++) {
1432  double LBG = m_input_LB_DIA_dose_av[i] + m_input_LC_DIA_dose_av[i];
1433  double HBG = m_input_HB_DIA_dose_av[i] + m_input_HC_DIA_dose_av[i];
1434  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1435  double To = ScaleFacTo_LER * m_input_LT_DIA_dose[i] + ScaleFacTo_HER * m_input_HT_DIA_dose[i];
1436  /*cout << "ch # " << i
1437  << " av LB " << m_input_LB_DIA_dose_av[i] << " LC " << m_input_LC_DIA_dose_av[i]
1438  << " HB " << m_input_HB_DIA_dose_av[i] << " HC " << m_input_HC_DIA_dose_av[i]
1439  << " LT " << m_input_LT_DIA_dose[i] << " HT " << m_input_HT_DIA_dose[i]
1440  << endl;*/
1441  //if (TMath::IsNaN(To)) To = 0;
1442  //if (TMath::IsNaN(BG)) BG = 0;
1443  m_beast.DIA_dose_av.push_back(BG + To);
1444  BG = 0; LBG = 0; HBG = 0;
1445  for (int j = 0; j < 12; j++) {
1446  LBG = 0; HBG = 0;
1447  if (m_input_LB_DIA_dose[j].size() > 0) {
1448  //LBG = m_input_LB_DIA_dose[j][i] + m_input_LC_DIA_dose[j][i];
1449  HBG = m_input_HB_DIA_dose[j][i] + m_input_HC_DIA_dose[j][i];
1450  LBG = m_input_LB_DIA_dose[j][i] * Zeff_LB + m_input_LC_DIA_dose[j][i] * Zeff_LC;
1451  /*cout << "section " << j
1452  << " LB " << m_input_LB_DIA_dose[j][i] << " LC " << m_input_LC_DIA_dose[j][i] << " HB " << m_input_HB_DIA_dose[j][i] << " HC " <<
1453  m_input_HC_DIA_dose[j][i] << endl;*/
1454  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1455  }
1456  }
1457  //if (TMath::IsNaN(To)) To = 0;
1458  //if (TMath::IsNaN(BG)) BG = 0;
1459  m_beast.DIA_dose.push_back(BG + To);
1460  }
1461 
1462  //Scale PIN
1463  for (int i = 0; i < (int)m_input_LT_PIN_dose.size(); i++) {
1464  double LBG = m_input_LB_PIN_dose_av[i] + m_input_LC_PIN_dose_av[i];
1465  double HBG = m_input_HB_PIN_dose_av[i] + m_input_HC_PIN_dose_av[i];
1466  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1467  double To = ScaleFacTo_LER * m_input_LT_PIN_dose[i] + ScaleFacTo_HER * m_input_HT_PIN_dose[i];
1468  //if (TMath::IsNaN(To)) To = 0;
1469  //if (TMath::IsNaN(BG)) BG = 0;
1470  m_beast.PIN_dose_av.push_back(BG + To);
1471  BG = 0; LBG = 0; HBG = 0;
1472  for (int j = 0; j < 12; j++) {
1473  LBG = 0; HBG = 0;
1474  if (m_input_LB_PIN_dose[j].size() > 0) {
1475  //LBG = m_input_LB_PIN_dose[j][i] + m_input_LC_PIN_dose[j][i];
1476  HBG = m_input_HB_PIN_dose[j][i] + m_input_HC_PIN_dose[j][i];
1477  LBG = m_input_LB_PIN_dose[j][i] * Zeff_LB + m_input_LC_PIN_dose[j][i] * Zeff_LC;
1478  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1479  }
1480  }
1481  //if (TMath::IsNaN(To)) To = 0;
1482  //if (TMath::IsNaN(BG)) BG = 0;
1483  m_beast.PIN_dose.push_back(BG + To);
1484  }
1485 
1486  //Scale DOSI
1487  for (int i = 0; i < (int)m_input_LT_DOSI.size(); i++) {
1488  //cout << "LB : " << m_input_LB_DOSI_av[i] << " LC " << m_input_LC_DOSI_av[i] << endl;
1489  //cout << "HB : " << m_input_HB_DOSI_av[i] << " HC " << m_input_HC_DOSI_av[i] << endl;
1490  //cout << "HT : " << m_input_HT_DOSI[i] << " LT " << m_input_LT_DOSI[i] << endl;
1491  double LBG = m_input_LB_DOSI_av[i] + m_input_LC_DOSI_av[i];
1492  double HBG = m_input_HB_DOSI_av[i] + m_input_HC_DOSI_av[i];
1493  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1494  double To = ScaleFacTo_LER * m_input_LT_DOSI[i] + ScaleFacTo_HER * m_input_HT_DOSI[i];
1495  //if (TMath::IsNaN(To)) To = 0;
1496  //if (TMath::IsNaN(BG)) BG = 0;
1497  m_beast.DOSI_av.push_back(BG + To);
1498  BG = 0; LBG = 0; HBG = 0;
1499  for (int j = 0; j < 12; j++) {
1500  LBG = 0; HBG = 0;
1501  if (m_input_LB_DOSI[j].size() > 0) {
1502  //LBG = m_input_LB_DOSI[j][i] + m_input_LC_DOSI[j][i];
1503  HBG = m_input_HB_DOSI[j][i] + m_input_HC_DOSI[j][i];
1504  LBG = m_input_LB_DOSI[j][i] * Zeff_LB + m_input_LC_DOSI[j][i] * Zeff_LC;
1505  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1506  }
1507  }
1508  //if (TMath::IsNaN(To)) To = 0;
1509  //if (TMath::IsNaN(BG)) BG = 0;
1510  m_beast.DOSI.push_back(BG + To);
1511  }
1512 
1513  //Scale BGO
1514  for (int i = 0; i < (int)m_input_LT_BGO_dose.size(); i++) {
1515  double LBG = m_input_LB_BGO_dose_av[i] + m_input_LC_BGO_dose_av[i];
1516  double HBG = m_input_HB_BGO_dose_av[i] + m_input_HC_BGO_dose_av[i];
1517  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1518  double To = ScaleFacTo_LER * m_input_LT_BGO_dose[i] + ScaleFacTo_HER * m_input_HT_BGO_dose[i];
1519  //if (TMath::IsNaN(To)) To = 0;
1520  //if (TMath::IsNaN(BG)) BG = 0;
1521  m_beast.BGO_energy_av.push_back(BG + To);
1522  BG = 0; LBG = 0; HBG = 0;
1523  for (int j = 0; j < 12; j++) {
1524  LBG = 0; HBG = 0;
1525  if (m_input_LB_BGO_dose[j].size() > 0) {
1526  //LBG = m_input_LB_BGO_dose[j][i] + m_input_LC_BGO_dose[j][i];
1527  HBG = m_input_HB_BGO_dose[j][i] + m_input_HC_BGO_dose[j][i];
1528  LBG = m_input_LB_BGO_dose[j][i] * Zeff_LB + m_input_LC_BGO_dose[j][i] * Zeff_LC;
1529  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1530  }
1531  }
1532  //if (TMath::IsNaN(To)) To = 0;
1533  //if (TMath::IsNaN(BG)) BG = 0;
1534  m_beast.BGO_energy.push_back(BG + To);
1535  }
1536  int he3order[4];
1537  if (m_beast.ts > 1464868800) {
1538  he3order[0] = 0;
1539  he3order[1] = 3;
1540  he3order[2] = 2;
1541  he3order[3] = 1;
1542  } else {
1543  he3order[0] = 0;
1544  he3order[1] = 1;
1545  he3order[2] = 2;
1546  he3order[3] = 3;
1547  }
1548  //Scale HE3
1549  for (int i = 0; i < (int)m_input_LT_HE3_rate.size(); i++) {
1550  double LBG = m_input_LB_HE3_rate_av[he3order[i]] + m_input_LC_HE3_rate_av[he3order[i]];
1551  double HBG = m_input_HB_HE3_rate_av[he3order[i]] + m_input_HC_HE3_rate_av[he3order[i]];
1552  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1553  double To = ScaleFacTo_LER * m_input_LT_HE3_rate[he3order[i]] + ScaleFacTo_HER * m_input_HT_HE3_rate[he3order[i]];
1554  //if (TMath::IsNaN(To)) To = 0;
1555  //if (TMath::IsNaN(BG)) BG = 0;
1556  m_beast.HE3_rate_av.push_back(BG + To);
1557  BG = 0; LBG = 0; HBG = 0;
1558  for (int j = 0; j < 12; j++) {
1559  LBG = 0; HBG = 0;
1560  if (m_input_LB_HE3_rate[j].size() > 0) {
1561  //LBG = m_input_LB_HE3_rate[j][he3order[i]] + m_input_LC_HE3_rate[j][he3order[i]];
1562  HBG = m_input_HB_HE3_rate[j][he3order[i]] + m_input_HC_HE3_rate[j][he3order[i]];
1563  LBG = m_input_LB_HE3_rate[j][he3order[i]] * Zeff_LB + m_input_LC_HE3_rate[j][he3order[i]] * Zeff_LC;
1564  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1565  }
1566  }
1567  //if (TMath::IsNaN(To)) To = 0;
1568  //if (TMath::IsNaN(BG)) BG = 0;
1569  m_beast.HE3_rate.push_back(BG + To);
1570  }
1571 
1572  //Scale TPC
1573  for (int i = 0; i < (int)m_input_LT_TPC_rate.size(); i++) {
1574  double LBG = m_input_LB_TPC_rate_av[i] + m_input_LC_TPC_rate_av[i];
1575  double HBG = m_input_HB_TPC_rate_av[i] + m_input_HC_TPC_rate_av[i];
1576  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1577  double To = ScaleFacTo_LER * m_input_LT_TPC_rate[i] + ScaleFacTo_HER * m_input_HT_TPC_rate[i];
1578  int tpc_ch = (int)(i / 5);
1579  int n_type = i - 5 * tpc_ch;
1580  m_beast.TPC_rate_av[tpc_ch][n_type] = (BG + To);
1581  BG = 0; LBG = 0; HBG = 0;
1582  for (int j = 0; j < 12; j++) {
1583  LBG = 0; HBG = 0;
1584  if (m_input_LB_TPC_rate[j].size() > 0) {
1585  //LBG = m_input_LB_TPC_rate[j][i] + m_input_LC_TPC_rate[j][i];
1586  HBG = m_input_HB_TPC_rate[j][i] + m_input_HC_TPC_rate[j][i];
1587  LBG = m_input_LB_TPC_rate[j][i] * Zeff_LB + m_input_LC_TPC_rate[j][i] * Zeff_LC;
1588  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1589  }
1590  }
1591  //m_beast.TPC_rate[tpc_ch][n_type] = (BG + To);
1592  }
1593 
1594 
1595  //Scale TPC_dose
1596  for (int i = 0; i < (int)m_input_LT_TPC_dose.size(); i++) {
1597  double LBG = m_input_LB_TPC_dose_av[i] + m_input_LC_TPC_dose_av[i];
1598  double HBG = m_input_HB_TPC_dose_av[i] + m_input_HC_TPC_dose_av[i];
1599  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1600  double To = ScaleFacTo_LER * m_input_LT_TPC_dose[i] + ScaleFacTo_HER * m_input_HT_TPC_dose[i];
1601  //int tpc_ch = (int)(i / 5);
1602  //int n_type = i - 5 * tpc_ch;
1603  //m_beast.TPC_dose_av[tpc_ch][n_type] = (BG + To);
1604  m_beast.TPC_dose_av.push_back(BG + To);
1605  BG = 0; LBG = 0; HBG = 0;
1606  for (int j = 0; j < 12; j++) {
1607  LBG = 0; HBG = 0;
1608  if (m_input_LB_TPC_dose[j].size() > 0) {
1609  //LBG = m_input_LB_TPC_dose[j][i] + m_input_LC_TPC_dose[j][i];
1610  HBG = m_input_HB_TPC_dose[j][i] + m_input_HC_TPC_dose[j][i];
1611  LBG = m_input_LB_TPC_dose[j][i] * Zeff_LB + m_input_LC_TPC_dose[j][i] * Zeff_LC;
1612  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1613  }
1614  }
1615  //m_beast.TPC_dose[tpc_ch][n_type] = (BG + To);
1616  }
1617 
1618  //Scale TPC_angular
1619  for (int i = 0; i < (int)m_input_LT_TPC_angular_rate.size(); i++) {
1620  double LBG = m_input_LB_TPC_angular_rate_av[i] + m_input_LC_TPC_angular_rate_av[i];
1621  double HBG = m_input_HB_TPC_angular_rate_av[i] + m_input_HC_TPC_angular_rate_av[i];
1622  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1623  double To = ScaleFacTo_LER * m_input_LT_TPC_angular_rate[i] + ScaleFacTo_HER * m_input_HT_TPC_angular_rate[i];
1624  int tpc_ch = (int)(i / (9 * 18));
1625  int angle = i - (9 * 18) * tpc_ch;
1626  int i_theta = (int)(angle / 18);
1627  int i_phi = angle - 9 * i_theta;
1628  m_beast.TPC_angular_rate_av[tpc_ch][i_theta][i_phi] = (BG + To);
1629  BG = 0; LBG = 0; HBG = 0;
1630  for (int j = 0; j < 12; j++) {
1631  LBG = 0; HBG = 0;
1632  if (m_input_LB_TPC_angular_rate[j].size() > 0) {
1633  //LBG = m_input_LB_TPC_angular_rate[j][i] + m_input_LC_TPC_angular_rate[j][i];
1634  HBG = m_input_HB_TPC_angular_rate[j][i] + m_input_HC_TPC_angular_rate[j][i];
1635  LBG = m_input_LB_TPC_angular_rate[j][i] * Zeff_LB + m_input_LC_TPC_angular_rate[j][i] * Zeff_LC;
1636  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1637  }
1638  }
1639  //m_beast.TPC_angular_rate[tpc_ch][n_type] = (BG + To);
1640  }
1641 
1642  //Scale TPC_angular_dose
1643  for (int i = 0; i < (int)m_input_LT_TPC_angular_dose.size(); i++) {
1644  double LBG = m_input_LB_TPC_angular_dose_av[i] + m_input_LC_TPC_angular_dose_av[i];
1645  double HBG = m_input_HB_TPC_angular_dose_av[i] + m_input_HC_TPC_angular_dose_av[i];
1646  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1647  double To = ScaleFacTo_LER * m_input_LT_TPC_angular_dose[i] + ScaleFacTo_HER * m_input_HT_TPC_angular_dose[i];
1648  int tpc_ch = (int)(i / (9 * 18));
1649  int angle = i - (9 * 18) * tpc_ch;
1650  int i_theta = (int)(angle / 18);
1651  int i_phi = angle - 9 * i_theta;
1652  m_beast.TPC_angular_dose_av[tpc_ch][i_theta][i_phi] = (BG + To);
1653  BG = 0; LBG = 0; HBG = 0;
1654  for (int j = 0; j < 12; j++) {
1655  LBG = 0; HBG = 0;
1656  if (m_input_LB_TPC_angular_dose[j].size() > 0) {
1657  //LBG = m_input_LB_TPC_angular_dose[j][i] + m_input_LC_TPC_angular_dose[j][i];
1658  HBG = m_input_HB_TPC_angular_dose[j][i] + m_input_HC_TPC_angular_dose[j][i];
1659  LBG = m_input_LB_TPC_angular_dose[j][i] * Zeff_LB + m_input_LC_TPC_angular_dose[j][i] * Zeff_LC;
1660  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1661  }
1662  }
1663  //m_beast.TPC_angular_dose[tpc_ch][n_type] = (BG + To);
1664  }
1665 
1666 
1667  //Scale CLAWS
1668  for (int i = 0; i < (int)m_input_LT_CLAWS_rate.size(); i++) {
1669  double LBG = m_input_LB_CLAWS_rate_av[i] + m_input_LC_CLAWS_rate_av[i];
1670  double HBG = m_input_HB_CLAWS_rate_av[i] + m_input_HC_CLAWS_rate_av[i];
1671  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1672  double To = ScaleFacTo_LER * m_input_LT_CLAWS_rate[i] + ScaleFacTo_HER * m_input_HT_CLAWS_rate[i];
1673  //if (TMath::IsNaN(To)) To = 0;
1674  //if (TMath::IsNaN(BG)) BG = 0;
1675  m_beast.CLAWS_rate_av.push_back(BG + To);
1676  BG = 0; LBG = 0; HBG = 0;
1677  for (int j = 0; j < 12; j++) {
1678  LBG = 0; HBG = 0;
1679  if (m_input_LB_CLAWS_rate[j].size() > 0) {
1680  //LBG = m_input_LB_CLAWS_rate[j][i] + m_input_LC_CLAWS_rate[j][i];
1681  HBG = m_input_HB_CLAWS_rate[j][i] + m_input_HC_CLAWS_rate[j][i];
1682  LBG = m_input_LB_CLAWS_rate[j][i] * Zeff_LB + m_input_LC_CLAWS_rate[j][i] * Zeff_LC;
1683  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1684  }
1685  }
1686  //if (TMath::IsNaN(To)) To = 0;
1687  //if (TMath::IsNaN(BG)) BG = 0;
1688  m_beast.CLAWS_rate.push_back(BG + To);
1689  }
1690 
1691  //Scale QCSS
1692  for (int i = 0; i < (int)m_input_LT_QCSS_rate.size(); i++) {
1693  double LBG = m_input_LB_QCSS_rate_av[i] + m_input_LC_QCSS_rate_av[i];
1694  double HBG = m_input_HB_QCSS_rate_av[i] + m_input_HC_QCSS_rate_av[i];
1695  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1696  double To = ScaleFacTo_LER * m_input_LT_QCSS_rate[i] + ScaleFacTo_HER * m_input_HT_QCSS_rate[i];
1697  //if (TMath::IsNaN(To)) To = 0;
1698  //if (TMath::IsNaN(BG)) BG = 0;
1699  m_beast.QCSS_rate_av.push_back(BG + To);
1700  BG = 0; LBG = 0; HBG = 0;
1701  for (int j = 0; j < 12; j++) {
1702  LBG = 0; HBG = 0;
1703  if (m_input_LB_QCSS_rate[j].size() > 0) {
1704  //LBG = m_input_LB_QCSS_rate[j][i] + m_input_LC_QCSS_rate[j][i];
1705  HBG = m_input_HB_QCSS_rate[j][i] + m_input_HC_QCSS_rate[j][i];
1706  LBG = m_input_LB_QCSS_rate[j][i] * Zeff_LB + m_input_LC_QCSS_rate[j][i] * Zeff_LC;
1707  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1708  }
1709  }
1710  //if (TMath::IsNaN(To)) To = 0;
1711  //if (TMath::IsNaN(BG)) BG = 0;
1712  m_beast.QCSS_rate.push_back(BG + To);
1713  }
1714 
1715  //Scale CSI
1716  for (int i = 0; i < (int)m_input_LT_CSI_dose.size(); i++) {
1717  double LBG = m_input_LB_CSI_dose_av[i] + m_input_LC_CSI_dose_av[i];
1718  double HBG = m_input_HB_CSI_dose_av[i] + m_input_HC_CSI_dose_av[i];
1719  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1720  double To = ScaleFacTo_LER * m_input_LT_CSI_dose[i] + ScaleFacTo_HER * m_input_HT_CSI_dose[i];
1721  //if (TMath::IsNaN(To)) To = 0;
1722  //if (TMath::IsNaN(BG)) BG = 0;
1723  m_beast.CSI_sumE_av.push_back(BG + To);
1724  BG = 0; LBG = 0; HBG = 0;
1725  for (int j = 0; j < 12; j++) {
1726  LBG = 0; HBG = 0;
1727  if (m_input_LB_CSI_dose[j].size() > 0) {
1728  //LBG = m_input_LB_CSI_dose[j][i] + m_input_LC_CSI_dose[j][i];
1729  HBG = m_input_HB_CSI_dose[j][i] + m_input_HC_CSI_dose[j][i];
1730  LBG = m_input_LB_CSI_dose[j][i] * Zeff_LB + m_input_LC_CSI_dose[j][i] * Zeff_LC;
1731  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1732  }
1733  }
1734  //if (TMath::IsNaN(To)) To = 0;
1735  //if (TMath::IsNaN(BG)) BG = 0;
1736  m_beast.CSI_sumE.push_back(BG + To);
1737  }
1738  for (int i = 0; i < (int)m_input_LT_CSI_dose_binE.size(); i++) {
1739  double LBG = m_input_LB_CSI_dose_binE_av[i] + m_input_LC_CSI_dose_binE_av[i];
1740  double HBG = m_input_HB_CSI_dose_binE_av[i] + m_input_HC_CSI_dose_binE_av[i];
1741  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1742  double To = ScaleFacTo_LER * m_input_LT_CSI_dose_binE[i] + ScaleFacTo_HER * m_input_HT_CSI_dose_binE[i];
1743  //if (TMath::IsNaN(To)) To = 0;
1744  //if (TMath::IsNaN(BG)) BG = 0;
1745  m_beast.CSI_Ebin_av.push_back(BG + To);
1746  BG = 0; LBG = 0; HBG = 0;
1747  for (int j = 0; j < 12; j++) {
1748  LBG = 0; HBG = 0;
1749  if (m_input_LB_CSI_dose_binE[j].size() > 0) {
1750  //LBG = m_input_LB_CSI_dose_binE[j][i] + m_input_LC_CSI_dose_binE[j][i];
1751  HBG = m_input_HB_CSI_dose_binE[j][i] + m_input_HC_CSI_dose_binE[j][i];
1752  LBG = m_input_LB_CSI_dose_binE[j][i] * Zeff_LB + m_input_LC_CSI_dose_binE[j][i] * Zeff_LC;
1753  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1754  }
1755  }
1756  //if (TMath::IsNaN(To)) To = 0;
1757  //if (TMath::IsNaN(BG)) BG = 0;
1758  m_beast.CSI_Ebin.push_back(BG + To);
1759  }
1760  for (int i = 0; i < (int)m_input_LT_CSI_rate.size(); i++) {
1761  double LBG = m_input_LB_CSI_rate_av[i] + m_input_LC_CSI_rate_av[i];
1762  double HBG = m_input_HB_CSI_rate_av[i] + m_input_HC_CSI_rate_av[i];
1763  double BG = LBG * ScaleFacBGav_LER + HBG * ScaleFacBGav_HER;
1764  double To = ScaleFacTo_LER * m_input_LT_CSI_rate[i] + ScaleFacTo_HER * m_input_HT_CSI_rate[i];
1765  //if (TMath::IsNaN(To)) To = 0;
1766  //if (TMath::IsNaN(BG)) BG = 0;
1767  m_beast.CSI_hitRate_av.push_back(BG + To);
1768  BG = 0; LBG = 0; HBG = 0;
1769  for (int j = 0; j < 12; j++) {
1770  LBG = 0; HBG = 0;
1771  if (m_input_LB_CSI_rate[j].size() > 0) {
1772  //LBG = m_input_LB_CSI_rate[j][i] + m_input_LC_CSI_rate[j][i];
1773  HBG = m_input_HB_CSI_rate[j][i] + m_input_HC_CSI_rate[j][i];
1774  LBG = m_input_LB_CSI_rate[j][i] * Zeff_LB + m_input_LC_CSI_rate[j][i] * Zeff_LC;
1775  BG += LBG * ScaleFacBG_LER[j] + HBG * ScaleFacBG_HER[j];
1776  }
1777  }
1778  //if (TMath::IsNaN(To)) To = 0;
1779  //if (TMath::IsNaN(BG)) BG = 0;
1780  m_beast.CSI_hitRate.push_back(BG + To);
1781  }
1782 
1783  m_treeBEAST->Fill();
1784 
1785  // set event metadata
1786  //evtMetaData->setEvent(m_eventCount);
1787  //evtMetaData->setRun(m_run);
1788  //evtMetaData->setExperiment(m_exp);
1789 
1790  m_eventCount++;
1791 
1792  }
1793 
1794 
1795  void NtuplePhase1_v6Module::endRun()
1796  {
1797  }
1798 
1799  void NtuplePhase1_v6Module::terminate()
1800  {
1801  delete m_tree;
1802  m_file->cd();
1803  m_treeBEAST->Write();
1804  m_treeTruth->Write();
1805  m_file->Close();
1806  }
1807 
1808  void NtuplePhase1_v6Module::printModuleParams() const
1809  {
1810  }
1811 
1813 } // end Belle2 namespace
Belle2::Module
Base class for Modules.
Definition: Module.h:72
Belle2::NtuplePhase1_v6Module
Read SKB PVs, simulated measurements of BEAST sensors, and write scaled simulated Ntuple in BEAST pha...
Definition: NtuplePhase1_v6Module.h:28
Belle2::StoreAccessorBase::registerInDataStore
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Definition: StoreAccessorBase.h:52
Belle2::StoreAccessorBase::create
bool create(bool replace=false)
Create a default object in the data store.
Definition: StoreAccessorBase.h:107
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2::RootIOUtilities::c_treeNames
const std::string c_treeNames[]
Names of trees.
Definition: RootIOUtilities.cc:27
Belle2::RootIOUtilities::expandWordExpansions
std::vector< std::string > expandWordExpansions(const std::vector< std::string > &filenames)
Performs wildcard expansion using wordexp(), returns matches.
Definition: RootIOUtilities.cc:114
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17