11 #include <ecl/modules/eclDigitizer/ECLDigitizerPureCsIModule.h>
20 #include <framework/gearbox/Unit.h>
21 #include <framework/logging/Logger.h>
22 #include <framework/utilities/FileSystem.h>
25 #include <ecl/digitization/ECLDspFitterPure.h>
26 #include <ecl/dataobjects/ECLHit.h>
27 #include <ecl/dataobjects/ECLDigit.h>
28 #include <ecl/dataobjects/ECLDsp.h>
29 #include <ecl/dataobjects/ECLPureCsIInfo.h>
30 #include <ecl/geometry/ECLGeometryPar.h>
31 #include <ecl/dbobjects/ECLWaveformData.h>
49 setDescription(
"Creates ECLDigiHits from ECLHits for Pure CsI.");
50 setPropertyFlags(c_ParallelProcessingCertified);
51 addParam(
"FirstRing", m_thetaIdMin,
"First ring (0-12)", 0);
52 addParam(
"LastRing", m_thetaIdMax,
"Last ring (0-12)", 12);
53 addParam(
"Background", m_background,
"Flag to use the DigitizerPureCsI configuration with backgrounds; Default is no background",
55 addParam(
"Calibration", m_calibration,
56 "Flag to use the DigitizerPureCsI for Waveform fit Covariance Matrix calibration; Default is false",
58 addParam(
"adcTickFactor", m_tickFactor,
"multiplication factor to get adc tick from trigger tick", 8);
59 addParam(
"sigmaTrigger", m_sigmaTrigger,
"Trigger resolution used", 0.);
60 addParam(
"elecNoise", m_elecNoise,
"Electronics noise energy equivalent in MeV", 0.5);
64 addParam(
"photostatresolution", m_photostatresolution,
"sigma for 1 MeV energy deposit", 0.22);
65 addParam(
"Debug", m_debug,
"debug mode on (default off)",
false);
66 addParam(
"debugtrgtime", m_testtrg,
"set fixed trigger time for testing purposes", 0);
67 addParam(
"debugsigtimeshift", m_testsig,
"shift signal arrival time for testing purposes (in microsec)", 0.);
68 addParam(
"debugenergydeposit", m_testenedep,
"energy deposit in all crystals for testing purposes", 0.);
69 addParam(
"NoCovMatrix", m_NoCovMatrix,
"Use a diagonal (neutral) Covariance matrix",
true);
73 ECLDigitizerPureCsIModule::~ECLDigitizerPureCsIModule()
77 void ECLDigitizerPureCsIModule::initialize()
80 EclConfigurationPure::m_tickPure = m_tickFactor * EclConfiguration::m_tick / EclConfiguration::m_ntrg;
82 m_ecldsps.registerInDataStore(eclDspArrayName());
84 m_ecldigits.registerInDataStore(eclDigitArrayName());
86 m_ecldigits.registerRelationTo(m_ecldsps);
89 m_eclpurecsiinfo.registerInDataStore(eclPureCsIInfoArrayName());
90 m_ecldigits.registerRelationTo(m_eclpurecsiinfo);
93 m_ecldigits.registerRelationTo(m_hitLists);
96 m_adc.resize(EclConfigurationPure::m_nch);
101 void ECLDigitizerPureCsIModule::beginRun()
105 void ECLDigitizerPureCsIModule::event()
112 double deltaT = m_sigmaTrigger == 0 ? 0 : gRandom->Gaus(0, m_sigmaTrigger);
118 vector< vector< const ECLHit*> > hitmap(EclConfigurationPure::m_nch);
120 for (
const auto& eclHit : eclHits) {
121 int j = eclHit.getCellId() - 1;
122 if (isPureCsI(j + 1)) {
123 assert(j < EclConfigurationPure::m_nch);
124 double hitE = eclHit.getEnergyDep() / Unit::GeV;
125 double hitTime = eclHit.getTimeAve() / Unit::us;
126 if (m_photostatresolution > 0) {
127 double nphotavg1MeV = 1 / (m_photostatresolution * m_photostatresolution);
128 int nphotavg = round((hitE / 0.001) * nphotavg1MeV);
129 int nphot = gRandom->Poisson(nphotavg);
130 hitE = (nphot / nphotavg1MeV) / 1000;
134 m_adc[j].AddHit(hitE , hitTime + deltaT, m_ss[m_tbl[j].iss]);
135 if (eclHit.getBackgroundTag() == BackgroundMetaData::bg_none) hitmap[j].push_back(&eclHit);
141 for (
int j = 0; j < EclConfigurationPure::m_nch; j++) {
142 if (! isPureCsI(j + 1))
continue;
145 m_adc[j].AddHit(m_testenedep, m_testsig, m_ss[m_tbl[j].iss]);
149 if (! m_calibration && a.total < 0.0001)
continue;
152 float adcNoise[EclConfigurationPure::m_nsmp];
153 memset(adcNoise, 0,
sizeof(adcNoise));
154 if (m_elecNoise > 0) {
155 float z[EclConfigurationPure::m_nsmp];
156 for (
int i = 0; i < EclConfigurationPure::m_nsmp; i++)
157 z[i] = gRandom->Gaus(0, 1);
158 m_noise[0].generateCorrelatedNoise(z, adcNoise);
161 int FitA[EclConfigurationPure::m_nsmp];
162 for (
int i = 0; i < EclConfigurationPure::m_nsmp; i++) {
163 FitA[i] = 20 * (1000 * a.c[i] + adcNoise[i]) + 3000;
170 if (! m_calibration) {
173 DSPFitterPure(m_fitparams[m_tbl[j].idn], FitA, m_testtrg, energyFit, tFit, fitChi2, qualityFit);
181 DSPFitterPure(m_fitparams[m_tbl[j].idn], FitA, 0, energyFit, tFit, fitChi2, qualityFit);
191 if (m_calibration || energyFit > 0) {
193 auto eclDsp = m_ecldsps.appendNew();
194 eclDsp->setCellId(CellId);
195 eclDsp->setDspA(FitA);
197 auto eclDigit = m_ecldigits.appendNew();
198 eclDigit->setCellId(CellId);
199 eclDigit->setAmp(energyFit);
200 eclDigit->setTimeFit(
int(tFit * 10));
201 eclDigit->setQuality(qualityFit);
203 auto AeclPureCsIInfo = m_eclpurecsiinfo.appendNew();
204 eclDigit->addRelationTo(AeclPureCsIInfo);
205 AeclPureCsIInfo->setPureCsI(1);
206 AeclPureCsIInfo->setCellId(CellId);
208 eclDigit->addRelationTo(eclDsp);
209 for (
const auto& hit : hitmap[j])
210 eclDigit->addRelationTo(hit);
219 for (
const auto& eclDigit : baselineDigits) {
220 int cellid = eclDigit.getCellId();
221 if (! isPureCsI(cellid)) {
222 auto eclDigitClone = m_ecldigits.appendNew();
223 eclDigitClone->setCellId(cellid);
224 eclDigitClone->setAmp(eclDigit.getAmp());
225 eclDigitClone->setTimeFit(eclDigit.getTimeFit());
226 eclDigitClone->setQuality(eclDigit.getQuality());
228 auto AeclPureCsIInfo = m_eclpurecsiinfo.appendNew();
229 eclDigitClone->addRelationTo(AeclPureCsIInfo);
230 AeclPureCsIInfo->setPureCsI(0);
231 AeclPureCsIInfo->setCellId(cellid);
239 void ECLDigitizerPureCsIModule::endRun()
243 void ECLDigitizerPureCsIModule::terminate()
247 void ECLDigitizerPureCsIModule::readDSPDB()
249 string dataFileName, dataFileName2;
251 dataFileName = FileSystem::findFile(
"/data/ecl/ECL-WF-Pure.root");
253 dataFileName2 = FileSystem::findFile(
"/data/ecl/ECL-WF-cov-Pure-BG.root");
254 B2INFO(
"ECLDigitizerPureCsI: Reading configuration data with background from: " << dataFileName);
255 B2INFO(
"ECLDigitizerPureCsI: Reading configuration data with background from: " << dataFileName2);
258 dataFileName = FileSystem::findFile(
"/data/ecl/ECL-WF-Pure.root");
260 dataFileName2 = FileSystem::findFile(
"/data/ecl/ECL-WF-cov-Pure.root");
261 B2INFO(
"ECLDigitizerPureCsI: Reading configuration data without background from: " << dataFileName);
262 B2INFO(
"ECLDigitizerPureCsI: Reading configuration data without background from: " << dataFileName2);
264 assert(! dataFileName.empty());
266 TFile rootfile(dataFileName.c_str());
267 const TH1F* sampledWF =
dynamic_cast<TH1F*
>(rootfile.Get(
"sampleddsp"));
268 assert(sampledWF !=
nullptr);
269 const TH1F* sampledWF1 =
dynamic_cast<TH1F*
>(rootfile.Get(
"sampleddsp1"));
270 assert(sampledWF1 !=
nullptr);
272 m_tbl.resize(EclConfigurationPure::m_nch);
277 m_ss[0].InitSample(sampledWF, sampledWF1);
279 for (
int i = 0; i < EclConfigurationPure::m_nch; i++) m_tbl[i].iss = 0;
280 B2INFO(
"ECLDigitizerPureCsI: " << m_ss.size() <<
" sampled signal templates were created.");
284 if (!(m_calibration || m_NoCovMatrix)) {
285 TFile rootfile2(dataFileName2.c_str());
286 TTree* tree = (TTree*) rootfile2.Get(
"EclWF");
288 const int maxncellid = 512;
290 vector<int> cellId(maxncellid);
292 tree->SetBranchAddress(
"ncellId", &ncellId);
293 tree->SetBranchAddress(
"cellId", cellId.data());
294 tree->SetBranchAddress(
"CovarianceM", &eclWFData);
295 for (Long64_t j = 0, jmax = tree->GetEntries(); j < jmax; j++) {
297 assert(ncellId <= maxncellid);
298 for (
int i = 0; i < ncellId; ++i)
299 m_tbl[cellId[i]].idn = m_fitparams.size();
302 m_fitparams.push_back(params);
305 B2INFO(
"ECLDigitizerPureCsI: parameters vector size : " << m_fitparams.size());
308 m_fitparams.resize(1);
309 for (
int i = 0; i < EclConfigurationPure::m_nch; i++)
311 for (
int i = 0; i < 16; i++)
312 for (
int j = 0; j < 16; j++)
313 if (i != j) m_fitparams[0].invC[i][j] = 0;
314 else m_fitparams[0].invC[i][j] = 1.0;
315 initParams(m_fitparams[0], m_ss[0]);
317 for (
auto& param : m_fitparams) {
318 initParams(param, m_ss[0]);
325 for (
int i = 0; i < EclConfigurationPure::m_nsmp; i++)
326 for (
int j = 0; j <= i; j++)
327 if (i == j) m_noise[0].setMatrixElement(index++, m_elecNoise);
328 else m_noise[0].setMatrixElement(index++, 0.);
331 m_noise[0].getMatrix(testM);
334 void ECLDigitizerPureCsIModule::mapGeometry()
337 for (
int cellId0 = 0; cellId0 < EclConfigurationPure::m_nch; cellId0++) {