Belle II Software development
eclLeakageCollectorModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9/* Own header. */
10#include <ecl/modules/eclLeakageCollector/eclLeakageCollectorModule.h>
11
12/* ECL headers. */
13#include <ecl/dataobjects/ECLShower.h>
14#include <ecl/dataobjects/ECLCalDigit.h>
15#include <ecl/geometry/ECLLeakagePosition.h>
16#include <ecl/dataobjects/ECLElementNumbers.h>
17
18/* Basf2 headers. */
19#include <framework/dataobjects/EventMetaData.h>
20#include <framework/geometry/VectorUtil.h>
21#include <mdst/dataobjects/MCParticle.h>
22
23/* Root headers. */
24#include <Math/Vector3D.h>
25#include <Math/VectorUtil.h>
26#include <TTree.h>
27
28/* C++ headers. */
29#include <cmath>
30#include <iostream>
31
32using namespace Belle2;
33using namespace ECL;
34
35//-----------------------------------------------------------------
36// Register the Module
37//-----------------------------------------------------------------
38REG_MODULE(eclLeakageCollector);
39
40//-----------------------------------------------------------------
41// Implementation
42//-----------------------------------------------------------------
43
44//-----------------------------------------------------------------
46 m_mcParticleArray("MCParticles"),
47 m_evtMetaData("EventMetaData")
48{
50 setDescription("Store quantities from single photon MC used to calculated ECL energy leakage corrections");
52 addParam("position_bins", m_position_bins, "number of crystal subdivisions in theta and phi", 29);
53 addParam("number_energies", m_number_energies, "number of generated energies", 8);
54 addParam("energies_forward", m_energies_forward, "generated photon energies, forward", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.483, 1.166, 2.816, 6.800});
55 addParam("energies_barrel", m_energies_barrel, "generated photon energies, barrel", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.458, 1.049, 2.402, 5.500});
56 addParam("energies_backward", m_energies_backward, "generated photon energies, backward", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.428, 0.917, 1.962, 4.200});
57 addParam("showerArrayName", m_showerArrayName, "name of ECLShower data object", std::string("ECLTrimmedShowers"));
58}
59
60
61//-----------------------------------------------------------------
63{
64 //-----------------------------------------------------------------
65 //..Parameters and other basic info
66 B2INFO("eclLeakageCollector: Experiment = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
67
68 //..Check that input parameters are consistent
69 const int n_e_forward = m_energies_forward.size();
70 const int n_e_barrel = m_energies_barrel.size();
71 const int n_e_backward = m_energies_backward.size();
72 if (n_e_forward != m_number_energies or n_e_barrel != m_number_energies or n_e_backward != m_number_energies) {
73 B2FATAL("eclLeakageCollector: length of energy vectors inconsistent with parameter number_energies: " << n_e_forward << " " <<
74 n_e_barrel << " " << n_e_backward << " " << m_number_energies);
75 }
76
77 //..Store generated energies as integers in MeV
78 i_energies.resize(nLeakReg, std::vector<int>(m_number_energies, 0));
79 for (int ie = 0; ie < m_number_energies; ie++) {
80 i_energies[0][ie] = (int)(1000.*m_energies_forward[ie] + 0.5);
81 i_energies[1][ie] = (int)(1000.*m_energies_barrel[ie] + 0.5);
82 i_energies[2][ie] = (int)(1000.*m_energies_backward[ie] + 0.5);
83 }
84
85 //..Require all energies are different, and that there are at least two
86 if (m_number_energies < 2) {
87 B2FATAL("eclLeakageCollector: require at least two energy points. m_number_energies = " << m_number_energies);
88 }
89 for (int ie = 0; ie < m_number_energies - 1; ie++) {
90 for (int ireg = 0; ireg < nLeakReg; ireg++) {
91 if (i_energies[ireg][ie] == i_energies[ireg][ie + 1]) {
92 B2FATAL("eclLeakageCollector: identical energies, ireg = " << ireg << ", " << i_energies[ireg][ie] << " MeV");
93 }
94 }
95 }
96
97 //-----------------------------------------------------------------
98 //..Write out the input parameters
99 B2INFO("eclLeakageCollector parameters: ");
100 B2INFO("position_bins " << m_position_bins);
101 B2INFO("number_energies " << m_number_energies);
102 std::cout << "energies_forward ";
103 for (int ie = 0; ie < m_number_energies; ie++) {std::cout << m_energies_forward[ie] << " ";}
104 std::cout << std::endl;
105 std::cout << "energies_barrel ";
106 for (int ie = 0; ie < m_number_energies; ie++) {std::cout << m_energies_barrel[ie] << " ";}
107 std::cout << std::endl;
108 std::cout << "energies_backward ";
109 for (int ie = 0; ie < m_number_energies; ie++) {std::cout << m_energies_backward[ie] << " ";}
110 std::cout << std::endl;
111 B2INFO("showerArrayName " << m_showerArrayName);
112
113 //-----------------------------------------------------------------
114 //..Define histogram to store parameters
115 const int nBinX = 3 + nLeakReg * m_number_energies;
116 auto inputParameters = new TH1F("inputParameters", "eclLeakageCollector job parameters", nBinX, 0, nBinX);
117 registerObject<TH1F>("inputParameters", inputParameters);
118
119 //..TTree stores required quantities for each photon
120 auto tree = new TTree("single_photon_leakage", "");
121 tree->Branch("cellID", &t_cellID, "cellID/I");
122 tree->Branch("thetaID", &t_thetaID, "thetaID/I");
123 tree->Branch("region", &t_region, "region/I");
124 tree->Branch("thetaBin", &t_thetaBin, "thetaBin/I");
125 tree->Branch("phiBin", &t_phiBin, "phiBin/I");
126 tree->Branch("phiMech", &t_phiMech, "phiMech/I");
127 tree->Branch("energyBin", &t_energyBin, "energyBin/I");
128 tree->Branch("nCrys", &t_nCrys, "nCrys/I");
129 tree->Branch("energyFrac", &t_energyFrac, "energyFrac/F");
130 tree->Branch("origEnergyFrac", &t_origEnergyFrac, "origEnergyFrac/F");
131 tree->Branch("locationError", &t_locationError, "locationError/F");
132 registerObject<TTree>("tree", tree);
133
134
135 //-----------------------------------------------------------------
136 //..Class to find cellID and position within crystal from theta and phi
137 std::cout << "creating leakagePosition object " << std::endl;
139
140 //-----------------------------------------------------------------
141 //..Required arrays
143 m_mcParticleArray.isRequired();
144
145 //-----------------------------------------------------------------
146 //..nOptimal payload. Get optimal number of crystals, and
147 // corresponding correction factors.
148 m_nOptimal2D = m_eclNOptimal->getNOptimal();
149 m_peakFracEnergy = m_eclNOptimal->getPeakFracEnergy();
150 m_bias = m_eclNOptimal->getBias();
151 m_logPeakEnergy = m_eclNOptimal->getLogPeakEnergy();
152 m_groupNumber = m_eclNOptimal->getGroupNumber();
153
154 //..Vectors of energy boundaries for each region
155 std::vector<float> eBoundariesFwd = m_eclNOptimal->getUpperBoundariesFwd();
156 std::vector<float> eBoundariesBrl = m_eclNOptimal->getUpperBoundariesBrl();
157 std::vector<float> eBoundariesBwd = m_eclNOptimal->getUpperBoundariesBwd();
158 m_nEnergyBins = eBoundariesBrl.size();
159
160 //..Copy values to m_eBoundaries
161 m_eBoundaries.resize(m_nLeakReg, std::vector<float>(m_nEnergyBins, 0.));
162 for (int ie = 0; ie < m_nEnergyBins; ie++) {
163 m_eBoundaries[0][ie] = eBoundariesFwd[ie];
164 m_eBoundaries[1][ie] = eBoundariesBrl[ie];
165 m_eBoundaries[2][ie] = eBoundariesBwd[ie];
166 B2INFO(" nOptimal upper boundaries for energy point " << ie << " " << m_eBoundaries[0][ie] << " " << m_eBoundaries[1][ie] << " " <<
167 m_eBoundaries[2][ie]);
168 }
169
170}
171
172//-----------------------------------------------------------------
174{
175
176 //-----------------------------------------------------------------
177 //..First time, store the job parameters
178 if (storeCalib) {
179 getObjectPtr<TH1F>("inputParameters")->Fill(0.01, m_position_bins);
180 getObjectPtr<TH1F>("inputParameters")->Fill(1.01, m_number_energies);
181 double firstBin = 2.01;
182 for (int ie = 0; ie < m_number_energies; ie++) {
183 getObjectPtr<TH1F>("inputParameters")->Fill(firstBin + ie, m_energies_forward[ie]);
184 }
185 firstBin += m_number_energies;
186 for (int ie = 0; ie < m_number_energies; ie++) {
187 getObjectPtr<TH1F>("inputParameters")->Fill(firstBin + ie, m_energies_barrel[ie]);
188 }
189 firstBin += m_number_energies;
190 for (int ie = 0; ie < m_number_energies; ie++) {
191 getObjectPtr<TH1F>("inputParameters")->Fill(firstBin + ie, m_energies_backward[ie]);
192 }
193
194 //..Keep track of how many times inputParameters was filled
195 int lastBin = getObjectPtr<TH1F>("inputParameters")->GetNbinsX();
196 getObjectPtr<TH1F>("inputParameters")->SetBinContent(lastBin, 1.);
197 storeCalib = false;
198 }
199
200 //-----------------------------------------------------------------
201 //..Generated MC particle. Should only be one entry, but check.
202 const int nMC = m_mcParticleArray.getEntries();
203 if (nMC != 1) {return;}
204 double mcLabE = m_mcParticleArray[0]->getEnergy();
205 ROOT::Math::XYZVector mcp3 = m_mcParticleArray[0]->getMomentum();
206 ROOT::Math::XYZVector vertex = m_mcParticleArray[0]->getProductionVertex();
207
208 //-----------------------------------------------------------------
209 //..Find the ECLShower (should only be one when using trimmed data object)
210 const int nShower = m_eclShowerArray.getEntries();
211 double minAngle = 999.;
212 int minShower = -1;
213 for (int is = 0; is < nShower; is++) {
214
215 //..Only interested in photon hypothesis showers
216 if (m_eclShowerArray[is]->getHypothesisId() == ECLShower::c_nPhotons) {
217
218 //..Make a position vector from theta, phi, and R
219 ROOT::Math::XYZVector position;
220 VectorUtil::setMagThetaPhi(
221 position, m_eclShowerArray[is]->getR(),
222 m_eclShowerArray[is]->getTheta(), m_eclShowerArray[is]->getPhi());
223
224 //..Direction is difference between position and vertex
225 ROOT::Math::XYZVector direction = ROOT::Math::XYZVector(position) - vertex;
226
227 double angle = ROOT::Math::VectorUtil::Angle(direction, mcp3);
228 if (angle < minAngle) {
229 minAngle = angle;
230 minShower = is;
231 }
232 }
233 }
234 if (minShower == -1) {return;}
235
236 //-----------------------------------------------------------------
237 //..Location of selected shower.
238 const int maxECellId = m_eclShowerArray[minShower]->getCentralCellId();
239 const float thetaLocation = m_eclShowerArray[minShower]->getTheta();
240 const float phiLocation = m_eclShowerArray[minShower]->getPhi();
241 std::vector<int> positionVector = leakagePosition->getLeakagePosition(maxECellId, thetaLocation, phiLocation, m_position_bins);
242
243 //..TTree items
244 t_cellID = positionVector[0];
245 t_thetaID = positionVector[1];
246 t_region = positionVector[2];
247 t_thetaBin = positionVector[3];
248 t_phiBin = positionVector[4];
249 t_phiMech = positionVector[5];
250
251 //-----------------------------------------------------------------
252 //..Generated and reconstructed energy quantities
253
254 //..Find the generated energy bin by requiring it be close enough to expected value
255 const float genEnergyMeV = 1000.*mcLabE;
256 const float tolerance = std::max(0.002 * genEnergyMeV, 1.0);
257 t_energyBin = -1;
258 for (int ie = 0; ie < m_number_energies; ie++) {
259 if (std::abs(genEnergyMeV - i_energies[t_region][ie]) < tolerance) {
260 t_energyBin = ie;
261 break;
262 }
263 }
264
265 //..Quit if the true energy is not equal to a generated one.
266 // This can happen if the cluster is reconstructed in the wrong region.
267 if (t_energyBin == -1) {return;}
268
269 //..Reconstructed energy after existing leakage correction, normalized to generated
270 t_origEnergyFrac = m_eclShowerArray[minShower]->getEnergy() / mcLabE;
271
272 //..Sum of nOptimal crystals (without leakage correction), when events were generated
273 const double eRaw = m_eclShowerArray[minShower]->getEnergyRaw();
274
275 //-----------------------------------------------------------------
276 //..Find nOptimal from cellID and 3x3 energy found by ECLSplitterN1Module.
277 // Note that payload may have been updated since events were generated, so
278 // we need to redo the sum of energies of the nOptimal crystals.
279 const int ig = m_groupNumber[maxECellId - 1];
280 const double e3x3Orig = m_eclShowerArray[minShower]->getNOptimalEnergy();
281 double e3x3 = e3x3Orig;
282
283 //..Alternate e3x3, for samples produced before nOptimal was introduced
284 const double eHighestCorr = m_eclShowerArray[minShower]->getEnergyHighestCrystal();
285 const double correction = m_eclShowerArray[minShower]->getEnergy() / m_eclShowerArray[minShower]->getEnergyRaw();
286 const double eHighestRaw = eHighestCorr / correction;
287 const double e3x3Alt = eHighestRaw / m_eclShowerArray[minShower]->getE1oE9();
288 if (e3x3 < 0.001) {e3x3 = e3x3Alt;}
289
290 //..Need the detector region to get the energy bin boundaries
291 int iRegion = 1; // barrel
292 if (ECLElementNumbers::isForward(maxECellId)) {iRegion = 0;}
293 if (ECLElementNumbers::isBackward(maxECellId)) {iRegion = 2;}
294
295 //..nOptimal energy bin for this energy.
296 int iEnergy = 0;
297 while (e3x3 > m_eBoundaries[iRegion][iEnergy] and iEnergy < m_nEnergyBins - 1) {iEnergy++;}
298
299 //..nOptimal
300 t_nCrys = m_nOptimal2D.GetBinContent(ig + 1, iEnergy + 1);
301
302 //-----------------------------------------------------------------
303 //..Find the sum of the nOptimal most-energetic digits
304
305 //..Get the ECLCalDigits and rank them by energy
306 std::vector<double> digitEnergies;
307 const auto showerDigitRelations = m_eclShowerArray[minShower]->getRelationsWith<ECLCalDigit>("ECLTrimmedDigits");
308 unsigned int nRelatedDigits = showerDigitRelations.size();
309 for (unsigned int ir = 0; ir < nRelatedDigits; ir++) {
310 const auto calDigit = showerDigitRelations.object(ir);
311 const auto weight = showerDigitRelations.weight(ir);
312 digitEnergies.push_back(calDigit->getEnergy() * weight);
313 }
314
315 //..Rank from lowest to highest
316 std::sort(digitEnergies.begin(), digitEnergies.end());
317
318 //..Sum the highest nOptimal digit energies (if there are that many)
319 double eSumOfN = 1.e-10; // cpp does not like 0
320 for (int isum = 0; isum < t_nCrys; isum++) {
321 const int i = (int)nRelatedDigits - 1 - isum;
322 if (i >= 0) {eSumOfN += digitEnergies[i];}
323 }
324
325 //-----------------------------------------------------------------
326 //..Correct the sum of nOptimal crystals for bias and nCrystal peak energy
327 // We need to do this because the nOptimal payload may have changed
328 // since the events were generated.
329
330 //..To allow for energy interpolation, there are three bins per group and energy:
331 // iy = iEnergy + 1 in all three cases
332 // ix = 3*ig + 1 = value for nominal energy, i.e logPeakEnergy(3*ig+1, iEnergy+1)
333 // ix = 3*ig + 2 = value for lower energy, i.e logPeakEnergy(3*ig+2, iEnergy+1)
334 // ix = 3*ig + 3 = value for higher energy, i.e logPeakEnergy(3*ig+3, iEnergy+1)
335 const int iy = iEnergy + 1;
336
337 const int ixNom = 3 * ig + 1;
338 const int ixLowerE = ixNom + 1;
339 const int ixHigherE = ixNom + 2;
340
341 const double logENom = m_logPeakEnergy.GetBinContent(ixNom, iy);
342 const double logELower = m_logPeakEnergy.GetBinContent(ixLowerE, iy);
343 const double logEHigher = m_logPeakEnergy.GetBinContent(ixHigherE, iy);
344
345 const double biasNom = m_bias.GetBinContent(ixNom, iy);
346 const double biasLower = m_bias.GetBinContent(ixLowerE, iy);
347 const double biasHigher = m_bias.GetBinContent(ixHigherE, iy);
348
349 const double peakNom = m_peakFracEnergy.GetBinContent(ixNom, iy);
350 const double peakLower = m_peakFracEnergy.GetBinContent(ixLowerE, iy);
351 const double peakHigher = m_peakFracEnergy.GetBinContent(ixHigherE, iy);
352
353 //..Interpolate in log energy
354 const double logESumN = log(eSumOfN);
355
356 double logEOther = logELower;
357 double biasOther = biasLower;
358 double peakOther = peakLower;
359 if (logESumN > logENom) {
360 logEOther = logEHigher;
361 biasOther = biasHigher;
362 peakOther = peakHigher;
363 }
364
365 //..The nominal and "other" energies may be identical if this is the first or last energy
366 double bias = biasNom;
367 double peak = peakNom;
368 if (std::abs(logEOther - logENom) > 0.0001) {
369 bias = biasNom + (biasOther - biasNom) * (logESumN - logENom) / (logEOther - logENom);
370 peak = peakNom + (peakOther - peakNom) * (logESumN - logENom) / (logEOther - logENom);
371 }
372
373 //..Normalized reconstructed energy after bias and nCrystal correction
374 t_energyFrac = (eSumOfN - bias) / peak / mcLabE;
375
376 //-----------------------------------------------------------------
377 //..Distance between generated and reconstructed positions
378 const double radius = m_eclShowerArray[minShower]->getR();
379 ROOT::Math::XYZVector measuredLocation;
380 VectorUtil::setMagThetaPhi(
381 measuredLocation, radius, thetaLocation, phiLocation);
382 ROOT::Math::XYZVector measuredDirection = ROOT::Math::XYZVector(measuredLocation) - vertex;
383 t_locationError = radius * ROOT::Math::VectorUtil::Angle(measuredDirection, mcp3);
384
385
386 //-----------------------------------------------------------------
387 //..Debug: dump some events
388 if (m_nDump < 100) {
389 m_nDump++;
390 B2INFO(" ");
391 B2INFO("cellID " << t_cellID << " ig " << ig << " iEnergy " << iEnergy << " ie " << t_energyBin << " nOpt " << t_nCrys);
392 B2INFO(" e3x3Orig " << e3x3Orig << " e3x3Alt " << e3x3Alt << " Eraw " << eRaw << " ESum " << eSumOfN << " eFrac " <<
394 B2INFO(" 3 log E " << logELower << " " << logENom << " " << logEHigher << " log E " << logESumN);
395 B2INFO(" 3 biases " << biasLower << " " << biasNom << " " << biasHigher << " bias " << bias);
396 B2INFO(" 3 peaks " << peakLower << " " << peakNom << " " << peakHigher << " peak " << peak);
397 }
398
399 //-----------------------------------------------------------------
400 //..Done
401 getObjectPtr<TTree>("tree")->Fill();
402}
403
void registerObject(std::string name, T *obj)
Register object with a name, takes ownership, do not access the pointer beyond prepare()
CalibrationCollectorModule()
Constructor. Sets the default prefix for calibration dataobjects.
T * getObjectPtr(std::string name)
Calls the CalibObjManager to get the requested stored collector data.
Class to store calibrated ECLDigits: ECLCalDigits.
Definition ECLCalDigit.h:23
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLShower.h:42
Class to get position information for a cluster for leakage corrections.
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
StoreArray< MCParticle > m_mcParticleArray
Required input array of MCParticles.
eclLeakageCollectorModule()
Constructor: Sets the description, the properties and the parameters of the module.
StoreArray< ECLShower > m_eclShowerArray
Required input array of ECLShowers.
TH2F m_logPeakEnergy
log of peak contained energy in GeV
float t_energyFrac
measured energy without leakage correction divided by generated
const int m_nLeakReg
3 ECL regions: 0 = forward, 1 = barrel, 2 = backward
float t_locationError
reconstructed minus generated position (cm)
std::vector< double > m_energies_forward
generated photon energies, forward
const int nLeakReg
Some other useful quantities.
int m_number_energies
number of generated energies (8)
std::vector< std::vector< int > > i_energies
Generated energies in MeV in each region.
bool storeCalib
store parameters first event
TH2F m_bias
2D histogram of bias = sum of ECLCalDigit energy minus true (GeV)
TH2F m_peakFracEnergy
2D histogram of peak fractional contained energy
virtual void collect() override
Accumulate TTree.
std::string m_showerArrayName
Required arrays.
StoreObjPtr< EventMetaData > m_evtMetaData
dataStore EventMetaData
DBObjPtr< ECLnOptimal > m_eclNOptimal
nOptimal payload
std::vector< double > m_energies_backward
generated photon energies, backward
ECL::ECLLeakagePosition * leakagePosition
location of position of cluster
virtual void prepare() override
Define histograms and read payloads from DB.
int t_phiBin
binned location in phi relative to crystal edge.
int t_phiMech
0: mechanical structure next to phi edge; 1: no mech structure
int t_region
region of photon 0=forward 1=barrel 2=backward
int m_nEnergyBins
number of energies bins in nOptimal payload
int t_thetaBin
binned location in theta relative to crystal edge
int t_nCrys
number of crystals used to calculate energy
std::vector< double > m_energies_barrel
generated photon energies, barrel
int m_position_bins
Parameters to pass along to the algorithm.
float t_origEnergyFrac
original leakage-corrected energy / generated
std::vector< std::vector< float > > m_eBoundaries
energy boundaries each region
TH2F m_nOptimal2D
2D hist of nOptimal for Ebin vs groupID
std::vector< int > m_groupNumber
group number for each crystal
int m_nDump
Number of events with diagnostic printouts.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
bool isForward(int cellId)
Check whether the crystal is in forward ECL.
bool isBackward(int cellId)
Check whether the crystal is in backward ECL.
Abstract base class for different kinds of events.