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