Belle II Software  release-05-01-25
SVDCrossTalkFinderModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: James Webb *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <svd/modules/svdCrossTalkFinder/SVDCrossTalkFinderModule.h>
12 #include <svd/modules/svdCrossTalkFinder/SVDCrossTalkFinderHelperFunctions.h>
13 
14 #include <vxd/geometry/GeoCache.h>
15 
16 using namespace std;
17 using namespace Belle2;
18 
19 REG_MODULE(SVDCrossTalkFinder)
20 
21 SVDCrossTalkFinderModule::SVDCrossTalkFinderModule() : Module()
22 
23 {
24  setDescription("Detect SVDRecoDigits created from cross-talk present in the origami sensors");
25 
27 
28  addParam("SVDRecoDigits", m_svdRecoDigitsName,
29  "SVDRecoDigit collection name", string(""));
30 
31  addParam("createCalibrationPayload", m_createCalibrationPayload,
32  "Create cross-talk strip channel payload", true);
33 
34  addParam("outputFilename", m_outputFilename,
35  "Filename of root file for calibration payload", std::string("crossTalkStripsCalibration.root"));
36 
37  addParam("uSideOccupancyFactor", m_uSideOccupancyFactor,
38  "Multiple of the average occupancy for high occupancy strip classification", 2);
39 
40  addParam("vSideOccupancyFactor", m_vSideOccupancyFactor,
41  "Multiple of the average occupancy for high occupancy strip classification", 2);
42 
43  addParam("nAPVFactor", m_nAPVFactor,
44  "Number of APV chips with at least one high occupancy strips in event required for cross talk flag", 30);
45 
46 }
47 
48 
50 {
51 
53  std::string sensorName;
54  TH1F* sensorHist;
55  //Prepare histograms for payload.
57  for (auto& layers : geo.getLayers(VXD::SensorInfoBase::SVD)) {
58  for (auto& ladders : geo.getLadders(layers)) {
59  for (auto& sensors : geo.getSensors(ladders)) {
60  for (int side = 0; side <= 1; side++) {
61  occupancyPDFName(sensors, side, sensorName);
62  if (m_sensorHistograms.count(sensorName) == 0) {
63  if (layers.getLayerNumber() == 3 or side == 1) {
64  sensorHist = new TH1F(sensorName.c_str(), "", 768, 0, 768);
65  } else {
66  sensorHist = new TH1F(sensorName.c_str(), "", 512, 0, 512);
67  }
68  m_sensorHistograms[sensorName] = sensorHist;
69 
70  }
71 
72  }
73  }
74  }
75  }
76  }
77 
78  // prepare storeArray
80  m_svdEventInfo.isRequired();
81 
82 }
83 
85 {
86 
87  vector<std::string> highOccChips_uSide ;
88  vector<std::string> highOccChips_vSide ;
89  vector<std::string> highOccChipsStripNum_uSide ;
90  vector<std::string> highOccChipsStripNum_vSide ;
91  vector<std::string> clusterStrips_uSide ;
92  vector<std::string> clusterStrips_vSide ;
93  vector<std::string> clusterChips_uSide ;
94  vector<int> strips_uSide;
95  vector<int> strips_vSide;
96 
97  for (auto& svdRecoDigit : m_svdRecoDigits) {
98  //Remove L3 and +fw sensors, not affected by cross-talk
99  if (svdRecoDigit.getSensorID().getLayerNumber() == 3 or svdRecoDigit.getSensorID().getSensorNumber() == 1) {
100  continue;
101  }
102 
103  int side = svdRecoDigit.isUStrip();
104  std::string sensorName;
105  occupancyPDFName(svdRecoDigit.getSensorID(), side, sensorName);
106 
107  double sensorAverage = 0.;
108  calculateAverage(svdRecoDigit.getSensorID(), sensorAverage, side);
109  int stripID = svdRecoDigit.getCellID();
110  std::string sensorStripNum = sensorName + "." + std::to_string(stripID);
111  double stripOccupancy = m_OccupancyCal.getOccupancy(svdRecoDigit.getSensorID(), side, stripID);
112  //Clustering only works assuming digits are ordered.//
113 
114  if (side == 1 && stripOccupancy > (m_uSideOccupancyFactor * sensorAverage)) {
115 
116  int adjacentStrip = 0;
117  if (!strips_uSide.empty()) {
118  adjacentStrip = stripID - strips_uSide.back();
119  }
120  strips_uSide.push_back(stripID);
121  if (!highOccChips_uSide.empty() && sensorName == highOccChips_uSide.back() && adjacentStrip == 1) {
122  clusterChips_uSide.push_back(sensorName);
123  if (clusterStrips_uSide.empty() || clusterStrips_uSide.back() != sensorStripNum) {
124  clusterStrips_uSide.push_back(highOccChipsStripNum_uSide.back());
125  }
126  clusterStrips_uSide.push_back(sensorStripNum);
127  }
128  highOccChipsStripNum_uSide.push_back(sensorStripNum);
129  highOccChips_uSide.push_back(sensorName);
130  }
131 
132  if (side == 0 && stripOccupancy > (m_vSideOccupancyFactor * sensorAverage)) {
133  int adjacentStrip = 0;
134  if (!strips_vSide.empty()) {
135  adjacentStrip = stripID - strips_vSide.back();
136  }
137  strips_vSide.push_back(stripID);
138  if (!highOccChips_vSide.empty() && sensorName == highOccChips_vSide.back() && adjacentStrip == 1) {
139  if (clusterStrips_vSide.empty() || clusterStrips_vSide.back() != sensorStripNum) {
140  clusterStrips_vSide.push_back(highOccChipsStripNum_vSide.back());
141  }
142  clusterStrips_vSide.push_back(sensorStripNum);
143  }
144  highOccChipsStripNum_vSide.push_back(sensorStripNum);
145  highOccChips_vSide.push_back(sensorName);
146  }
147 
148  }
149 
150  std::sort(clusterChips_uSide.begin(), clusterChips_uSide.end());
151  clusterChips_uSide.erase(unique(clusterChips_uSide.begin(), clusterChips_uSide.end()), clusterChips_uSide.end());
152  int numberOfClusterChips = std::distance(clusterChips_uSide.begin(), std::unique(clusterChips_uSide.begin(),
153  clusterChips_uSide.end()));
154 
155 // Crosstalk events flagged using u-side
156  if (numberOfClusterChips > m_nAPVFactor) {
157  m_svdEventInfo->setCrossTalk(true);
158  for (auto& svdRecoDigit : m_svdRecoDigits) {
159  std::string sensorID = svdRecoDigit.getSensorID();
160  std::string digitID = sensorID + "." + std::to_string(svdRecoDigit.isUStrip());
161  std::string stripID = digitID + "." + std::to_string(svdRecoDigit.getCellID());
163  if (std::find(clusterStrips_uSide.begin(), clusterStrips_uSide.end(), stripID) != clusterStrips_uSide.end()) {
164  std::string sensorName;
165  occupancyPDFName(svdRecoDigit.getSensorID(), svdRecoDigit.isUStrip(), sensorName);
166  auto xTalkStrip = m_sensorHistograms.at(sensorName);
167  //Only fill bin once
168  if (xTalkStrip->GetBinContent(svdRecoDigit.getCellID()) < 1.) xTalkStrip->Fill(svdRecoDigit.getCellID(), true);
169  }
170  if (std::find(clusterStrips_vSide.begin(), clusterStrips_vSide.end(), stripID) != clusterStrips_vSide.end()) {
171  std::string sensorName;
172  occupancyPDFName(svdRecoDigit.getSensorID(), svdRecoDigit.isUStrip(), sensorName);
173  auto xTalkStrip = m_sensorHistograms.at(sensorName);
174  if (xTalkStrip->GetBinContent(svdRecoDigit.getCellID()) < 1.) xTalkStrip->Fill(svdRecoDigit.getCellID(), true);
175  }
176 
177  }
178 
179  }//reco digit loop
180  }
181 
182 }
183 
185 {
186  B2INFO("SVDCrossTalkFinderModule::terminate");
187 
189  m_histogramFile = new TFile(m_outputFilename.c_str(), "RECREATE");
190  m_histogramFile->cd();
191  std::string sensorName;
193  for (auto& layers : geo.getLayers(VXD::SensorInfoBase::SVD)) {
194  for (auto& ladders : geo.getLadders(layers)) {
195  for (auto& sensors : geo.getSensors(ladders)) {
196  for (int side = 0; side <= 1; side++) {
197 
198  occupancyPDFName(sensors, side, sensorName);
199  auto sensorOnMap = m_sensorHistograms.at(sensorName);
200  sensorOnMap->Write();
201 
202  }
203  }
204  }
205  }
206  m_histogramFile->Close();
207  }
208 }
209 
210 
211 void SVDCrossTalkFinderModule::calculateAverage(const VxdID& sensorID, double& mean, int side)
212 {
213  double nBins = 0;
214  if (side == 1) nBins = 768; //U-side 768 channels
215  else nBins = 512;
216  double count = 0;
217  for (int i = 0; i < nBins; i++) {
218  count += m_OccupancyCal.getOccupancy(sensorID, side, i);
219  }
220  mean = count / nBins;
221 }
Belle2::SVDCrossTalkFinderModule::m_svdRecoDigitsName
std::string m_svdRecoDigitsName
Function to calculate sensor average occupancy.
Definition: SVDCrossTalkFinderModule.h:73
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:216
Belle2::SVDCrossTalkFinderModule::m_histogramFile
TFile * m_histogramFile
Filename of root file containing cross-talk strip calibration payload.
Definition: SVDCrossTalkFinderModule.h:94
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::Module::c_ParallelProcessingCertified
@ 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:82
Belle2::VXD::GeoCache::getLayers
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition: GeoCache.cc:177
Belle2::SVDCrossTalkFinderModule::m_OccupancyCal
SVDOccupancyCalibrations m_OccupancyCal
SVDOccupancy calibrations db object.
Definition: SVDCrossTalkFinderModule.h:98
Belle2::SVDCrossTalkFinderModule::terminate
virtual void terminate() override
Final output.
Definition: SVDCrossTalkFinderModule.cc:184
Belle2::VXD::GeoCache::getSensors
const std::set< Belle2::VxdID > & getSensors(Belle2::VxdID ladder) const
Return a set of all sensor IDs belonging to a given ladder.
Definition: GeoCache.cc:205
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::Module::setPropertyFlags
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:210
Belle2::SVDOccupancyCalibrations::getOccupancy
float getOccupancy(const VxdID &sensorID, const bool &isU, const unsigned short &strip) const
This is the method for getting the occupancy, or the deviation from the average, still to be decided.
Definition: SVDOccupancyCalibrations.h:71
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:215
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDCrossTalkFinderModule::m_svdEventInfo
StoreObjPtr< SVDEventInfo > m_svdEventInfo
The storeObject for svdEventInfo.
Definition: SVDCrossTalkFinderModule.h:82
Belle2::occupancyPDFName
void occupancyPDFName(const VxdID &sensor, int side, std::string &PDFName)
Function to maintain common naming convention between calibration occupancy file generation and occup...
Definition: SVDCrossTalkFinderHelperFunctions.h:37
Belle2::VXD::SensorInfoBase::SVD
@ SVD
SVD Sensor.
Definition: SensorInfoBase.h:45
Belle2::SVDCrossTalkFinderModule::m_vSideOccupancyFactor
int m_vSideOccupancyFactor
Parameter to define high occupancy strips (some multiple above sensor average occupancy)
Definition: SVDCrossTalkFinderModule.h:86
Belle2::SVDCrossTalkFinderModule::m_svdRecoDigits
StoreArray< SVDRecoDigit > m_svdRecoDigits
The storeArray for svdRecoDigits.
Definition: SVDCrossTalkFinderModule.h:76
Belle2::SVDCrossTalkFinderModule::m_nAPVFactor
int m_nAPVFactor
Parameter to define high occupancy strips (some multiple above sensor average occupancy)
Definition: SVDCrossTalkFinderModule.h:88
Belle2::SVDCrossTalkFinderModule::event
virtual void event() override
Event.
Definition: SVDCrossTalkFinderModule.cc:84
Belle2::Module::addParam
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:562
Belle2::SVDCrossTalkFinderModule::m_sensorHistograms
std::map< std::string, TH1F * > m_sensorHistograms
Pointer to root TFile containing histograms for calibration payload.
Definition: SVDCrossTalkFinderModule.h:95
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:41
Belle2::SVDCrossTalkFinderModule::m_createCalibrationPayload
bool m_createCalibrationPayload
Parameter to set number of sensors with possible cross-talk clusters required for event flagging.
Definition: SVDCrossTalkFinderModule.h:90
Belle2::VXD::GeoCache::getLadders
const std::set< Belle2::VxdID > & getLadders(Belle2::VxdID layer) const
Return a set of all ladder IDs belonging to a given layer.
Definition: GeoCache.cc:194
Belle2::SVDCrossTalkFinderModule::m_outputFilename
std::string m_outputFilename
If true module will produce and write-out payload for SVDCrossTalkStripsCalibrations.
Definition: SVDCrossTalkFinderModule.h:92
Belle2::SVDCrossTalkFinderModule::initialize
virtual void initialize() override
Init the module.
Definition: SVDCrossTalkFinderModule.cc:49