Belle II Software  release-08-01-10
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 
28 using namespace Belle2;
29 
30 REG_MODULE(NoKickCutsEval);
31 
32 NoKickCutsEvalModule::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 
44 void 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 
99 void 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 
141 void 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 
454 void NoKickCutsEvalModule::terminate() {}
455 
456 
457 double 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 
494 double 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:95
Belle2::NoKickCutsEvalModule::m_outputFile
TFile * m_outputFile
output file of cuts
Definition: NoKickCutsEvalModule.h:114
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:109
Belle2::NoKickCutsEvalModule::c_tmax
const double c_tmax
150 degrees.
Definition: NoKickCutsEvalModule.h:97
Belle2::NoKickCutsEvalModule::c_multLimit
const double c_multLimit
multiplier of the range limit of the histograms of DeltaX
Definition: NoKickCutsEvalModule.h:105
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:110
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:115
Belle2::NoKickCutsEvalModule::m_namePar
std::vector< TString > m_namePar
name of track parameters
Definition: NoKickCutsEvalModule.h:119
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:94
Belle2::NoKickCutsEvalModule::c_nbinp
const int c_nbinp
number of momentum bins
Definition: NoKickCutsEvalModule.h:99
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:116
Belle2::NoKickCutsEvalModule::c_nbinlay
const int c_nbinlay
present IP too.
Definition: NoKickCutsEvalModule.h:101
Belle2::NoKickCutsEvalModule::c_over
const double c_over
escape flag of some methods
Definition: NoKickCutsEvalModule.h:106
Belle2::NoKickCutsEvalModule::c_nbinpar
const int c_nbinpar
number of track parameters
Definition: NoKickCutsEvalModule.h:100
Belle2::NoKickCutsEvalModule::c_nbint
const int c_nbint
number of theta parameters
Definition: NoKickCutsEvalModule.h:102
Belle2::NoKickCutsEvalModule::c_twidth
double c_twidth
width of theta bin
Definition: NoKickCutsEvalModule.h:104
Belle2::NoKickCutsEvalModule::c_fitMethod
bool c_fitMethod
flag to activate the fit method to evaluate the cuts
Definition: NoKickCutsEvalModule.h:111
Belle2::NoKickCutsEvalModule::c_tmin
const double c_tmin
17 degrees.
Definition: NoKickCutsEvalModule.h:96
Belle2::NoKickCutsEvalModule::m_pCounter
int m_pCounter
conter of hit out of range in momentum
Definition: NoKickCutsEvalModule.h:107
Belle2::NoKickCutsEvalModule::m_tCounter
int m_tCounter
counter of hit out of range in theta
Definition: NoKickCutsEvalModule.h:108
Belle2::NoKickCutsEvalModule::c_pwidth
double c_pwidth
width of momentum bin
Definition: NoKickCutsEvalModule.h:103
Belle2::NoKickCutsEvalModule::m_trackSel
NoKickRTSel m_trackSel
auxiliary variable to use methods of NoKickRTSel
Definition: NoKickCutsEvalModule.h:113
Belle2::NoKickCutsEvalModule::c_nbin
const int c_nbin
number of bins of histogram of DeltaX
Definition: NoKickCutsEvalModule.h:98
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::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::REG_MODULE
REG_MODULE(arichBtest)
Register the Module.
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
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