Belle II Software development
NoKickCutsEvalModule.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 <framework/datastore/StoreArray.h>
10#include <framework/datastore/RelationArray.h>
11#include <TFile.h>
12#include <tracking/modules/vxdtfRedesign/NoKickCutsEvalModule.h>
13
14#include <mdst/dataobjects/MCParticle.h>
15#include <svd/dataobjects/SVDTrueHit.h>
16
17#include <stdlib.h>
18#include <fstream>
19#include <string>
20#include "TH1.h"
21#include "TH3.h"
22#include "TF1.h"
23#include "TH2.h"
24#include "TF2.h"
25#include "TString.h"
26#include <algorithm>
27
28using namespace Belle2;
29
30REG_MODULE(NoKickCutsEval);
31
32NoKickCutsEvalModule::NoKickCutsEvalModule() : Module()
33{
34 setDescription("This module evaluate cuts necessary for the selction of reco tracks based on Multiple Scattering, NoKickRTSel");
35
37 addParam("useValidation", c_validationON,
38 "print in output file validation plot: track parameters distributions and cuts distributions", false);
39
41 addParam("useFitMethod", c_fitMethod, "apply the method of double-gaussian fit to evaluate the cuts", false);
42}
43
44void NoKickCutsEvalModule::initialize()
45{
46 m_histoLim.push_back(0.4 * c_multLimit);
47 m_histoLim.push_back(1. * c_multLimit);
48 m_histoLim.push_back(0.3 * c_multLimit);
49 m_histoLim.push_back(1. * c_multLimit);
50 m_histoLim.push_back(0.3 * c_multLimit);
51
52
53 for (int par = 0; par < c_nbinpar; par++) {
54 std::vector<std::vector<std::vector<std::vector<TH1F*>>>> histo_par;
55 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
56 std::vector<std::vector<std::vector<TH1F*>>> histo_lay1;
57 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
58 std::vector<std::vector<TH1F*>> histo_lay2;
59 for (int theta = 0; theta < c_nbint; theta++) {
60 std::vector<TH1F*> histo_theta;
61 for (int p = 0; p < c_nbinp; p++) {
62 histo_theta.push_back(new TH1F("histo_" + m_namePar.at(par) + Form("_layer%d-%d_theta%d_p%d", lay1, lay2, theta, p),
63 "histo_" + m_namePar.at(par) + Form("_layer%d-%d_theta%d_p%d", lay1, lay2, theta, p), c_nbin, -m_histoLim.at(par),
64 m_histoLim.at(par)));
65 }
66 histo_lay2.push_back(histo_theta);
67 histo_theta.clear();
68 }
69 histo_lay1.push_back(histo_lay2);
70 histo_lay2.clear();
71 }
72 histo_par.push_back(histo_lay1);
73 histo_lay1.clear();
74 }
75 m_histo.push_back(histo_par);
76 histo_par.clear();
77 }
78
80 StoreArray<SVDCluster> storeClusters("");
81 StoreArray<SVDTrueHit> storeTrueHits("");
82 StoreArray<MCParticle> storeMCParticles("");
83 StoreArray<RecoTrack> recoTracks("");
84
85 storeClusters.isRequired();
86 storeTrueHits.isRequired();
87 storeMCParticles.isRequired();
88 recoTracks.isRequired();
89
90 RelationArray relClusterTrueHits(storeClusters, storeTrueHits);
91 RelationArray relClusterMCParticles(storeClusters, storeMCParticles);
92 RelationArray recoTracksToMCParticles(recoTracks, storeMCParticles);
93
95 m_outputFile = new TFile("NoKickCuts.root", "RECREATE");
96}
97
98
99void NoKickCutsEvalModule::event()
100{
101 StoreArray<RecoTrack> recoTracks;
102
103 for (const RecoTrack& track : recoTracks) {
104 m_trackSel.hit8TrackBuilder(track);
105 std::vector<hitXP> XP8 = m_trackSel.m_8hitTrack;
106 bool PriorCut = m_trackSel.globalCut(XP8);
107 m_trackSel.m_8hitTrack.clear();
108 m_trackSel.m_hitXP.clear();
109 m_trackSel.m_setHitXP.clear();
110 if (!PriorCut) {m_globCounter++; continue;}
111
112 if (XP8.size() > 0) {
113 for (int i = 0; (i + 1) < (int)XP8.size(); i++) {
114 for (int par = 0; par < c_nbinpar; par++) {
115 int p = (int)((XP8.at(i).m_momentum0.R() - c_pmin) / c_pwidth);
116 if (p > c_nbinp - 1 || p < 0) {
117 m_pCounter++;
118 continue;
119 }
120 double sinTheta = abs(XP8.at(i).m_momentum0.Y()) / sqrt(pow(XP8.at(i).m_momentum0.Y(), 2) + pow(XP8.at(i).m_momentum0.Z(), 2));
121 int t = (int)((asin(sinTheta) - c_tmin) / c_pwidth);
122 if (t > c_nbint - 1 || t < 0) {
123 m_tCounter++;
124 continue;
125 }
126 double deltaPar = deltaParEval(XP8.at(i), XP8.at(i + 1), (NoKickCuts::EParameters)par);
127 if (deltaPar == c_over) continue;
128 m_histo.at(par).at(XP8.at(i).m_sensorLayer).at(XP8.at(i + 1).m_sensorLayer).at(t).at(p)->Fill(deltaPar);
129 if (i == 0) {
130 deltaPar = deltaParEval(XP8.at(i), XP8.at(i), (NoKickCuts::EParameters)par, true);
131 if (deltaPar == c_over)continue;
132 m_histo.at(par).at(0).at(XP8.at(i).m_sensorLayer).at(t).at(p)->Fill(deltaPar);
133 }
134 }
135 }
136 }
137 }
138}
139
140
141void NoKickCutsEvalModule::endRun()
142{
143 //-------------------------------FIT-EVALUATE THE CUTS---------------------------------------------------//
144
145 std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> cut_m;
146 std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> cut_M;
147
148 for (int par = 0; par < c_nbinpar; par++) {
149 std::vector<std::vector<std::vector<std::vector<double>>>> cut_M_par;
150 std::vector<std::vector<std::vector<std::vector<double>>>> cut_m_par;
151 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
152 std::vector<std::vector<std::vector<double>>> cut_M_lay1;
153 std::vector<std::vector<std::vector<double>>> cut_m_lay1;
154 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
155 std::vector<std::vector<double>> cut_M_lay2;
156 std::vector<std::vector<double>> cut_m_lay2;
157 for (int theta = 0; theta < c_nbint; theta++) {
158 std::vector<double> cut_M_theta;
159 std::vector<double> cut_m_theta;
160 for (int p = 0; p < c_nbinp; p++) {
161
162 //--------------first method to evaluate cuts, not used -------------------------//
163 if (c_fitMethod) {
164 TF1* fit_2gaus = new TF1("fit_2gaus", "[0]*TMath::Gaus(x,[1],[2], kTRUE)+[3]*TMath::Gaus(x,[4],[5],kTRUE)+[6]",
165 -m_histoLim.at(par), m_histoLim.at(par));
166
167 int bin0 = c_nbin / 2;
168 int nbin0 = m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin0);
169 double moltSigma1 = 0.5;
170 double moltSigma2 = 2;
171 double sigma1 = (double)moltSigma1 * (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetStdDev());
172 double sigma2 = (double)moltSigma2 * (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetStdDev());
173 double norm1 = sigma1 * sqrt(2 * M_PI) * 0.9 * nbin0;
174 double norm2 = sigma2 * sqrt(2 * M_PI) * 0.1 * nbin0;
175 double mean1 = (double)m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetMean();
176 double mean2 = (double)m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetMean();
177 double bkg = (double)m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(2);
178
179 fit_2gaus->SetParameters(norm1, mean1, sigma1, norm2, mean2, sigma2, bkg);
180 m_histo.at(par).at(lay1).at(lay2).at(theta).at(p) -> Fit(fit_2gaus, "", "", -m_histoLim.at(par), m_histoLim.at(par));
181 cut_M_theta.push_back(3 * (sqrt(fit_2gaus->GetParameter(2)*fit_2gaus->GetParameter(2))));
182 cut_m_theta.push_back(3 * (-sqrt(fit_2gaus->GetParameter(2)*fit_2gaus->GetParameter(2))));
183 }
184 //---------------END of first method ------------------------//
185
186 else {
187 //--------second method to evaluate (without fit), used-----------------//
188 double integral = m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->Integral();
189 double sum_M = m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(c_nbin + 1);
190 double sum_m = m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(0);
191 double percent = 1 - cutFunction(p, c_pwidth);
192
193 int bin_m = 0;
194 int bin_M = c_nbin + 1;
195 while (sum_m < integral * percent / 2) {
196 bin_m++;
197 sum_m = sum_m + m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin_m);
198 }
199 while (sum_M < integral * percent / 2) {
200 bin_M--;
201 sum_M = sum_M + m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(bin_M);
202 }
203 if (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetEntries() < 100) {
204 int filledBin_m = 0;
205 int filledBin_M = c_nbin + 1;
206 while (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(filledBin_m) == 0 && filledBin_m < (double) c_nbin / 2) {
207 filledBin_m++;
208 }
209 while (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinContent(filledBin_M) == 0 && filledBin_M > (double) c_nbin / 2) {
210 filledBin_M--;
211 }
212 bin_m = filledBin_m;
213 bin_M = filledBin_M;
214 }
215 cut_M_theta.push_back(m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinCenter(bin_M));
216 cut_m_theta.push_back(m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetBinCenter(bin_m));
217 }
218 //-----------------------END of second method ------------------------//
219
220 }
221 cut_M_lay2.push_back(cut_M_theta);
222 cut_M_theta.clear();
223 cut_m_lay2.push_back(cut_m_theta);
224 cut_m_theta.clear();
225 }
226 cut_M_lay1.push_back(cut_M_lay2);
227 cut_M_lay2.clear();
228 cut_m_lay1.push_back(cut_m_lay2);
229 cut_m_lay2.clear();
230 }
231 cut_M_par.push_back(cut_M_lay1);
232 cut_M_lay1.clear();
233 cut_m_par.push_back(cut_m_lay1);
234 cut_m_lay1.clear();
235 }
236 cut_M.push_back(cut_M_par);
237 cut_M_par.clear();
238 cut_m.push_back(cut_m_par);
239 cut_m_par.clear();
240 }
241
242 //------------------------------------------FIT THE CUTS --------------------------------//
243
245 std::vector<std::vector<std::vector<TH2F*>>> cut_M_histo;
246 std::vector<std::vector<std::vector<TH2F*>>> cut_m_histo;
247
248 for (int par = 0; par < c_nbinpar; par++) {
249 std::vector<std::vector<TH2F*>> cut_M_histo_par;
250 std::vector<std::vector<TH2F*>> cut_m_histo_par;
251 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
252 std::vector<TH2F*> cut_M_histo_lay1;
253 std::vector<TH2F*> cut_m_histo_lay1;
254 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
255 cut_M_histo_lay1.push_back(new TH2F("CUTS_M_" + m_namePar.at(par) + Form("layer%d_%d", lay1, lay2),
256 "CUTS_M_" + m_namePar.at(par) + Form("_layer%d_%d", lay1, lay2), c_nbinp, c_pmin, c_pmax, c_nbint, c_tmin, c_tmax));
257 cut_m_histo_lay1.push_back(new TH2F("CUTS_m_" + m_namePar.at(par) + Form("layer%d_%d", lay1, lay2),
258 "CUTS_m_" + m_namePar.at(par) + Form("_layer%d_%d", lay1, lay2), c_nbinp, c_pmin, c_pmax, c_nbint, c_tmin, c_tmax));
259 for (int theta = 1; theta <= c_nbint; theta++) {
260 for (int p = 1; p <= c_nbinp; p++) {
261 cut_M_histo_lay1.at(lay2)->SetBinContent(p, theta, cut_M.at(par).at(lay1).at(lay2).at(theta - 1).at(p - 1));
262 cut_m_histo_lay1.at(lay2)->SetBinContent(p, theta, cut_m.at(par).at(lay1).at(lay2).at(theta - 1).at(p - 1));
263 }
264 }
265 }
266 cut_M_histo_par.push_back(cut_M_histo_lay1);
267 cut_M_histo_lay1.clear();
268 cut_m_histo_par.push_back(cut_m_histo_lay1);
269 cut_m_histo_lay1.clear();
270 }
271 cut_M_histo.push_back(cut_M_histo_par);
272 cut_M_histo_par.clear();
273 cut_m_histo.push_back(cut_m_histo_par);
274 cut_m_histo_par.clear();
275 }
276
278 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_norm;
279 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_pow;
280 std::vector<std::vector<std::vector<std::vector<double>>>> cut_out_bkg;
281
282 for (int minmax = 0; minmax < 2; minmax++) {
283 std::vector<std::vector<std::vector<double>>> cut_out_norm_minmax;
284 std::vector<std::vector<std::vector<double>>> cut_out_pow_minmax;
285 std::vector<std::vector<std::vector<double>>> cut_out_bkg_minmax;
286 for (int par = 0; par < c_nbinpar; par++) {
287 std::vector<std::vector<double>> cut_out_norm_par;
288 std::vector<std::vector<double>> cut_out_pow_par;
289 std::vector<std::vector<double>> cut_out_bkg_par;
290 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
291 std::vector<double> cut_out_norm_lay1;
292 std::vector<double> cut_out_pow_lay1;
293 std::vector<double> cut_out_bkg_lay1;
294 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
295
296 TF2* fit_MS = new TF2("fit_MS", "[0]*1/(TMath::Power(x,[1])*TMath::Sqrt(TMath::Sin(y)))+[2]", c_pmin, c_pmax, c_tmin, c_tmax);
297
298 double norm = cut_M.at(par).at(lay1).at(lay2).at(1).at(1);
299 double Pow = 1;
300 double bkg = cut_M.at(par).at(lay1).at(lay2).at(1).at(c_nbinp - 1);
301 fit_MS->SetParameters(norm, Pow, bkg);
302
303 if (minmax == 0) {
304 cut_m_histo.at(par).at(lay1).at(lay2) -> Fit("fit_MS");
305
306 }
307 if (minmax == 1) {
308 cut_M_histo.at(par).at(lay1).at(lay2) -> Fit("fit_MS");
309
310 }
311 cut_out_norm_lay1.push_back(fit_MS->GetParameter(0));
312 cut_out_pow_lay1.push_back(fit_MS->GetParameter(1));
313 cut_out_bkg_lay1.push_back(fit_MS->GetParameter(2));
314 }
315 cut_out_norm_par.push_back(cut_out_norm_lay1);
316 cut_out_norm_lay1.clear();
317 cut_out_pow_par.push_back(cut_out_pow_lay1);
318 cut_out_pow_lay1.clear();
319 cut_out_bkg_par.push_back(cut_out_bkg_lay1);
320 cut_out_bkg_lay1.clear();
321 }
322 cut_out_norm_minmax.push_back(cut_out_norm_par);
323 cut_out_norm_par.clear();
324 cut_out_pow_minmax.push_back(cut_out_pow_par);
325 cut_out_pow_par.clear();
326 cut_out_bkg_minmax.push_back(cut_out_bkg_par);
327 cut_out_bkg_par.clear();
328 }
329 cut_out_norm.push_back(cut_out_norm_minmax);
330 cut_out_norm_minmax.clear();
331 cut_out_pow.push_back(cut_out_pow_minmax);
332 cut_out_pow_minmax.clear();
333 cut_out_bkg.push_back(cut_out_bkg_minmax);
334 cut_out_bkg_minmax.clear();
335 }
336
337//----------------------------------------------VALIDATION PLOTS ------------------------------//
338
339//-------------Some debugs lines-----------------//
340 if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::LogConfig::c_Debug, 30, PACKAGENAME())) {
341 for (int g = 0; g < c_nbinp; g++) {
342 double p_out_mom = g * c_pwidth + c_pmin;
343 double t_out_theta = c_twidth + c_tmin;
344 B2DEBUG(30, "momentum=" << p_out_mom);
345 B2DEBUG(30, "d0, 3-4, Min: " << cut_m.at(1).at(3).at(4).at(1).at(g));
346 B2DEBUG(30, "min cut (TH2F):" << cut_m_histo.at(1).at(3).at(4)->GetBinContent(g + 1, 2));
347 double norm_min = cut_out_norm.at(0).at(1).at(3).at(4);
348 double pow_min = cut_out_pow.at(0).at(1).at(3).at(4);
349 double bkg_min = cut_out_bkg.at(0).at(1).at(3).at(4);
350 B2DEBUG(30, "norm par min:" << norm_min);
351 B2DEBUG(30, "pow par min:" << pow_min);
352 B2DEBUG(30, "bkg par min:" << bkg_min);
353 B2DEBUG(30, "evaluate min cut:" << norm_min / (sqrt(sin(t_out_theta)) * pow(p_out_mom, pow_min)) + bkg_min);
354 B2DEBUG(30, "d0, 3-4, Max: " << cut_M.at(1).at(3).at(4).at(1).at(g));
355 B2DEBUG(30, "max cut (TH2F):" << cut_M_histo.at(1).at(3).at(4)->GetBinContent(g + 1, 2));
356 double norm_max = cut_out_norm.at(1).at(1).at(3).at(4);
357 double pow_max = cut_out_pow.at(1).at(1).at(3).at(4);
358 double bkg_max = cut_out_bkg.at(1).at(1).at(3).at(4);
359 B2DEBUG(30, "norm par max:" << norm_max);
360 B2DEBUG(30, "pow par max:" << pow_max);
361 B2DEBUG(30, "bkg par max:" << bkg_max);
362 B2DEBUG(30, "evaluate max cut:" << norm_max / (sqrt(sin(t_out_theta)) * pow(p_out_mom, pow_max)) + bkg_max);
363 B2DEBUG(30, "----------------------------------------");
364 }
365 }
366 //-----------end of debug lines ------//
367
368 if (c_validationON == 1) {
370 m_outputFile->cd();
371 for (int par = 0; par < c_nbinpar; par++) {
372 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
373 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
374 for (int theta = 0; theta < c_nbint; theta++) {
375 for (int p = 0; p < c_nbinp; p++) {
376 double layerdiff = lay2 - lay1;
377 if (layerdiff >= 0 && layerdiff < 3) {
378 if (m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->GetEntries() > 0) {
379 m_histo.at(par).at(lay1).at(lay2).at(theta).at(p)->Write();
380 }
381 }
382 }
383 }
384 }
385 }
386 }
387
389 for (int par = 0; par < c_nbinpar; par++) {
390 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
391 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
392 for (int minmax = 0; minmax < 2; minmax++) {
393 if (minmax == 0) {
394 int layerdiff = lay2 - lay1;
395 if (layerdiff >= 0 && layerdiff < 3) {
396 cut_m_histo.at(par).at(lay1).at(lay2)->Write();
397 }
398 }
399 if (minmax == 1) {
400 int layerdiff = lay2 - lay1;
401 if (layerdiff >= 0 && layerdiff < 3) {
402 cut_M_histo.at(par).at(lay1).at(lay2)->Write();
403 }
404 }
405 }
406 }
407 }
408 }
409 }
410
411 //--------------------------------OUTPUT: histogram booking, filling -------------------//
412
413 TH3F* output_norm_m = new TH3F("output_norm_m", "output_norm_m", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
414 TH3F* output_pow_m = new TH3F("output_pow_m", "output_pow_m", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
415 TH3F* output_bkg_m = new TH3F("output_bkg_m", "output_bkg_m", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
416
417 TH3F* output_norm_M = new TH3F("output_norm_M", "output_norm_M", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
418 TH3F* output_pow_M = new TH3F("output_pow_M", "output_pow_M", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
419 TH3F* output_bkg_M = new TH3F("output_bkg_M", "output_bkg_M", c_nbinpar, 0, 4, c_nbinlay, 0, 4, c_nbinlay, 0, 4);
420
421
422 for (int par = 0; par < c_nbinpar; par++) {
423 for (int lay1 = 0; lay1 < c_nbinlay; lay1++) {
424 for (int lay2 = 0; lay2 < c_nbinlay; lay2++) {
425 output_norm_m->SetBinContent(par, lay1, lay2, cut_out_norm.at(0).at(par).at(lay1).at(lay2));
426 output_norm_M->SetBinContent(par, lay1, lay2, cut_out_norm.at(1).at(par).at(lay1).at(lay2));
427
428 output_pow_m->SetBinContent(par, lay1, lay2, cut_out_pow.at(0).at(par).at(lay1).at(lay2));
429 output_pow_M->SetBinContent(par, lay1, lay2, cut_out_pow.at(1).at(par).at(lay1).at(lay2));
430
431 output_bkg_m->SetBinContent(par, lay1, lay2, cut_out_bkg.at(0).at(par).at(lay1).at(lay2));
432 output_bkg_M->SetBinContent(par, lay1, lay2, cut_out_bkg.at(1).at(par).at(lay1).at(lay2));
433 }
434 }
435 }
436
437 m_outputFile->cd();
438 output_norm_m->Write();
439 output_norm_M->Write();
440 output_pow_m->Write();
441 output_pow_M->Write();
442 output_bkg_m->Write();
443 output_bkg_M->Write();
444 m_outputFile->Close();
445 delete m_outputFile;
446
447 B2INFO("number of spacepoint with theta out of limits=" << m_tCounter);
448 B2INFO("number of spacepoint with momentum out of limits=" << m_pCounter);
449 B2INFO("number of tracks cutted by global cuts=" << m_globCounter);
450
451
452}
453
454void NoKickCutsEvalModule::terminate() {}
455
456
457double NoKickCutsEvalModule::deltaParEval(hitXP hit1, hitXP hit2, NoKickCuts::EParameters par, bool is0)
458{
459 double out = c_over;
460 int layer1 = hit1.m_sensorLayer;
461 int layer2 = hit2.m_sensorLayer;
462 double layerdiff = layer2 - layer1;
463 if (layerdiff >= 0 && layerdiff < 3) {
464 switch (par) {
465 case NoKickCuts::c_Omega:
466 out = abs(hit1.getOmegaEntry() - hit2.getOmegaEntry());
467 if (is0) out = abs(hit1.getOmega0() - hit2.getOmegaEntry());
468 break;
469
470 case NoKickCuts::c_D0:
471 out = hit1.getD0Entry() - hit2.getD0Entry();
472 if (is0) out = hit1.getD00() - hit2.getD0Entry();
473 break;
474
475 case NoKickCuts::c_Phi0:
476 out = asin(sin(hit1.getPhi0Entry())) - asin(sin(hit2.getPhi0Entry()));
477 if (is0) out = asin(sin(hit1.getPhi00())) - asin(sin(hit2.getPhi0Entry()));
478 break;
479
480 case NoKickCuts::c_Z0:
481 out = hit1.getZ0Entry() - hit2.getZ0Entry();
482 if (is0) out = hit1.getZ00() - hit2.getZ0Entry();
483 break;
484
485 case NoKickCuts::c_Tanlambda:
486 out = hit1.getTanLambdaEntry() - hit2.getTanLambdaEntry();
487 if (is0) out = hit1.getTanLambda0() - hit2.getTanLambdaEntry();
488 break;
489 }
490 }
491 return out;
492}
493
494double NoKickCutsEvalModule::cutFunction(int p, double pwidth)
495{
496 double out;
497 double mom = p * pwidth;
498 if (mom > 0.04)
499 out = -7.5 * pow(10, -7) / pow(mom, 3.88) + 1;
500 else out = 6.3 * mom + 0.57;
501 return out;
502}
Belle2::LogConfig::c_Debug
@ c_Debug
Debug: for code development.
Definition: LogConfig.h:26
Belle2::LogSystem::Instance
static LogSystem & Instance()
Static method to get a reference to the LogSystem instance.
Definition: LogSystem.cc:31
Belle2::Module
Base class for Modules.
Definition: Module.h:72
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Belle2::NoKickCutsEvalModule::c_pmax
const double c_pmax
maximum momentum evaluated
Definition: NoKickCutsEvalModule.h:92
Belle2::NoKickCutsEvalModule::m_outputFile
TFile * m_outputFile
output file of cuts
Definition: NoKickCutsEvalModule.h:111
Belle2::NoKickCutsEvalModule::cutFunction
double cutFunction(int p, double pwidth)
This is the funcion that select the percentage that has to be cut away from deltaPar distributions (f...
Definition: NoKickCutsEvalModule.cc:494
Belle2::NoKickCutsEvalModule::m_globCounter
int m_globCounter
counter of tracks cutted from global cuts
Definition: NoKickCutsEvalModule.h:106
Belle2::NoKickCutsEvalModule::c_tmax
const double c_tmax
150 degrees.
Definition: NoKickCutsEvalModule.h:94
Belle2::NoKickCutsEvalModule::c_multLimit
const double c_multLimit
multiplier of the range limit of the histograms of DeltaX
Definition: NoKickCutsEvalModule.h:102
Belle2::NoKickCutsEvalModule::initialize
void initialize() override
Initialize the Module.
Definition: NoKickCutsEvalModule.cc:44
Belle2::NoKickCutsEvalModule::c_validationON
bool c_validationON
flag to activate some validation plots
Definition: NoKickCutsEvalModule.h:107
Belle2::NoKickCutsEvalModule::deltaParEval
double deltaParEval(hitXP hit1, hitXP hit2, NoKickCuts::EParameters par, bool is0=false)
enum for the track-parameters
Definition: NoKickCutsEvalModule.cc:457
Belle2::NoKickCutsEvalModule::event
void event() override
This method is the core of the module.
Definition: NoKickCutsEvalModule.cc:99
Belle2::NoKickCutsEvalModule::m_histoLim
std::vector< double > m_histoLim
limits of DeltaX histograms
Definition: NoKickCutsEvalModule.h:112
Belle2::NoKickCutsEvalModule::m_namePar
std::vector< TString > m_namePar
name of track parameters
Definition: NoKickCutsEvalModule.h:116
Belle2::NoKickCutsEvalModule::endRun
void endRun() override
This method is called if the current run ends.
Definition: NoKickCutsEvalModule.cc:141
Belle2::NoKickCutsEvalModule::c_pmin
const double c_pmin
alternative cut function (not used, wider cuts)
Definition: NoKickCutsEvalModule.h:91
Belle2::NoKickCutsEvalModule::c_nbinp
const int c_nbinp
number of momentum bins
Definition: NoKickCutsEvalModule.h:96
Belle2::NoKickCutsEvalModule::terminate
void terminate() override
This method is called at the end of the event processing.
Definition: NoKickCutsEvalModule.cc:454
Belle2::NoKickCutsEvalModule::m_histo
std::vector< std::vector< std::vector< std::vector< std::vector< TH1F * > > > > > m_histo
DeltaX histograms.
Definition: NoKickCutsEvalModule.h:113
Belle2::NoKickCutsEvalModule::c_nbinlay
const int c_nbinlay
present IP too.
Definition: NoKickCutsEvalModule.h:98
Belle2::NoKickCutsEvalModule::c_over
const double c_over
escape flag of some methods
Definition: NoKickCutsEvalModule.h:103
Belle2::NoKickCutsEvalModule::c_nbinpar
const int c_nbinpar
number of track parameters
Definition: NoKickCutsEvalModule.h:97
Belle2::NoKickCutsEvalModule::c_nbint
const int c_nbint
number of theta parameters
Definition: NoKickCutsEvalModule.h:99
Belle2::NoKickCutsEvalModule::c_twidth
double c_twidth
width of theta bin
Definition: NoKickCutsEvalModule.h:101
Belle2::NoKickCutsEvalModule::c_fitMethod
bool c_fitMethod
flag to activate the fit method to evaluate the cuts
Definition: NoKickCutsEvalModule.h:108
Belle2::NoKickCutsEvalModule::c_tmin
const double c_tmin
17 degrees.
Definition: NoKickCutsEvalModule.h:93
Belle2::NoKickCutsEvalModule::m_pCounter
int m_pCounter
conter of hit out of range in momentum
Definition: NoKickCutsEvalModule.h:104
Belle2::NoKickCutsEvalModule::m_tCounter
int m_tCounter
counter of hit out of range in theta
Definition: NoKickCutsEvalModule.h:105
Belle2::NoKickCutsEvalModule::c_pwidth
double c_pwidth
width of momentum bin
Definition: NoKickCutsEvalModule.h:100
Belle2::NoKickCutsEvalModule::m_trackSel
NoKickRTSel m_trackSel
auxiliary variable to use methods of NoKickRTSel
Definition: NoKickCutsEvalModule.h:110
Belle2::NoKickCutsEvalModule::c_nbin
const int c_nbin
number of bins of histogram of DeltaX
Definition: NoKickCutsEvalModule.h:95
Belle2::NoKickCuts::EParameters
EParameters
enum for parameters name
Definition: NoKickCuts.h:44
Belle2::NoKickRTSel::m_8hitTrack
std::vector< hitXP > m_8hitTrack
vector of selected hit
Definition: NoKickRTSel.h:38
Belle2::NoKickRTSel::globalCut
bool globalCut(const std::vector< hitXP > &track8)
This method make some global cuts on the tracks (layer 3 and 6 required, d0 and z0 inside beam pipe).
Definition: NoKickRTSel.cc:111
Belle2::NoKickRTSel::m_setHitXP
std::set< hitXP, hitXP::timeCompare > m_setHitXP
set of hit to order the hit in time
Definition: NoKickRTSel.h:37
Belle2::NoKickRTSel::m_hitXP
std::vector< hitXP > m_hitXP
vector of hit, to convert the track
Definition: NoKickRTSel.h:36
Belle2::NoKickRTSel::hit8TrackBuilder
void hit8TrackBuilder(const RecoTrack &track)
this metod build a vector of hitXP from a track selecting the first hit on each layer of VXD (8 hit f...
Definition: NoKickRTSel.cc:96
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:79
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
Belle2::StoreAccessorBase::isRequired
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Definition: StoreAccessorBase.h:78
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
Belle2::hitXP
This class collects some information of a TrueHit, using SVDCLuster and MCParticle information too.
Definition: hitXP.h:32
Belle2::hitXP::getOmegaEntry
double getOmegaEntry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:325
Belle2::hitXP::getPhi00
double getPhi00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:367
Belle2::hitXP::getZ0Entry
double getZ0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:373
Belle2::hitXP::getOmega0
double getOmega0() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:331
Belle2::hitXP::m_sensorLayer
int m_sensorLayer
layer of the hit
Definition: hitXP.h:52
Belle2::hitXP::getD00
double getD00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:355
Belle2::hitXP::getPhi0Entry
double getPhi0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:361
Belle2::hitXP::getZ00
double getZ00() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:379
Belle2::hitXP::getTanLambda0
double getTanLambda0() const
evaluate relative parameter using IP position and momentum
Definition: hitXP.h:343
Belle2::hitXP::getD0Entry
double getD0Entry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:349
Belle2::hitXP::getTanLambdaEntry
double getTanLambdaEntry() const
evaluate relative parameter using entrypoint position and momentum
Definition: hitXP.h:337
Belle2::Module::addParam
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17