Belle II Software  release-05-02-19
ECLDigitizerPureCsIModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Guglielmo De Nardo *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 //This module
11 #include <ecl/modules/eclDigitizer/ECLDigitizerPureCsIModule.h>
12 
13 // ROOT
14 #include <TRandom.h>
15 #include <TFile.h>
16 #include <TTree.h>
17 #include <TH1.h>
18 
19 //Framework
20 #include <framework/gearbox/Unit.h>
21 #include <framework/logging/Logger.h>
22 #include <framework/utilities/FileSystem.h>
23 
24 //ECL
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>
32 
33 using namespace std;
34 using namespace Belle2;
35 using namespace ECL;
36 
37 //-----------------------------------------------------------------
38 // Register the Module
39 //-----------------------------------------------------------------
40 REG_MODULE(ECLDigitizerPureCsI)
41 
42 //-----------------------------------------------------------------
43 // Implementation
44 //-----------------------------------------------------------------
45 
47 {
48  //Set module properties
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",
54  false);
55  addParam("Calibration", m_calibration,
56  "Flag to use the DigitizerPureCsI for Waveform fit Covariance Matrix calibration; Default is false",
57  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);
61  /* photo statistics resolution measurement at LNF sigma = 55 % at 1 MeV
62  Csi(Tl) is 12%
63  */
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);
70 
71 }
72 
73 ECLDigitizerPureCsIModule::~ECLDigitizerPureCsIModule()
74 {
75 }
76 
77 void ECLDigitizerPureCsIModule::initialize()
78 {
79  m_nEvent = 0 ;
80  EclConfigurationPure::m_tickPure = m_tickFactor * EclConfiguration::m_tick / EclConfiguration::m_ntrg;
81 
82  m_ecldsps.registerInDataStore(eclDspArrayName());
83 
84  m_ecldigits.registerInDataStore(eclDigitArrayName());
85 
86  m_ecldigits.registerRelationTo(m_ecldsps);
87 
88 
89  m_eclpurecsiinfo.registerInDataStore(eclPureCsIInfoArrayName());
90  m_ecldigits.registerRelationTo(m_eclpurecsiinfo);
91 
92 
93  m_ecldigits.registerRelationTo(m_hitLists);
94  readDSPDB();
95 
96  m_adc.resize(EclConfigurationPure::m_nch);
97 
98  mapGeometry();
99 }
100 
101 void ECLDigitizerPureCsIModule::beginRun()
102 {
103 }
104 
105 void ECLDigitizerPureCsIModule::event()
106 {
107  StoreArray<ECLHit> eclHits;
108 
109  /* add trigger resolution defined in a module paramer
110  shifting the waveform starting time by a random deltaT,
111  assuming that t0=0 adc channel is determined by the trigger */
112  double deltaT = m_sigmaTrigger == 0 ? 0 : gRandom->Gaus(0, m_sigmaTrigger);
113 
114  // clear the storage for the event
115  memset(m_adc.data(), 0, m_adc.size()*sizeof(adccounts_type));
116 
117  // emulate response for ECL hits after ADC measurements
118  vector< vector< const ECLHit*> > hitmap(EclConfigurationPure::m_nch);
119 
120  for (const auto& eclHit : eclHits) {
121  int j = eclHit.getCellId() - 1; //0~8735
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;
131 
132  // hitE = gRandom->Gaus(hitE, 0.001 * m_photostatresolution * sqrt(hitE * 1000));
133  }
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);
136  }
137  }
138 
139  // loop over entire calorimeter
140 
141  for (int j = 0; j < EclConfigurationPure::m_nch; j++) {
142  if (! isPureCsI(j + 1)) continue;
143 
144  if (m_debug) {
145  m_adc[j].AddHit(m_testenedep, m_testsig, m_ss[m_tbl[j].iss]);
146  //cout << "Adding enedep = " << m_testenedep << " time: " << m_testsig << endl;
147  }
148  adccounts_type& a = m_adc[j];
149  if (! m_calibration && a.total < 0.0001) continue;
150 
151  //Noise generation
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);
159  }
160 
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;
164  }
165 
166  int energyFit = 0; // fit output : Amplitude 18 bits
167  double tFit = 0; // fit output : T_ave 12 bits
168  int qualityFit = 0; // fit output : quality 2 bits
169 
170  if (! m_calibration) {
171  double fitChi2 = 0;
172  if (m_debug) {
173  DSPFitterPure(m_fitparams[m_tbl[j].idn], FitA, m_testtrg, energyFit, tFit, fitChi2, qualityFit);
174  /*
175  cout << "energy: " << energyFit
176  << " tFit: " << tFit
177  << " qualityfit: " << qualityFit
178  << endl;
179  */
180  } else {
181  DSPFitterPure(m_fitparams[m_tbl[j].idn], FitA, 0, energyFit, tFit, fitChi2, qualityFit);
182  /*
183  cout << "energy: " << energyFit
184  << " tFit: " << tFit
185  << " qualityfit: " << qualityFit
186  << endl;
187  */
188  }
189  }
190 
191  if (m_calibration || energyFit > 0) {
192  int CellId = j + 1;
193  auto eclDsp = m_ecldsps.appendNew();
194  eclDsp->setCellId(CellId);
195  eclDsp->setDspA(FitA);
196 
197  auto eclDigit = m_ecldigits.appendNew();
198  eclDigit->setCellId(CellId); // cellId in range from 1 to 8736
199  eclDigit->setAmp(energyFit); // E (GeV) = energyFit/20000;
200  eclDigit->setTimeFit(int(tFit * 10)); // time is in 0.1 ns units
201  eclDigit->setQuality(qualityFit);
202 
203  auto AeclPureCsIInfo = m_eclpurecsiinfo.appendNew();
204  eclDigit->addRelationTo(AeclPureCsIInfo);
205  AeclPureCsIInfo->setPureCsI(1);
206  AeclPureCsIInfo->setCellId(CellId);
207 
208  eclDigit->addRelationTo(eclDsp);
209  for (const auto& hit : hitmap[j])
210  eclDigit->addRelationTo(hit);
211  }
212  } //store each crystal hit
213 
214  // temporary solution to run Pure CsI reconstruction
215  // and baseline independently and simultaneously
216  // cloning barrel and bwd digits
217 
218  StoreArray<ECLDigit> baselineDigits("ECLDigits");
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());
227  //eclDigitClone->setPureCsI(0);
228  auto AeclPureCsIInfo = m_eclpurecsiinfo.appendNew();
229  eclDigitClone->addRelationTo(AeclPureCsIInfo);
230  AeclPureCsIInfo->setPureCsI(0);
231  AeclPureCsIInfo->setCellId(cellid);
232  }
233  }
234 
235 
236  m_nEvent++;
237 }
238 
239 void ECLDigitizerPureCsIModule::endRun()
240 {
241 }
242 
243 void ECLDigitizerPureCsIModule::terminate()
244 {
245 }
246 
247 void ECLDigitizerPureCsIModule::readDSPDB()
248 {
249  string dataFileName, dataFileName2;
250  if (m_background) {
251  dataFileName = FileSystem::findFile("/data/ecl/ECL-WF-Pure.root");
252  if (! m_calibration)
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);
256 
257  } else {
258  dataFileName = FileSystem::findFile("/data/ecl/ECL-WF-Pure.root");
259  if (! m_calibration)
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);
263  }
264  assert(! dataFileName.empty());
265 
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);
271 
272  m_tbl.resize(EclConfigurationPure::m_nch);
273 
274  // at the moment there is only one sampled signal shape in the pool
275  // since all shaper parameters are the same for all crystals
276  m_ss.resize(1);
277  m_ss[0].InitSample(sampledWF, sampledWF1);
278 
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.");
281 
282  rootfile.Close();
283 
284  if (!(m_calibration || m_NoCovMatrix)) {
285  TFile rootfile2(dataFileName2.c_str());
286  TTree* tree = (TTree*) rootfile2.Get("EclWF");
287  ECLWaveformData* eclWFData = new ECLWaveformData;
288  const int maxncellid = 512;
289  int ncellId;
290  vector<int> cellId(maxncellid);//[ncellId] buffer for crystal identification number
291 
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++) {
296  tree->GetEntry(j);
297  assert(ncellId <= maxncellid);
298  for (int i = 0; i < ncellId; ++i)
299  m_tbl[cellId[i]].idn = m_fitparams.size();
300  fitparams_type params;
301  eclWFData->getMatrix(params.invC);
302  m_fitparams.push_back(params);
303  }
304  }
305  B2INFO("ECLDigitizerPureCsI: parameters vector size : " << m_fitparams.size());
306  // at the moment there is one set of fitparams
307  if (m_NoCovMatrix) {
308  m_fitparams.resize(1);
309  for (int i = 0; i < EclConfigurationPure::m_nch; i++)
310  m_tbl[i].idn = 0;
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]);
316  } else {
317  for (auto& param : m_fitparams) {
318  initParams(param, m_ss[0]);
319  }
320  }
321 
322  // at the moment same noise for all crystals
323  m_noise.resize(1);
324  int index = 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); // units are MeV energy noise eq from electronics
328  else m_noise[0].setMatrixElement(index++, 0.); //uncorrelated
329 
330  float testM[31][31];
331  m_noise[0].getMatrix(testM);
332 }
333 
334 void ECLDigitizerPureCsIModule::mapGeometry()
335 {
336  ECLGeometryPar* eclgeo = ECLGeometryPar::Instance();
337  for (int cellId0 = 0; cellId0 < EclConfigurationPure::m_nch; cellId0++) {
338  eclgeo->Mapping(cellId0);
339  m_thetaID[cellId0] = eclgeo->GetThetaID();
340  }
341 }
Belle2::ECL::EclConfigurationPure::fitparamspure_t
a struct for the fit parameters for the pure CsI calorimeter
Definition: EclConfigurationPure.h:57
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECL::EclConfigurationPure::adccountspure_t
a struct for the fit parameters for the pure CsI calorimeter
Definition: EclConfigurationPure.h:50
Belle2::ECLDigitizerPureCsIModule
The ECLDigitizerPureCsI module.
Definition: ECLDigitizerPureCsIModule.h:57
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::ECLGeometryPar::Mapping
void Mapping(int cid)
Mapping theta, phi Id.
Definition: ECLGeometryPar.cc:387
Belle2::ECLWaveformData
ECLWaveformData - container for inverse covariant matrix and shape parameters for time and amplitude ...
Definition: ECLWaveformData.h:39
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition: ECLGeometryPar.h:45
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::ECLWaveformData::getMatrix
void getMatrix(float M[16][16]) const
Getter method for all matrix as two dimentional array (floats)
Definition: ECLWaveformData.h:66
Belle2::ECL::ECLGeometryPar::GetThetaID
int GetThetaID()
Get Theta Id.
Definition: ECLGeometryPar.h:78