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