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