Belle II Software  release-05-02-19
TOPAlignmentCollectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <top/modules/collectors/TOPAlignmentCollectorModule.h>
12 #include <top/geometry/TOPGeometryPar.h>
13 #include <top/reconstruction/TOPconfigure.h>
14 #include <top/reconstruction/TOPtrack.h>
15 
16 // framework aux
17 #include <framework/logging/Logger.h>
18 
19 // root
20 #include <TTree.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TRandom.h>
24 
25 using namespace std;
26 
27 namespace Belle2 {
33  using namespace TOP;
34 
35  //-----------------------------------------------------------------
36  // Register module
37  //-----------------------------------------------------------------
38 
39  REG_MODULE(TOPAlignmentCollector)
40 
41  //-----------------------------------------------------------------
42  // Implementation
43  //-----------------------------------------------------------------
44 
46  {
47  // set module description and processing properties
48  setDescription("A collector for geometrical alignment of a TOP module with dimuons or Bhabhas. Iterative alignment procedure (NIMA 876 (2017) 260-264) is run here, algorithm just collects the results.");
49  setPropertyFlags(c_ParallelProcessingCertified);
50 
51  // module parameters
52  addParam("minBkgPerBar", m_minBkgPerBar,
53  "Minimal number of background photons per module", 0.0);
54  addParam("scaleN0", m_scaleN0,
55  "Scale factor for figure-of-merit N0", 1.0);
56  addParam("targetModule", m_targetMid,
57  "Module to be aligned. Must be 1 <= Mid <= 16.");
58  addParam("maxFails", m_maxFails,
59  "Maximum number of consecutive failed iterations before resetting the procedure", 100);
60  addParam("sample", m_sample,
61  "sample type: one of dimuon or bhabha", std::string("dimuon"));
62  addParam("deltaEcms", m_deltaEcms,
63  "c.m.s energy window (half size) if sample is dimuon or bhabha", 0.1);
64  addParam("dr", m_dr, "cut on POCA in r", 1.0);
65  addParam("dz", m_dz, "cut on POCA in abs(z)", 2.0);
66  addParam("minZ", m_minZ, "minimal local z of extrapolated hit", -130.0);
67  addParam("maxZ", m_maxZ, "maximal local z of extrapolated hit", 130.0);
68  addParam("stepPosition", m_stepPosition, "step size for translations", 1.0);
69  addParam("stepAngle", m_stepAngle, "step size for rotations", 0.01);
70  addParam("stepTime", m_stepTime, "step size for t0", 0.05);
71  addParam("stepRefind", m_stepRefind,
72  "step size for scaling of refractive index (dn/n)", 0.005);
73  addParam("gridSize", m_gridSize,
74  "size of a 2D grid for time-of-propagation averaging in analytic PDF: "
75  "[number of emission points along track, number of Cerenkov angles]. "
76  "No grid used if list is empty.", m_gridSize);
77  std::string names;
78  auto align = TOPalign();
79  for (const auto& parName : align.getParameterNames()) names += parName + ", ";
80  names.pop_back();
81  names.pop_back();
82  addParam("parInit", m_parInit,
83  "initial parameter values in the order [" + names + "]. "
84  "If list is too short, the missing ones are set to 0.", m_parInit);
85  auto parFixed = m_parFixed;
86  addParam("parFixed", m_parFixed, "list of names of parameters to be fixed. "
87  "Valid names are: " + names, parFixed);
88 
89  }
90 
91 
92  void TOPAlignmentCollectorModule::prepare()
93  {
94  // input collections
95 
96  m_digits.isRequired();
97  m_tracks.isRequired();
98  m_extHits.isRequired();
99  m_recBunch.isRequired();
100 
101  // check if target module ID is valid
102 
103  const auto* geo = TOPGeometryPar::Instance()->getGeometry();
104  if (!geo->isModuleIDValid(m_targetMid)) {
105  B2FATAL("Target module ID = " << m_targetMid << " is invalid. Exiting...");
106  }
107 
108  // set track selector
109 
110  if (m_sample == "dimuon" or m_sample == "bhabha") {
111  m_selector = TrackSelector(m_sample);
112  m_selector.setDeltaEcms(m_deltaEcms);
113  m_selector.setCutOnPOCA(m_dr, m_dz);
114  m_selector.setCutOnLocalZ(m_minZ, m_maxZ);
115  } else {
116  B2ERROR("Invalid sample type '" << m_sample << "'");
117  }
118 
119  // set alignment objects
120 
121  for (int set = 0; set < c_numSets; set++) {
122  auto align = TOPalign();
123  align.setModuleID(m_targetMid);
124  align.setSteps(m_stepPosition, m_stepAngle, m_stepTime, m_stepRefind);
125  if (m_gridSize.size() == 2) {
126  align.setGrid(m_gridSize[0], m_gridSize[1]);
127  B2INFO("TOPAligner: grid for time-of-propagation averaging is set");
128  }
129  align.setParameters(m_parInit);
130  for (const auto& parName : m_parFixed) {
131  align.fixParameter(parName);
132  }
133  m_align.push_back(align);
134  m_countFails.push_back(0);
135  }
136 
137  // configure detector in reconstruction code
138 
139  TOPconfigure config;
140 
141  // create and register output histograms and ntuples
142 
143  int numModules = geo->getNumModules();
144  auto h1 = new TH2F("tracks_per_slot", "Number of tracks per slot and sample",
145  numModules, 0.5, numModules + 0.5, c_numSets, 0, c_numSets);
146  h1->SetXTitle("slot number");
147  h1->SetYTitle("sample number");
148  registerObject<TH2F>("tracks_per_slot", h1);
149 
150  for (int slot = 1; slot <= numModules; slot++) {
151  std::string slotName = "_s" + to_string(slot);
152  std::string slotTitle = "(slot " + to_string(slot) + ")";
153 
154  std::string hname = "local_z" + slotName;
155  std::string title = "Distribution of tracks along bar " + slotTitle;
156  auto h2 = new TH1F(hname.c_str(), title.c_str(), 100, -150.0, 150.0);
157  h2->SetXTitle("local z [cm]");
158  registerObject<TH1F>(hname, h2);
159 
160  hname = "cth_vs_p" + slotName;
161  title = "Track momentum " + slotTitle;
162  auto h3 = new TH2F(hname.c_str(), title.c_str(), 100, 0.0, 7.0, 100, -1.0, 1.0);
163  h3->SetXTitle("p [GeV/c]");
164  h3->SetYTitle("cos #theta");
165  registerObject<TH2F>(hname, h3);
166 
167  hname = "poca_xy" + slotName;
168  title = "Track POCA in x-y " + slotTitle;
169  auto h4 = new TH2F(hname.c_str(), title.c_str(), 100, -m_dr, m_dr, 100, -m_dr, m_dr);
170  h4->SetXTitle("x [cm]");
171  h4->SetYTitle("y [cm]");
172  registerObject<TH2F>(hname, h4);
173 
174  hname = "poca_z" + slotName;
175  title = "Track POCA in z " + slotTitle;
176  auto h5 = new TH1F(hname.c_str(), title.c_str(), 100, -m_dz, m_dz);
177  h5->SetXTitle("z [cm]");
178  registerObject<TH1F>(hname, h5);
179 
180  hname = "Ecms" + slotName;
181  title = "Track c.m.s. energy " + slotTitle;
182  auto h6 = new TH1F(hname.c_str(), title.c_str(), 100, 5.1, 5.4);
183  h6->SetXTitle("E_{cms} [GeV]");
184  registerObject<TH1F>(hname, h6);
185 
186  hname = "charge" + slotName;
187  title = "Charge of track " + slotTitle;
188  auto h7 = new TH1F(hname.c_str(), title.c_str(), 3, -1.5, 1.5);
189  h7->SetXTitle("charge");
190  registerObject<TH1F>(hname, h7);
191 
192  hname = "timeHits" + slotName;
193  title = "Photon time distribution " + slotTitle;
194  auto h8 = new TH2F(hname.c_str(), title.c_str(), 512, 0, 512, 200, 0, 50);
195  h8->SetXTitle("channel number");
196  h8->SetYTitle("time [ns]");
197  registerObject<TH2F>(hname, h8);
198 
199  hname = "numPhot" + slotName;
200  title = "Number of photons " + slotTitle;
201  auto h9 = new TH1F(hname.c_str(), title.c_str(), 100, 0, 100);
202  h9->SetXTitle("photons per track");
203  registerObject<TH1F>(hname, h9);
204  }
205 
206  for (int set = 0; set < c_numSets; set++) {
207  std::string name = "alignTree" + to_string(set);
208  m_treeNames.push_back(name);
209  auto alignTree = new TTree(name.c_str(), "TOP alignment results");
210  alignTree->Branch("ModuleId", &m_targetMid);
211  alignTree->Branch("iter", &m_iter);
212  alignTree->Branch("ntrk", &m_ntrk);
213  alignTree->Branch("errorCode", &m_errorCode);
214  alignTree->Branch("iterPars", &m_vAlignPars);
215  alignTree->Branch("iterParsErr", &m_vAlignParsErr);
216  alignTree->Branch("valid", &m_valid);
217  alignTree->Branch("numPhot", &m_numPhot);
218  alignTree->Branch("x", &m_x);
219  alignTree->Branch("y", &m_y);
220  alignTree->Branch("z", &m_z);
221  alignTree->Branch("p", &m_p);
222  alignTree->Branch("theta", &m_theta);
223  alignTree->Branch("phi", &m_phi);
224  alignTree->Branch("r_poca", &m_pocaR);
225  alignTree->Branch("z_poca", &m_pocaZ);
226  alignTree->Branch("x_poca", &m_pocaX);
227  alignTree->Branch("y_poca", &m_pocaY);
228  alignTree->Branch("Ecms", &m_cmsE);
229  alignTree->Branch("charge", &m_charge);
230  alignTree->Branch("PDG", &m_PDG);
231  registerObject<TTree>(name, alignTree);
232  }
233 
234  }
235 
236 
237  void TOPAlignmentCollectorModule::collect()
238  {
239  // bunch must be reconstructed
240 
241  if (not m_recBunch.isValid()) return;
242  if (not m_recBunch->isReconstructed()) return;
243 
244  // add photons
245 
246  TOPalign::clearData();
247 
248  for (const auto& digit : m_digits) {
249  if (digit.getHitQuality() == TOPDigit::c_Good)
250  TOPalign::addData(digit.getModuleID(), digit.getPixelID(), digit.getTime(),
251  digit.getTimeError());
252  }
253 
254  TOPalign::setPhotonYields(m_minBkgPerBar, m_scaleN0);
255 
256  // track-by-track iterations
257 
258  for (const auto& track : m_tracks) {
259 
260  // construct TOPtrack from mdst track
261  TOPtrack trk(&track);
262  if (not trk.isValid()) continue;
263 
264  // skip if track not hitting target module
265  if (trk.getModuleID() != m_targetMid) continue;
266 
267  // track selection
268  if (not m_selector.isSelected(trk)) continue;
269 
270  // generate sub-sample number
271  int sub = gRandom->Integer(c_numSets);
272  auto& align = m_align[sub];
273  auto& countFails = m_countFails[sub];
274  const auto& name = m_treeNames[sub];
275  auto h1 = getObjectPtr<TH2F>("tracks_per_slot");
276  h1->Fill(trk.getModuleID(), sub);
277 
278  // do an iteration
279  int err = align.iterate(trk, m_selector.getChargedStable());
280  m_iter++;
281 
282  // check number of consecutive failures, and in case reset
283  if (err == 0) {
284  countFails = 0;
285  } else if (countFails <= m_maxFails) {
286  countFails++;
287  } else {
288  B2INFO("Reached maximum allowed number of failed iterations. "
289  "Resetting TOPalign object of set = " << sub << " at iter = " << m_iter);
290  align.reset();
291  countFails = 0;
292  }
293 
294  // get new parameter values and estimated errors
295  m_vAlignPars = align.getParameters();
296  m_vAlignParsErr = align.getErrors();
297  m_ntrk = align.getNumUsedTracks();
298  m_errorCode = err;
299  m_valid = align.isValid();
300  m_numPhot = align.getNumOfPhotons();
301 
302  // set other ntuple variables
303  const auto& localPosition = m_selector.getLocalPosition();
304  m_x = localPosition.X();
305  m_y = localPosition.Y();
306  m_z = localPosition.Z();
307  const auto& localMomentum = m_selector.getLocalMomentum();
308  m_p = localMomentum.Mag();
309  m_theta = localMomentum.Theta();
310  m_phi = localMomentum.Phi();
311  const auto& pocaPosition = m_selector.getPOCAPosition();
312  m_pocaR = pocaPosition.Perp();
313  m_pocaZ = pocaPosition.Z();
314  m_pocaX = pocaPosition.X();
315  m_pocaY = pocaPosition.Y();
316  m_cmsE = m_selector.getCMSEnergy();
317  m_charge = trk.getCharge();
318  m_PDG = trk.getPDGcode();
319 
320  // fill output tree
321  auto alignTree = getObjectPtr<TTree>(name);
322  alignTree->Fill();
323 
324  // fill control histograms
325  std::string slotName = "_s" + to_string(m_targetMid);
326  auto h2 = getObjectPtr<TH1F>("local_z" + slotName);
327  h2->Fill(m_z);
328  auto h3 = getObjectPtr<TH2F>("cth_vs_p" + slotName);
329  h3->Fill(m_p, cos(m_theta));
330  auto h4 = getObjectPtr<TH2F>("poca_xy" + slotName);
331  h4->Fill(m_pocaX, m_pocaY);
332  auto h5 = getObjectPtr<TH1F>("poca_z" + slotName);
333  h5->Fill(m_pocaZ);
334  auto h6 = getObjectPtr<TH1F>("Ecms" + slotName);
335  h6->Fill(m_cmsE);
336  auto h7 = getObjectPtr<TH1F>("charge" + slotName);
337  h7->Fill(m_charge);
338  auto h8 = getObjectPtr<TH2F>("timeHits" + slotName);
339  for (const auto& digit : m_digits) {
340  if (digit.getHitQuality() != TOPDigit::c_Good) continue;
341  if (digit.getModuleID() != m_targetMid) continue;
342  h8->Fill(digit.getChannel(), digit.getTime());
343  }
344  auto h9 = getObjectPtr<TH1F>("numPhot" + slotName);
345  h9->Fill(m_numPhot);
346 
347  } // tracks
348 
349  }
350 
352 } // Belle2
Belle2::TOP::TOPtrack::getPDGcode
int getPDGcode() const
Return PDG code.
Definition: TOPtrack.h:199
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::TOP::TOPtrack
Class to hold reconstructed track, interface to fortran.
Definition: TOPtrack.h:42
Belle2::TOP::TOPtrack::isValid
bool isValid() const
Check if track is properly defined.
Definition: TOPtrack.h:87
Belle2::TOP::TOPtrack::getModuleID
int getModuleID() const
Return module ID.
Definition: TOPtrack.h:217
Belle2::TOPAlignmentCollectorModule
Collector for geometrical alignment of a TOP module with dimuons or Bhabhas.
Definition: TOPAlignmentCollectorModule.h:48
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TOP::TOPalign
Alignment of a TOP module: core of the method is coded in fortran.
Definition: TOPalign.h:38
Belle2::TOP::TOPconfigure
Configure TOP geometry for reconstruction: provides interface to fortran code.
Definition: TOPconfigure.h:36
Belle2::TOP::TrackSelector
Utility for the track selection - used in various calibration modules.
Definition: TrackSelector.h:37
Belle2::TOP::TOPtrack::getCharge
int getCharge() const
Return charge.
Definition: TOPtrack.h:205