Belle II Software prerelease-10-00-00a
WireEfficiencyAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <cdc/calibration/WireEfficiencyAlgorithm.h>
10#include <calibration/CalibrationAlgorithm.h>
11#include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
12#include <tracking/trackFindingCDC/topology/CDCWireLayer.h>
13#include <cdc/geometry/CDCGeometryPar.h>
14#include <cdc/dbobjects/CDCBadWires.h>
15#include <framework/database/IntervalOfValidity.h>
16#include <framework/logging/Logger.h>
17#include <TH2F.h>
18#include <TFitResult.h>
19#include <TH1F.h>
20#include <TF1.h>
21#include <TGraphAsymmErrors.h>
22#include <TMath.h>
23#include <math.h>
24
25using namespace Belle2;
26using namespace CDC;
27using namespace TrackFindingCDC;
29{
31 " -------------------------- Wire Efficiency Estimation Algorithm -------------------------\n"
32 );
33}
34
36{
37 B2INFO("Creating efficiencies for every layer");
38
39 unsigned short wireID;
40 unsigned short layerID;
41 float z;
42 bool isFound;
43
44 auto efftree = getObjectPtr<TTree>("efftree");
45 efftree->SetBranchAddress("wireID", &wireID);
46 efftree->SetBranchAddress("layerID", &layerID);
47 efftree->SetBranchAddress("z", &z);
48 efftree->SetBranchAddress("isFound", &isFound);
49
50 // Set ranges and names for the TEfficiency objects
51 B2INFO("Building empty efficiency objects");
52 static const CDCWireTopology& wireTopology = CDCWireTopology::getInstance();
53 B2INFO("Wire Topology found");
54 static const CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
55 B2INFO("CDC Geometry found");
56
57 for (const CDCWireLayer& wireLayer : wireTopology.getWireLayers()) {
58 unsigned short layerNo = wireLayer.getICLayer();
59 //B2INFO("Got layer " << layerNo);
60 //std::string m_nameOfLayer = std::string("effLayer_").append(std::to_string(layerNo));
61 //std::string m_titleOfLayer = std::string("Efficiency of wires in layer ").append(std::to_string(layerNo));
62 //B2INFO("Built names for " << layerNo);
63
64 unsigned short nzbins = 30 - layerNo / 7;
65 unsigned short nwidbins = cdcgeo.nWiresInLayer(layerNo);
66 double widbins[2] = { -0.5, cdcgeo.nWiresInLayer(layerNo) - 0.5};
67 double zbins[2] = {wireLayer.getBackwardZ(), wireLayer.getForwardZ()};
68
69 TEfficiency* effInLayer = new TEfficiency(Form("effLayer_%d", layerNo), Form("Hit efficiency pf L%d ; z (cm) ; IWire ", layerNo),
70 nzbins, zbins[0], zbins[1], nwidbins, widbins[0], widbins[1]);
71 m_efficiencyList->Add(effInLayer);
72 B2INFO("Teff for layer " << layerNo << " was successfully listed.");
73 }
74 TFile* outputCollection = new TFile("LayerEfficiencies.root", "RECREATE");
75
76 // loop over entries in the tree and build the TEfficiencies
77 B2INFO("Filling the efficiencies");
78 const Long64_t nEntries = efftree->GetEntries();
79 B2INFO("Number of entries in tree: " << nEntries);
80 for (Long64_t i = 0; i < nEntries; i++) {
81 efftree->GetEntry(i);
82 TEfficiency* efficiencyInLayer = (TEfficiency*)m_efficiencyList->At(layerID);
83 efficiencyInLayer->Fill(isFound, z, wireID);
84 }
85 m_efficiencyList->Write();
86 B2INFO("TEfficiencies successfully filled.");
87 outputCollection->Close();
88 B2INFO("TEfficiencies successfully saved");
89 // setting 1000 entries as the minimum value to consider "enough data", but this might need to be tailored in the future.
90 if (nEntries > 1000) return true;
91 else return false;
92}
93
95{
96 B2INFO("Beginning detection of bad wires");
97 static const CDCWireTopology& wireTopology = CDCWireTopology::getInstance();
98 B2INFO("Wire Topology found");
99
100 for (const CDCWireLayer& wireLayer : wireTopology.getWireLayers()) {
101 unsigned short layerNo = wireLayer.getICLayer();
102 B2INFO("Checking layer " << layerNo);
103
104 // need to use casting here because GetPassedHistogram assumes it returns TH1
105 auto efficiencyInLayer = (TEfficiency*)m_efficiencyList->At(layerNo);
106 auto passed = (TH2F*)efficiencyInLayer->GetPassedHistogram();
107 auto total = (TH2F*)efficiencyInLayer->GetTotalHistogram();
108
109 // Ignoring layers that have no hits at all
110 if (!total->GetEntries()) continue;
111
112 // prepare fitting functions
113 double minFitRange = wireLayer.getBackwardZ() + 30;
114 double maxFitRange = wireLayer.getForwardZ() - 30;
115 TF1* constantFunction = new TF1("constantFunction", "pol0", minFitRange, maxFitRange);
116
117 // Estimate average efficiency of the layer as function of z
118 auto passedProjectedX = passed->ProjectionX("projectionx1");
119 auto totalProjectedX = total->ProjectionX("projectionx2");
120 TEfficiency* efficiencyProjectedX = new TEfficiency(*passedProjectedX, *totalProjectedX);
121
122 unsigned short minFitBin = passedProjectedX->FindBin(minFitRange) + 1;
123 unsigned short maxFitBin = passedProjectedX->FindBin(maxFitRange) - 1;
124
125 // Estimate average efficiency of each wire in the layer
126 auto passedProjectedY = passed->ProjectionY("projectiony1", minFitRange, maxFitRange);
127 auto totalProjectedY = total->ProjectionY("projectiony2", minFitRange, maxFitRange);
128 TEfficiency* efficiencyProjectedY = new TEfficiency(*passedProjectedY, *totalProjectedY);
129
130 // Estimate average efficiency of the whole layer by summing up averages of individual wires
131 float totalAverage = 0;
132 int nonZeroWires = 0;
133 for (int i = 0; i <= passedProjectedY->GetNbinsX(); ++i) {
134 float efficiencyAtBin = efficiencyProjectedY->GetEfficiency(i);
135 if (efficiencyAtBin > 0.2) {
136 totalAverage += efficiencyAtBin;
137 nonZeroWires++;
138 }
139 }
140 // safeguard in case all wires are 0
141 nonZeroWires > 0 ? totalAverage /= nonZeroWires : totalAverage == 0;
142
143 TGraphAsymmErrors* graphEfficiencyProjected = efficiencyProjectedX->CreateGraph();
144
145 B2INFO("Successfully determined average properties of layer " << layerNo);
146 B2INFO("Looping over wires");
147
148 // Due to the way histograms are binned, every bin center should correspond a wire ID.
149 for (int i = 0; i <= passed->GetNbinsY(); ++i) {
150
151 // project out 1 wire
152 auto singleWirePassed = passed->ProjectionX("single wire projection passed", i, i);
153 auto singleWireTotal = total->ProjectionX("single wire projection total", i, i);
154
155 // if not a single hit within the central area of wire then mark it as bad
156 if (!singleWirePassed->Integral(minFitBin, maxFitBin)) {
157 double wireID = passed->GetYaxis()->GetBinCenter(i);
158 unsigned short maxWires = CDCGeometryPar::Instance().nWiresInLayer(layerNo);
159 if (wireID < 0 || wireID >= maxWires) {
160 B2ERROR("Invalid wireID: " << wireID << " for LayerID: " << layerNo
161 << ". Max wires: " << maxWires);
162 continue; // remove invalid wireID
163 }
164 m_badWireList->setWire(layerNo, round(wireID), 0);
165 continue;
166 }
167 // constructs efficiency of one wire
168 TEfficiency* singleWireEfficiency = new TEfficiency(*singleWirePassed, *singleWireTotal);
169 TGraphAsymmErrors* graphSingleWireEfficiency = singleWireEfficiency->CreateGraph();
170
171 // fits the efficiency with a constant
172 TFitResultPtr singleWireFitResults = graphSingleWireEfficiency->Fit(constantFunction, "SQR");
173 double singleWireEfficiencyFromFit = (singleWireFitResults.Get())->Value(0);
174 double singleWireUpErrorFromFit = (singleWireFitResults.Get())->UpperError(0);
175
176 // applies a chi test on the wire
177 float p_value = chiTest(graphSingleWireEfficiency, graphEfficiencyProjected, minFitRange, maxFitRange);
178 // check three conditions and decide if wire should be marked as bad
179 bool averageCondition = (0.95 * totalAverage) > (singleWireEfficiencyFromFit + 3 * singleWireUpErrorFromFit);
180 bool pvalueCondition = p_value < 0.01;
181 bool generalCondition = singleWireEfficiencyFromFit < 0.4;
182 if (generalCondition || (averageCondition && pvalueCondition)) {
183 double wireID = passed->GetYaxis()->GetBinCenter(i);
184 m_badWireList->setWire(layerNo, round(wireID), singleWireEfficiencyFromFit);
185 }
186 }
187 B2INFO("Bad wires for " << layerNo << " recorded");
188 }
189 m_badWireList->outputToFile("wireFile.txt");
190 saveCalibration(m_badWireList, "CDCBadWires");
191 B2INFO("Bad wire list successfully saved.");
192}
193
194double WireEfficiencyAlgorithm::chiTest(TGraphAsymmErrors* graph1, TGraphAsymmErrors* graph2, double minValue, double maxValue)
195{
196
197 // Vars to perform the chi test
198 double chi = 0;
199 unsigned short ndof = 0;
200
201
202 int numOfEntries1 = graph1->GetN();
203 int numOfEntries2 = graph2->GetN();
204
205 // loop over entries in both graphs, index by index.
206
207 for (int index1 = 0; index1 < numOfEntries1; ++index1) {
208 // TGraph values are usually not listed in increasing order. Need to check that they are within min/max range for comparison.
209 if (graph1->GetX()[index1] < minValue) continue;
210 if (graph1->GetX()[index1] > maxValue) continue;
211 for (int index2 = 0; index2 < numOfEntries2; ++index2) {
212 if (graph1->GetX()[index1] == graph2->GetX()[index2]) {
213 // this is broken up just for readability
214 double chiNumerator = pow(graph1->GetY()[index1] - graph2->GetY()[index2], 2);
215 double chiDenominator = pow(graph1->GetErrorYhigh(index1), 2) + pow(graph1->GetErrorYlow(index1), 2);
216 chi += chiNumerator / chiDenominator;
217 ndof++;
218 continue;
219 }
220 }
221 }
222 return TMath::Prob(chi, ndof);
223}
224
226{
227 const auto exprun = getRunList()[0];
228 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
229 updateDBObjPtrs(1, exprun.second, exprun.first);
230
231 B2INFO("Creating CDCGeometryPar object");
233
234 bool enoughEventsInInput;
235 enoughEventsInInput = buildEfficiencies();
236
238
239 if (enoughEventsInInput) return c_OK;
240 else return c_NotEnoughData;
241}
The Class for CDC Geometry Parameters.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
TList * m_efficiencyList
TList of efficiencies.
double chiTest(TGraphAsymmErrors *graph1, TGraphAsymmErrors *graph2, double minVale, double maxValue)
chitest
CDCBadWires * m_badWireList
BadWireList that willbe built.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
EResult calibrate() override
Run algo on data.
bool buildEfficiencies()
create 2D TEfficiency for each wire and return True if more than 1000 entries
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Class representing a sense wire layer in the central drift chamber.
Class representing the sense wire arrangement in the whole of the central drift chamber.
const std::vector< Belle2::TrackFindingCDC::CDCWireLayer > & getWireLayers() const
Getter for the underlying storing layer vector.
static CDCWireTopology & getInstance()
Getter for the singleton instance of the wire topology.
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.