Belle II Software  release-05-01-25
ARICHDigitizerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Luka Santelj *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 // Own include
11 #include <arich/modules/arichDigitizer/ARICHDigitizerModule.h>
12 
13 // Hit classes
14 #include <arich/dataobjects/ARICHSimHit.h>
15 #include <arich/dataobjects/ARICHDigit.h>
16 
17 // framework - DataStore
18 #include <framework/datastore/DataStore.h>
19 #include <framework/datastore/StoreArray.h>
20 
21 // framework aux
22 #include <framework/logging/Logger.h>
23 
24 // magnetic field manager
25 #include <framework/geometry/BFieldManager.h>
26 
27 #include <framework/dataobjects/BackgroundInfo.h>
28 
29 
30 // ROOT
31 #include <TVector2.h>
32 #include <TVector3.h>
33 #include <TRandom3.h>
34 
35 using namespace std;
36 using namespace boost;
37 
38 namespace Belle2 {
44  //-----------------------------------------------------------------
45  // Register the Module
46  //-----------------------------------------------------------------
47 
48  REG_MODULE(ARICHDigitizer)
49 
50 
51  //-----------------------------------------------------------------
52  // Implementation
53  //-----------------------------------------------------------------
54 
56  m_maxQE(0)
57  {
58 
59  // Set description()
60  setDescription("This module creates ARICHDigits from ARICHSimHits. Here spatial digitization is done, channel-by-channel QE is applied, and readout time window cut is applied.");
61 
62  // Set property flags
63  setPropertyFlags(c_ParallelProcessingCertified);
64 
65  // Add parameters
66  addParam("TimeWindow", m_timeWindow, "Readout time window width in ns", 250.0);
67  addParam("BackgroundHits", m_bkgLevel, "Number of background hits per hapd per readout (electronics noise)", 0.066);
68  addParam("MagneticFieldDistorsion", m_bdistort, "Apply image distortion due to non-perp. magnetic field", 0);
69  }
70 
71  ARICHDigitizerModule::~ARICHDigitizerModule()
72  {
73 
74  }
75 
76  void ARICHDigitizerModule::initialize()
77  {
78  // Print set parameters
79  printModuleParams();
80 
81  // QE at 400nm (3.1eV) applied in SensitiveDetector
82  m_maxQE = m_simPar->getQE(3.1);
83 
85  digits.registerInDataStore();
86 
87  m_bgOverlay = false;
88  StoreObjPtr<BackgroundInfo> bgInfo("", DataStore::c_Persistent);
89  if (bgInfo.isValid()) {
90  if (bgInfo->getMethod() == BackgroundInfo::c_Overlay) m_bgOverlay = true;
91  }
92 
93  }
94 
95  void ARICHDigitizerModule::beginRun()
96  {
97  if (m_simPar->getNBkgHits() > 0) m_bkgLevel = m_simPar->getNBkgHits();
98  }
99 
100  void ARICHDigitizerModule::event()
101  {
102 
103  // Get the collection of ARICHSimHits from the Data store.
104  //------------------------------------------------------
105  StoreArray<ARICHSimHit> arichSimHits;
106  //-----------------------------------------------------
107 
108  // Get the collection of arichDigits from the Data store,
109  // (or have one created)
110  //-----------------------------------------------------
111  StoreArray<ARICHDigit> arichDigits;
112 
113  //---------------------------------------------------------------------
114  // Convert SimHits one by one to digitizer hits.
115  //---------------------------------------------------------------------
116 
117  // We try to include the effect of opposite-polarity crosstalk among channels
118  // on each chip, which depend on number of p.e. on chip
119 
120  std::map<pair<int, int>, int> photoElectrons; // this contains number of photoelectrons falling on each channel
121  std::map<pair<int, int>, int> chipHits; // this contains number of photoelectrons on each chip
122 
123  // Get number of photon hits in this event
124  int nHits = arichSimHits.getEntries();
125 
126  // Loop over all photon hits
127  for (int iHit = 0; iHit < nHits; ++iHit) {
128 
129  ARICHSimHit* aSimHit = arichSimHits[iHit];
130 
131  // check for time window
132  double globaltime = aSimHit->getGlobalTime();
133 
134  if (globaltime < 0. || globaltime > m_timeWindow) continue;
135 
136  TVector2 locpos(aSimHit->getLocalPosition().X(), aSimHit->getLocalPosition().Y());
137 
138  // Get id of module
139  int moduleID = aSimHit->getModuleID();
140 
141  if (m_bdistort) magFieldDistorsion(locpos, moduleID);
142 
143  // skip if not active
144  if (!m_modInfo->isActive(moduleID)) continue;
145 
146  // Get id number of hit channel
147  int chX, chY;
148  m_geoPar->getHAPDGeometry().getXYChannel(locpos.X(), locpos.Y(), chX, chY);
149 
150  if (chX < 0 && chY < 0) continue;
151 
152  int asicChannel = m_chnMap->getAsicFromXY(chX, chY);
153 
154 
155  // eliminate un-active channels
156  if (asicChannel < 0 || !m_chnMask->isActive(moduleID, asicChannel)) continue;
157 
158  // apply channel dependent QE scale factor
159  //double qe_scale = 0.27 / m_maxQE;
160  //double qe_scale = m_modInfo->getChannelQE(moduleID, asicChannel) * m_simPar->getColEff() / m_maxQE; // eventually move collection efficiency to here!
161  double qe_scale = m_modInfo->getChannelQE(moduleID, asicChannel) / m_maxQE;
162 
163  if (qe_scale > 1.) B2ERROR("Channel QE is higher than QE applied in SensitiveDetector");
164  if (gRandom->Uniform(1.) > qe_scale) continue;
165 
166  // photon was converted to photoelectron
167  chipHits[make_pair(moduleID, asicChannel / 36)] += 1;
168  photoElectrons[make_pair(moduleID, asicChannel)] += 1;
169 
170  }
171 
172  // loop over produced photoelectrons. Apply suppression due to the reverse polarization crosstalk
173  // among channels on the same chip, and produce hit bitmap (4 bits).
174 
175  for (std::map<std::pair<int, int> , int>::iterator it = photoElectrons.begin(); it != photoElectrons.end(); ++it) {
176 
177  std::pair<int, int> modch = it->first;
178  double npe = double(it->second);
179 
180  // reduce efficiency
181  npe /= (1.0 + m_simPar->getChipNegativeCrosstalk() * (double(chipHits[make_pair(modch.first, modch.second / 36)]) - 1.0));
182  if (npe < 1.0 && gRandom->Uniform(1) > npe) continue;
183 
184  // Make hit bitmap (depends on number of p.e. on channel). For now bitmap is 0001 for single p.e., 0011 for 2 p.e., ...
185  // More proper implementation is to be done ...
186  uint8_t bitmap = 0;
187  for (int i = 0; i < npe; i++) {
188  bitmap |= 1 << i;
189  if (i == 3) break;
190  }
191 
192  // make new digit!
193  arichDigits.appendNew(modch.first, modch.second, bitmap);
194 
195  }
196 
197  //--- if background not overlayed add electronic noise hits
198  if (m_bgOverlay) return;
199  uint8_t bitmap = 1;
200  unsigned nSlots = m_geoPar->getDetectorPlane().getNSlots();
201  for (unsigned id = 1; id < nSlots + 1; id++) {
202  if (!m_modInfo->isActive(id)) continue;
203  int nbkg = gRandom->Poisson(m_bkgLevel);
204  for (int i = 0; i < nbkg; i++) {
205  arichDigits.appendNew(id, gRandom->Integer(144), bitmap);
206  }
207  }
208 
209  }
210 
211  void ARICHDigitizerModule::magFieldDistorsion(TVector2& hit, int copyno)
212  {
213 
214  TVector2 hitGlob;
215  double phi = m_geoPar->getDetectorPlane().getSlotPhi(copyno);
216  double r = m_geoPar->getDetectorPlane().getSlotR(copyno);
217  double z = m_geoPar->getDetectorZPosition() + m_geoPar->getHAPDGeometry().getWinThickness();
218  hitGlob.SetMagPhi(r, phi);
219  TVector2 shift = hit;
220  hitGlob += shift.Rotate(phi);
221  TVector3 Bfield = BFieldManager::getField(m_geoPar->getMasterVolume().pointToGlobal(TVector3(hitGlob.X(), hitGlob.Y(), z)));
222  double cc = m_geoPar->getHAPDGeometry().getPhotocathodeApdDistance() / abs(Bfield.Z());
223  shift.SetX(cc * Bfield.X());
224  shift.SetY(cc * Bfield.Y());
225  shift = shift.Rotate(-phi);
226  hit.SetX(hit.X() + shift.X());
227  hit.SetY(hit.Y() + shift.Y());
228  }
229 
230  void ARICHDigitizerModule::endRun()
231  {
232  }
233 
234  void ARICHDigitizerModule::terminate()
235  {
236 
237  }
238 
240 } // end Belle2 namespace
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::ARICHSimHit::getLocalPosition
TVector2 getLocalPosition() const
Get local position of hit (in module coordinates)
Definition: ARICHSimHit.h:79
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ARICHSimHit::getModuleID
int getModuleID() const
Get ID number of module that registered hit.
Definition: ARICHSimHit.h:76
Belle2::Module
Base class for Modules.
Definition: Module.h:74
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
Belle2::ARICHSimHit
Class ARICHSimHit - Geant4 simulated hit for ARICH.
Definition: ARICHSimHit.h:39
Belle2::asicChannel
record to be used to store ASIC info
Definition: CDCCrossTalkClasses.h:23
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::ARICHDigitizerModule
ARICH digitizer module.
Definition: ARICHDigitizerModule.h:43
Belle2::ARICHSimHit::getGlobalTime
float getGlobalTime() const override
Get global time of hit.
Definition: ARICHSimHit.h:82
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120