Belle II Software  release-05-02-19
SVDSimpleClusterizerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Giulia Casarosa *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <svd/modules/svdReconstruction/SVDSimpleClusterizerModule.h>
12 
13 #include <framework/datastore/DataStore.h>
14 #include <framework/datastore/RelationArray.h>
15 #include <framework/datastore/RelationIndex.h>
16 #include <framework/logging/Logger.h>
17 
18 #include <svd/geometry/SensorInfo.h>
19 #include <svd/dataobjects/SVDEventInfo.h>
20 
21 using namespace std;
22 using namespace Belle2;
23 using namespace Belle2::SVD;
24 
25 
26 //-----------------------------------------------------------------
27 // Register the Module
28 //-----------------------------------------------------------------
29 REG_MODULE(SVDSimpleClusterizer)
30 
31 //-----------------------------------------------------------------
32 // Implementation
33 //-----------------------------------------------------------------
34 
36  m_cutSeed(5.0), m_cutAdjacent(3.0), m_sizeHeadTail(3), m_cutCluster(0), m_useDB(true)
37 {
38  //Set module properties
39  setDescription("Clusterize SVDRecoDigits fitted by the Center of Gravity estimator");
40  setPropertyFlags(c_ParallelProcessingCertified);
41 
42  // 1. Collections.
43  addParam("RecoDigits", m_storeRecoDigitsName,
44  "SVDRecoDigits collection name", string(""));
45  addParam("Clusters", m_storeClustersName,
46  "SVDCluster collection name", string(""));
47  addParam("SVDTrueHits", m_storeTrueHitsName,
48  "TrueHit collection name", string(""));
49  addParam("MCParticles", m_storeMCParticlesName,
50  "MCParticles collection name", string(""));
51  addParam("ShaperDigits", m_storeShaperDigitsName,
52  "SVDShaperDigits collection name",
53  string(""));//NOTE: This collection is not directly accessed in this module, but indirectly accessed through SimpleClusterCandidate to get clustered samples.
54 
55  // 2. Clustering
56  addParam("AdjacentSN", m_cutAdjacent,
57  "SN for digits to be considered for clustering", m_cutAdjacent);
58  addParam("SeedSN", m_cutSeed,
59  "SN for digits to be considered as seed", m_cutSeed);
60  addParam("HeadTailSize", m_sizeHeadTail,
61  "Cluster size at which to switch to Analog head tail algorithm", m_sizeHeadTail);
62  addParam("ClusterSN", m_cutCluster,
63  "minimum value of the SNR of the cluster", m_cutCluster);
64  addParam("timeAlgorithm", m_timeAlgorithm,
65  " int to choose time algorithm: 0 = 6-sample CoG (default for 6-sample acquisition mode), 1 = 3-sample CoG (default for 3-sample acquisition mode), 2 = 3-sample ELS",
66  m_timeAlgorithm);
67  addParam("Calibrate3SampleWithEventT0", m_calibrate3SampleWithEventT0,
68  " if true returns the calibrated time instead of the raw time for 3-sample time algorithms",
69  m_calibrate3SampleWithEventT0);
70  addParam("useDB", m_useDB,
71  "if false use clustering module parameters", m_useDB);
72 
73 
74 }
75 
76 void SVDSimpleClusterizerModule::initialize()
77 {
78  //Register collections
79  m_storeClusters.registerInDataStore(m_storeClustersName, DataStore::c_ErrorIfAlreadyRegistered);
80  m_storeDigits.isRequired(m_storeRecoDigitsName);
81  m_storeTrueHits.isOptional(m_storeTrueHitsName);
82  m_storeMCParticles.isOptional(m_storeMCParticlesName);
83 
84  RelationArray relClusterDigits(m_storeClusters, m_storeDigits);
85  RelationArray relClusterTrueHits(m_storeClusters, m_storeTrueHits);
86  RelationArray relClusterMCParticles(m_storeClusters, m_storeMCParticles);
87  RelationArray relDigitTrueHits(m_storeDigits, m_storeTrueHits);
88  RelationArray relDigitMCParticles(m_storeDigits, m_storeMCParticles);
89 
90  relClusterDigits.registerInDataStore();
91  //Relations to simulation objects only if the ancestor relations exist
92  if (relDigitTrueHits.isOptional())
93  relClusterTrueHits.registerInDataStore();
94  if (relDigitMCParticles.isOptional())
95  relClusterMCParticles.registerInDataStore();
96 
97  //Store names to speed up creation later
98  m_storeClustersName = m_storeClusters.getName();
99  m_storeRecoDigitsName = m_storeDigits.getName();
100  m_storeTrueHitsName = m_storeTrueHits.getName();
101  m_storeMCParticlesName = m_storeMCParticles.getName();
102 
103  m_relClusterRecoDigitName = relClusterDigits.getName();
104  m_relClusterTrueHitName = relClusterTrueHits.getName();
105  m_relClusterMCParticleName = relClusterMCParticles.getName();
106  m_relRecoDigitTrueHitName = relDigitTrueHits.getName();
107  m_relRecoDigitMCParticleName = relDigitMCParticles.getName();
108 
109  // Report:
110  B2DEBUG(1, "SVDSimpleClusterizer Parameters (in default system unit, *=cannot be set directly):");
111 
112  B2DEBUG(1, " 1. COLLECTIONS:");
113  B2DEBUG(1, " --> MCParticles: " << DataStore::arrayName<MCParticle>(m_storeMCParticlesName));
114  B2DEBUG(1, " --> SVDRecoDigits: " << DataStore::arrayName<SVDRecoDigit>(m_storeRecoDigitsName));
115  B2DEBUG(1, " --> SVDClusters: " << DataStore::arrayName<SVDCluster>(m_storeClustersName));
116  B2DEBUG(1, " --> SVDTrueHits: " << DataStore::arrayName<SVDTrueHit>(m_storeTrueHitsName));
117  B2DEBUG(1, " --> DigitMCRel: " << m_relRecoDigitMCParticleName);
118  B2DEBUG(1, " --> ClusterMCRel: " << m_relClusterMCParticleName);
119  B2DEBUG(1, " --> ClusterDigitRel: " << m_relClusterRecoDigitName);
120  B2DEBUG(1, " --> DigitTrueRel: " << m_relRecoDigitTrueHitName);
121  B2DEBUG(1, " --> ClusterTrueRel: " << m_relClusterTrueHitName);
122  B2DEBUG(1, " 2. CLUSTERING:");
123  B2DEBUG(1, " --> Neighbour cut: " << m_cutAdjacent);
124  B2DEBUG(1, " --> Seed cut: " << m_cutSeed);
125  B2DEBUG(1, " --> Size HeadTail: " << m_sizeHeadTail);
126 }
127 
128 
129 
130 void SVDSimpleClusterizerModule::event()
131 {
132  int nDigits = m_storeDigits.getEntries();
133  if (nDigits == 0)
134  return;
135 
136  m_storeClusters.clear();
137 
138  RelationArray relClusterMCParticle(m_storeClusters, m_storeMCParticles,
139  m_relClusterMCParticleName);
140  if (relClusterMCParticle) relClusterMCParticle.clear();
141 
142  RelationArray relClusterDigit(m_storeClusters, m_storeDigits,
143  m_relClusterRecoDigitName);
144  if (relClusterDigit) relClusterDigit.clear();
145 
146  RelationArray relClusterTrueHit(m_storeClusters, m_storeTrueHits,
147  m_relClusterTrueHitName);
148  if (relClusterTrueHit) relClusterTrueHit.clear();
149 
150 
151  if (m_useDB) {
152  m_cutSeed = m_ClusterCal.getMinSeedSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
153  m_cutAdjacent = m_ClusterCal.getMinAdjSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
154  m_cutCluster = m_ClusterCal.getMinClusterSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
155  }
156  //create a dummy cluster just to start
157  SimpleClusterCandidate clusterCandidate(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip(),
158  m_sizeHeadTail, m_cutSeed, m_cutAdjacent, m_cutCluster, m_timeAlgorithm, m_storeShaperDigitsName, m_storeRecoDigitsName);
159 
160  //loop over the SVDRecoDigits
161  int i = 0;
162  while (i < nDigits) {
163 
164  //retrieve the VxdID, sensor and cellID of the current RecoDigit
165  VxdID thisSensorID = m_storeDigits[i]->getSensorID();
166  bool thisSide = m_storeDigits[i]->isUStrip();
167  int thisCellID = m_storeDigits[i]->getCellID();
168 
169  if (m_useDB) {
170  m_cutSeed = m_ClusterCal.getMinSeedSNR(thisSensorID, thisSide);
171  m_cutAdjacent = m_ClusterCal.getMinAdjSNR(thisSensorID, thisSide);
172  m_cutCluster = m_ClusterCal.getMinClusterSNR(thisSensorID, thisSide);
173  }
174 
175  //Ignore digits with insufficient signal
176  float thisNoise = m_NoiseCal.getNoiseInElectrons(thisSensorID, thisSide, thisCellID);
177  float thisCharge = m_storeDigits[i]->getCharge();
178  B2DEBUG(10, "Noise = " << thisNoise << " e-, Charge = " << thisCharge);
179 
180  if ((float)thisCharge / thisNoise < m_cutAdjacent) {
181  i++;
182  continue;
183  }
184 
185  //this strip has a sufficient S/N
186  stripInCluster aStrip;
187  aStrip.recoDigitIndex = i;
188  aStrip.charge = thisCharge;
189  aStrip.cellID = thisCellID;
190  aStrip.noise = thisNoise;
191  //this is the 6-sample CoG time and will be used to compute the 6-sample CoG cluster time, it will not be used for in 3-sample time algorithms:
192  aStrip.time = m_storeDigits[i]->getTime();
193  aStrip.timeError = m_storeDigits[i]->getTimeError();
194 
195  //try to add the strip to the existing cluster
196  if (! clusterCandidate.add(thisSensorID, thisSide, aStrip)) {
197 
198  //if the strip is not added, write the cluster, if present and good:
199  if (clusterCandidate.size() > 0) {
200  clusterCandidate.finalizeCluster();
201  if (clusterCandidate.isGoodCluster()) {
202  writeClusters(clusterCandidate);
203  }
204  }
205 
206  //prepare for the next cluster:
207  clusterCandidate = SimpleClusterCandidate(thisSensorID, thisSide, m_sizeHeadTail, m_cutSeed, m_cutAdjacent, m_cutCluster,
208  m_timeAlgorithm,
209  m_storeShaperDigitsName, m_storeRecoDigitsName);
210 
211  //start another cluster:
212  if (! clusterCandidate.add(thisSensorID, thisSide, aStrip))
213  B2WARNING("this state is forbidden!!");
214 
215  }
216  i++;
217  } //exit loop on RecoDigits
218 
219  //write the last cluster, if good
220  if (clusterCandidate.size() > 0) {
221  clusterCandidate.finalizeCluster();
222  if (clusterCandidate.isGoodCluster())
223  writeClusters(clusterCandidate);
224  }
225 
226  B2DEBUG(1, "Number of clusters: " << m_storeClusters.getEntries());
227 }
228 
229 
230 void SVDSimpleClusterizerModule::writeClusters(SimpleClusterCandidate cluster)
231 {
232 
233  RelationArray relClusterDigit(m_storeClusters, m_storeDigits, m_relClusterRecoDigitName);
234 
235  RelationArray relClusterMCParticle(m_storeClusters, m_storeMCParticles, m_relClusterMCParticleName);
236  RelationArray relClusterTrueHit(m_storeClusters, m_storeTrueHits, m_relClusterTrueHitName);
237 
238  RelationIndex<SVDRecoDigit, MCParticle> relDigitMCParticle(m_storeDigits, m_storeMCParticles, m_relRecoDigitMCParticleName);
239  RelationIndex<SVDRecoDigit, SVDTrueHit> relDigitTrueHit(m_storeDigits, m_storeTrueHits, m_relRecoDigitTrueHitName);
240 
241 
242  VxdID sensorID = cluster.getSensorID();
243  bool isU = cluster.isUSide();
244  float seedCharge = cluster.getSeedCharge();
245  float charge = cluster.getCharge();
246  float size = cluster.size();
247  float SNR = cluster.getSNR();
248  float position = cluster.getPosition();
249  float positionError = m_ClusterCal.getCorrectedClusterPositionError(sensorID, isU, size, cluster.getPositionError());
250  //this is the 6-sample CoG time, it will not be used for in 3-sample time algorithms:
251  float time = cluster.getTime();
252  float timeError = cluster.getTimeError();
253  int firstFrame = cluster.getFirstFrame();
254 
255  //first check SVDEventInfo name
256  StoreObjPtr<SVDEventInfo> temp_eventinfo("SVDEventInfo");
257  std::string m_svdEventInfoName = "SVDEventInfo";
258  if (!temp_eventinfo.isValid())
259  m_svdEventInfoName = "SVDEventInfoSim";
260  StoreObjPtr<SVDEventInfo> eventinfo(m_svdEventInfoName);
261  if (!eventinfo) B2ERROR("No SVDEventInfo!");
262 
263  //depending on the algorithm, time contains different information:
264  //6-sample CoG (0): this is the calibrated time already
265  //3-sample CoG (1) or ELS (2) this is the raw time, you need to calibrate:
266  //It is possile to get the uncalibrated 3-sample raw time here
267  //to get the 6-sample raw time there is an option in SVDCoGTimeEstimatorModule
268  float caltime = time;
269  if (m_timeAlgorithm == 1 and m_calibrate3SampleWithEventT0)
270  caltime = m_3CoGTimeCal.getCorrectedTime(sensorID, isU, -1, time, -1);
271  else if (m_timeAlgorithm == 2 and m_calibrate3SampleWithEventT0)
272  caltime = m_3ELSTimeCal.getCorrectedTime(sensorID, isU, -1, time, -1);
273 
274  // last step:
275  // shift cluster time by TB time AND by FirstFrame ( FF = 0 for the 6-sample CoG Time)
276  // the relative shift between 3- and 6-sample DAQ is also corrected
277  // NOTE: this shift is removed in the SVDTimeCalibrationCollector in the CAF
278  time = eventinfo->getTimeInFTSWReference(caltime, firstFrame);
279 
280  // Store Cluster into Datastore
281  m_storeClusters.appendNew(SVDCluster(
282  sensorID, isU, position, positionError, time, timeError, charge, seedCharge, size, SNR, -1, firstFrame
283  ));
284 
285  //register relation between RecoDigit and Cluster
286  int clsIndex = m_storeClusters.getEntries() - 1;
287 
288  map<int, float> mc_relations;
289  map<int, float> truehit_relations;
290 
291  vector<pair<int, float> > digit_weights;
292  digit_weights.reserve(size);
293 
294  std::vector<stripInCluster> strips = cluster.getStripsInCluster();
295 
296  for (auto strip : strips) {
297 
298  //Fill map with MCParticle relations
299  if (relDigitMCParticle) {
300  typedef const RelationIndex<SVDRecoDigit, MCParticle>::Element relMC_type;
301  for (relMC_type& mcRel : relDigitMCParticle.getElementsFrom(m_storeDigits[strip.recoDigitIndex])) {
302  //negative weights are from ignored particles, we don't like them and
303  //thus ignore them :D
304  if (mcRel.weight < 0) continue;
305  mc_relations[mcRel.indexTo] += mcRel.weight;
306  };
307  };
308  //Fill map with SVDTrueHit relations
309  if (relDigitTrueHit) {
310  typedef const RelationIndex<SVDRecoDigit, SVDTrueHit>::Element relTrueHit_type;
311  for (relTrueHit_type& trueRel : relDigitTrueHit.getElementsFrom(m_storeDigits[strip.recoDigitIndex])) {
312  //negative weights are from ignored particles, we don't like them and
313  //thus ignore them :D
314  if (trueRel.weight < 0) continue;
315  truehit_relations[trueRel.indexTo] += trueRel.weight;
316  };
317  };
318 
319  digit_weights.push_back(make_pair(strip.recoDigitIndex, strip.charge));
320  }
321 
322 
323  //Create Relations to this Digit
324  if (!mc_relations.empty()) {
325  relClusterMCParticle.add(clsIndex, mc_relations.begin(), mc_relations.end());
326  }
327  if (!truehit_relations.empty()) {
328  relClusterTrueHit.add(clsIndex, truehit_relations.begin(), truehit_relations.end());
329  }
330 
331  relClusterDigit.add(clsIndex, digit_weights.begin(), digit_weights.end());
332 }
Belle2::SVD::SimpleClusterCandidate::size
int size() const
return the cluster size (number of strips of the cluster
Definition: SimpleClusterCandidate.h:185
Belle2::RelationArray::clear
void clear() override
Clear all elements from the relation.
Definition: RelationArray.h:279
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::SVD::stripInCluster
structure containing the relevant informations of eachstrip of the cluster
Definition: SimpleClusterCandidate.h:38
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::SVD::stripInCluster::charge
float charge
strip charge
Definition: SimpleClusterCandidate.h:40
Belle2::SVD::stripInCluster::recoDigitIndex
int recoDigitIndex
index of the reco digit
Definition: SimpleClusterCandidate.h:39
Belle2::SVD::SimpleClusterCandidate
Class representing a cluster candidate during simple clustering of the SVD.
Definition: SimpleClusterCandidate.h:50
Belle2::RelationIndex
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
Definition: RelationIndex.h:84
Belle2::StoreAccessorBase::isOptional
bool isOptional(const std::string &name="")
Tell the DataStore about an optional input.
Definition: StoreAccessorBase.h:95
Belle2::StoreAccessorBase::getName
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Definition: StoreAccessorBase.h:130
Belle2::StoreAccessorBase::registerInDataStore
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Definition: StoreAccessorBase.h:54
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RelationIndexContainer::Element
Element type for the index.
Definition: RelationIndexContainer.h:62
Belle2::RelationIndex::getElementsFrom
range_from getElementsFrom(const FROM *from) const
Return a range of all elements pointing from the given object.
Definition: RelationIndex.h:157
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::SVD
Namespace to encapsulate code needed for simulation and reconstrucion of the SVD.
Definition: GeoSVDCreator.h:35
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::SVD::SVDSimpleClusterizerModule
SVDSimpleClusterizerModule: The SVD SimpleClusterizer.
Definition: SVDSimpleClusterizerModule.h:52
Belle2::SVD::SimpleClusterCandidate::isGoodCluster
bool isGoodCluster()
return true if the cluster candidate can be promoted to cluster
Definition: SimpleClusterCandidate.cc:209
Belle2::SVD::SimpleClusterCandidate::finalizeCluster
void finalizeCluster()
compute the position, time and their error of the cluster
Definition: SimpleClusterCandidate.cc:108
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::SVD::stripInCluster::time
float time
6-sample CoG strip time
Definition: SimpleClusterCandidate.h:43
Belle2::SVD::stripInCluster::timeError
float timeError
6-sample CoG strip time error
Definition: SimpleClusterCandidate.h:44
Belle2::SVD::stripInCluster::cellID
int cellID
strip cellID
Definition: SimpleClusterCandidate.h:42
Belle2::SVD::SimpleClusterCandidate::add
bool add(VxdID vxdID, bool isUside, struct stripInCluster &aStrip)
Add a Strip to the current cluster.
Definition: SimpleClusterCandidate.cc:76
Belle2::SVD::stripInCluster::noise
float noise
strip noise
Definition: SimpleClusterCandidate.h:41
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120