Belle II Software development
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/logging/Logger.h>
16#include <TH2F.h>
17#include <TFitResult.h>
18#include <TF1.h>
19#include <TGraphAsymmErrors.h>
20#include <TMath.h>
21
22using namespace Belle2;
23using namespace CDC;
24using namespace TrackFindingCDC;
26{
28 " -------------------------- Wire Efficiency Estimation Algorithm -------------------------\n"
29 );
30}
31
33{
34 B2INFO("Creating efficiencies for every layer");
35
36 unsigned short wireID;
37 unsigned short layerID;
38 float z;
39 bool isFound;
40
41 auto efftree = getObjectPtr<TTree>("efftree");
42 efftree->SetBranchAddress("wireID", &wireID);
43 efftree->SetBranchAddress("layerID", &layerID);
44 efftree->SetBranchAddress("z", &z);
45 efftree->SetBranchAddress("isFound", &isFound);
46
47 // Set ranges and names for the TEfficiency objects
48 B2INFO("Building empty efficiency objects");
49 static const CDCWireTopology& wireTopology = CDCWireTopology::getInstance();
50 B2INFO("Wire Topology found");
51 static const CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
52 B2INFO("CDC Geometry found");
53
54 for (const CDCWireLayer& wireLayer : wireTopology.getWireLayers()) {
55 unsigned short layerNo = wireLayer.getICLayer();
56 //B2INFO("Got layer " << layerNo);
57 //std::string m_nameOfLayer = std::string("effLayer_").append(std::to_string(layerNo));
58 //std::string m_titleOfLayer = std::string("Efficiency of wires in layer ").append(std::to_string(layerNo));
59 //B2INFO("Built names for " << layerNo);
60
61 unsigned short nzbins = 30 - layerNo / 7;
62 unsigned short nwidbins = cdcgeo.nWiresInLayer(layerNo);
63 double widbins[2] = { -0.5, cdcgeo.nWiresInLayer(layerNo) - 0.5};
64 double zbins[2] = {wireLayer.getBackwardZ(), wireLayer.getForwardZ()};
65
66 TEfficiency* effInLayer = new TEfficiency(Form("effLayer_%d", layerNo), Form("Hit efficiency pf L%d ; z (cm) ; IWire ", layerNo),
67 nzbins, zbins[0], zbins[1], nwidbins, widbins[0], widbins[1]);
68 m_efficiencyList->Add(effInLayer);
69 B2INFO("Teff for layer " << layerNo << " was successfully listed.");
70 }
71 TFile* outputCollection = new TFile("LayerEfficiencies.root", "RECREATE");
72
73 // loop over entries in the tree and build the TEfficiencies
74 B2INFO("Filling the efficiencies");
75 const Long64_t nEntries = efftree->GetEntries();
76 B2INFO("Number of entries in tree: " << nEntries);
77 for (Long64_t i = 0; i < nEntries; i++) {
78 efftree->GetEntry(i);
79 TEfficiency* efficiencyInLayer = (TEfficiency*)m_efficiencyList->At(layerID);
80 efficiencyInLayer->Fill(isFound, z, wireID);
81 }
82 m_efficiencyList->Write();
83 B2INFO("TEfficiencies successfully filled.");
84 outputCollection->Close();
85 B2INFO("TEfficiencies successfully saved");
86 // setting 1000 entries as the minimum value to consider "enough data", but this might need to be tailored in the future.
87 if (nEntries > 1000) return true;
88 else return false;
89}
90
92{
93 B2INFO("Beginning detection of bad wires");
94 static const CDCWireTopology& wireTopology = CDCWireTopology::getInstance();
95 B2INFO("Wire Topology found");
96
97 for (const CDCWireLayer& wireLayer : wireTopology.getWireLayers()) {
98 unsigned short layerNo = wireLayer.getICLayer();
99 B2INFO("Checking layer " << layerNo);
100
101 // need to use casting here because GetPassedHistogram assumes it returns TH1
102 auto efficiencyInLayer = (TEfficiency*)m_efficiencyList->At(layerNo);
103 auto passed = (TH2F*)efficiencyInLayer->GetPassedHistogram();
104 auto total = (TH2F*)efficiencyInLayer->GetTotalHistogram();
105
106 // Ignoring layers that have no hits at all
107 if (!total->GetEntries()) continue;
108
109 // prepare fitting functions
110 double minFitRange = wireLayer.getBackwardZ() + 30;
111 double maxFitRange = wireLayer.getForwardZ() - 30;
112 TF1* constantFunction = new TF1("constantFunction", "pol0", minFitRange, maxFitRange);
113
114 // Estimate average efficiency of the layer as function of z
115 auto passedProjectedX = passed->ProjectionX("projectionx1");
116 auto totalProjectedX = total->ProjectionX("projectionx2");
117 TEfficiency* efficiencyProjectedX = new TEfficiency(*passedProjectedX, *totalProjectedX);
118
119 unsigned short minFitBin = passedProjectedX->FindBin(minFitRange) + 1;
120 unsigned short maxFitBin = passedProjectedX->FindBin(maxFitRange) - 1;
121
122 // Estimate average efficiency of each wire in the layer
123 auto passedProjectedY = passed->ProjectionY("projectiony1", minFitRange, maxFitRange);
124 auto totalProjectedY = total->ProjectionY("projectiony2", minFitRange, maxFitRange);
125 TEfficiency* efficiencyProjectedY = new TEfficiency(*passedProjectedY, *totalProjectedY);
126
127 // Estimate average efficiency of the whole layer by summing up averages of individual wires
128 float totalAverage = 0;
129 int nonZeroWires = 0;
130 for (int i = 0; i <= passedProjectedY->GetNbinsX(); ++i) {
131 float efficiencyAtBin = efficiencyProjectedY->GetEfficiency(i);
132 if (efficiencyAtBin > 0.2) {
133 totalAverage += efficiencyAtBin;
134 nonZeroWires++;
135 }
136 }
137 // safeguard in case all wires are 0
138 nonZeroWires > 0 ? totalAverage /= nonZeroWires : totalAverage == 0;
139
140 TGraphAsymmErrors* graphEfficiencyProjected = efficiencyProjectedX->CreateGraph();
141
142 B2INFO("Successfully determined average properties of layer " << layerNo);
143 B2INFO("Looping over wires");
144
145 // Due to the way histograms are binned, every bin center should correspond a wire ID.
146 for (int i = 0; i <= passed->GetNbinsY(); ++i) {
147
148 // project out 1 wire
149 auto singleWirePassed = passed->ProjectionX("single wire projection passed", i, i);
150 auto singleWireTotal = total->ProjectionX("single wire projection total", i, i);
151
152 // if not a single hit within the central area of wire then mark it as bad
153 if (!singleWirePassed->Integral(minFitBin, maxFitBin)) {
154 double wireID = passed->GetYaxis()->GetBinCenter(i);
155 unsigned short maxWires = CDCGeometryPar::Instance().nWiresInLayer(layerNo);
156 if (wireID < 0 || wireID >= maxWires) {
157 B2ERROR("Invalid wireID: " << wireID << " for LayerID: " << layerNo
158 << ". Max wires: " << maxWires);
159 continue; // remove invalid wireID
160 }
161 m_badWireList->setWire(layerNo, round(wireID), 0);
162 continue;
163 }
164 // constructs efficiency of one wire
165 TEfficiency* singleWireEfficiency = new TEfficiency(*singleWirePassed, *singleWireTotal);
166 TGraphAsymmErrors* graphSingleWireEfficiency = singleWireEfficiency->CreateGraph();
167
168 // fits the efficiency with a constant
169 TFitResultPtr singleWireFitResults = graphSingleWireEfficiency->Fit(constantFunction, "SQR");
170 double singleWireEfficiencyFromFit = (singleWireFitResults.Get())->Value(0);
171 double singleWireUpErrorFromFit = (singleWireFitResults.Get())->UpperError(0);
172
173 // applies a chi test on the wire
174 float p_value = chiTest(graphSingleWireEfficiency, graphEfficiencyProjected, minFitRange, maxFitRange);
175 // check three conditions and decide if wire should be marked as bad
176 bool averageCondition = (0.95 * totalAverage) > (singleWireEfficiencyFromFit + 3 * singleWireUpErrorFromFit);
177 bool pvalueCondition = p_value < 0.01;
178 bool generalCondition = singleWireEfficiencyFromFit < 0.4;
179 if (generalCondition || (averageCondition && pvalueCondition)) {
180 double wireID = passed->GetYaxis()->GetBinCenter(i);
181 m_badWireList->setWire(layerNo, round(wireID), singleWireEfficiencyFromFit);
182 }
183 }
184 B2INFO("Bad wires for " << layerNo << " recorded");
185 }
186 m_badWireList->outputToFile("wireFile.txt");
187 saveCalibration(m_badWireList, "CDCBadWires");
188 B2INFO("Bad wire list successfully saved.");
189}
190
191double WireEfficiencyAlgorithm::chiTest(TGraphAsymmErrors* graph1, TGraphAsymmErrors* graph2, double minValue, double maxValue)
192{
193
194 // Vars to perform the chi test
195 double chi = 0;
196 unsigned short ndof = 0;
197
198
199 int numOfEntries1 = graph1->GetN();
200 int numOfEntries2 = graph2->GetN();
201
202 // loop over entries in both graphs, index by index.
203
204 for (int index1 = 0; index1 < numOfEntries1; ++index1) {
205 // TGraph values are usually not listed in increasing order. Need to check that they are within min/max range for comparison.
206 if (graph1->GetX()[index1] < minValue) continue;
207 if (graph1->GetX()[index1] > maxValue) continue;
208 for (int index2 = 0; index2 < numOfEntries2; ++index2) {
209 if (graph1->GetX()[index1] == graph2->GetX()[index2]) {
210 // this is broken up just for readability
211 double chiNumerator = pow(graph1->GetY()[index1] - graph2->GetY()[index2], 2);
212 double chiDenominator = pow(graph1->GetErrorYhigh(index1), 2) + pow(graph1->GetErrorYlow(index1), 2);
213 chi += chiNumerator / chiDenominator;
214 ndof++;
215 continue;
216 }
217 }
218 }
219 return TMath::Prob(chi, ndof);
220}
221
223{
224 const auto exprun = getRunList()[0];
225 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
226 updateDBObjPtrs(1, exprun.second, exprun.first);
227
228 B2INFO("Creating CDCGeometryPar object");
230
231 bool enoughEventsInInput;
232 enoughEventsInInput = buildEfficiencies();
233
235
236 if (enoughEventsInInput) return c_OK;
237 else return c_NotEnoughData;
238}
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.