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>
17 #include <framework/logging/Logger.h>
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);
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);
79 for (
const auto& parName : align.getParameterNames()) names += parName +
", ";
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);
92 void TOPAlignmentCollectorModule::prepare()
96 m_digits.isRequired();
97 m_tracks.isRequired();
98 m_extHits.isRequired();
99 m_recBunch.isRequired();
103 const auto* geo = TOPGeometryPar::Instance()->getGeometry();
104 if (!geo->isModuleIDValid(m_targetMid)) {
105 B2FATAL(
"Target module ID = " << m_targetMid <<
" is invalid. Exiting...");
110 if (m_sample ==
"dimuon" or m_sample ==
"bhabha") {
112 m_selector.setDeltaEcms(m_deltaEcms);
113 m_selector.setCutOnPOCA(m_dr, m_dz);
114 m_selector.setCutOnLocalZ(m_minZ, m_maxZ);
116 B2ERROR(
"Invalid sample type '" << m_sample <<
"'");
121 for (
int set = 0; set < c_numSets; set++) {
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");
129 align.setParameters(m_parInit);
130 for (
const auto& parName : m_parFixed) {
131 align.fixParameter(parName);
133 m_align.push_back(align);
134 m_countFails.push_back(0);
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);
150 for (
int slot = 1; slot <= numModules; slot++) {
151 std::string slotName =
"_s" + to_string(slot);
152 std::string slotTitle =
"(slot " + to_string(slot) +
")";
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);
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);
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);
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);
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);
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);
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);
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);
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);
237 void TOPAlignmentCollectorModule::collect()
241 if (not m_recBunch.isValid())
return;
242 if (not m_recBunch->isReconstructed())
return;
246 TOPalign::clearData();
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());
254 TOPalign::setPhotonYields(m_minBkgPerBar, m_scaleN0);
258 for (
const auto& track : m_tracks) {
262 if (not trk.
isValid())
continue;
268 if (not m_selector.isSelected(trk))
continue;
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");
279 int err = align.iterate(trk, m_selector.getChargedStable());
285 }
else if (countFails <= m_maxFails) {
288 B2INFO(
"Reached maximum allowed number of failed iterations. "
289 "Resetting TOPalign object of set = " << sub <<
" at iter = " << m_iter);
295 m_vAlignPars = align.getParameters();
296 m_vAlignParsErr = align.getErrors();
297 m_ntrk = align.getNumUsedTracks();
299 m_valid = align.isValid();
300 m_numPhot = align.getNumOfPhotons();
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();
321 auto alignTree = getObjectPtr<TTree>(name);
325 std::string slotName =
"_s" + to_string(m_targetMid);
326 auto h2 = getObjectPtr<TH1F>(
"local_z" + slotName);
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);
334 auto h6 = getObjectPtr<TH1F>(
"Ecms" + slotName);
336 auto h7 = getObjectPtr<TH1F>(
"charge" + slotName);
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());
344 auto h9 = getObjectPtr<TH1F>(
"numPhot" + slotName);