Belle II Software  release-05-01-25
B2BIIFixMdstModule_ecl.cc
1 //
2 // $Id: B2BIIFixMdst_ecl.cc 11174 2010-09-29 06:29:44Z hitoshi $
3 //
4 // $Log$
5 //
6 // Revision 2.0 2015/03/11 tkeck
7 // Conversion to Belle II
8 //
9 // Revision 1.30 2005/07/28 09:56:02 hitoshi
10 // for exp43 run700-1149 (by Miyabayashi)
11 //
12 // Revision 1.29 2005/07/02 08:04:37 hitoshi
13 // for exp43 run300-699 (by Miyabayashi)
14 //
15 // Revision 1.28 2005/06/06 04:28:01 katayama
16 // New for exp. 43 run - 300 from Miyabayasi and Senyo san
17 //
18 // Revision 1.27 2005/04/17 09:38:42 hitoshi
19 // for Exp41 r1200-end (by Miyabayashi)
20 //
21 // Revision 1.26 2005/04/15 08:29:56 hitoshi
22 // for exp41 r800-1199 (by Miyabayashi)
23 //
24 // Revision 1.25 2005/03/29 16:24:45 hitoshi
25 // for Exp41 runno<=799 (by Miyabayashi)
26 //
27 // Revision 1.24 2005/03/10 15:09:52 hitoshi
28 // for Exp39 (by Miyabayashi san)
29 //
30 // Revision 1.23 2004/07/14 02:13:57 katayama
31 // for exp. 37-1913
32 //
33 // Revision 1.22 2004/07/10 00:24:04 katayama
34 // exp. 37 run < 1400
35 //
36 // Revision 1.21 2004/06/25 05:22:40 hitoshi
37 // *** empty log message ***
38 //
39 // Revision 1.20 2004/06/19 17:09:13 hitoshi
40 // for exp35 (by Miyabayashi san)
41 //
42 // Revision 1.19 2004/06/07 05:46:28 hitoshi
43 // for exp33. also mass window of make_pi0 is widened (by Miyabaya).
44 //
45 // Revision 1.18 2003/07/11 09:13:44 hitoshi
46 // added a comment (by Miyabayashi; no change in function itself).
47 //
48 // Revision 1.17 2003/06/30 22:16:29 hitoshi
49 // update for exp27 (by Miyabayashi, Sanjay, Senyo).
50 //
51 // Revision 1.16 2003/05/31 09:22:54 hitoshi
52 // updae for exp25 (ny Miyabayashi and Senyo).
53 //
54 // Revision 1.15 2003/03/19 05:02:53 hitoshi
55 // minor change (by Miyabayashi).
56 //
57 // Revision 1.14 2003/02/19 10:37:29 hitoshi
58 // incorporated an asymmetric pi0 resol function (by Miyabayashi).
59 //
60 // Revision 1.13 2002/12/28 23:10:08 katayama
61 // ecl and benergy for exp. 21 and 32 from Sanjay and Miyavayashi san
62 //
63 // Revision 1.12 2002/07/03 05:18:18 hitoshi
64 // ECL corrections for exp19 run900 - 1643 (by Miyabayashi san).
65 //
66 // Revision 1.11 2002/06/22 07:37:58 katayama
67 // From miyabayashi san. exp. 9 and exp. 19 ~899
68 //
69 // Revision 1.10 2002/06/19 12:10:36 katayama
70 // New for exp. 7 from Senyo and Miyabayashi sans
71 //
72 // Revision 1.9 2002/06/10 17:33:48 hitoshi
73 // added new corections for e11 (by Miyabayashi).
74 //
75 // Revision 1.8 2002/06/09 15:28:29 hitoshi
76 // added corrections for run<460 in e19.
77 //
78 // Revision 1.7 2002/06/01 12:09:32 hitoshi
79 // added corrections for runs1-299 in exp19 (by Miyabayashi).
80 //
81 // Revision 1.6 2002/05/24 08:03:30 hitoshi
82 // added corrections for exp13 reprocessed data (by Miyabayashi san).
83 //
84 // Revision 1.5 2002/05/09 08:57:55 hitoshi
85 // added corrections for exp17 (by Miyabayashi san).
86 //
87 // Revision 1.4 2002/04/23 04:57:46 katayama
88 // New correction from Miyabayashi san
89 //
90 // Revision 1.3 2002/04/19 07:24:39 katayama
91 // Added pi0/vertex, calibdedx
92 //
93 // Revision 1.2 2002/04/05 01:19:32 katayama
94 // templates for ecl_primary_vertex
95 //
96 // Revision 1.1 2002/03/13 02:55:20 katayama
97 // First version
98 //
99 //
100 #include <iostream>
101 #include <cmath>
102 #include <cfloat>
103 
104 #include <b2bii/modules/B2BIIMdstInput/B2BIIFixMdstModule.h>
105 #include "belle_legacy/panther/panther.h"
106 
107 #include "belle_legacy/tables/mdst.h"
108 #include "belle_legacy/tables/belletdf.h"
109 
110 #include "CLHEP/Vector/ThreeVector.h"
111 #include "CLHEP/Vector/LorentzVector.h"
112 #include "CLHEP/Matrix/Vector.h"
113 #include "CLHEP/Matrix/Matrix.h"
114 
115 namespace Belle2 {
122 //=====================================================================
123 //**** correct_ecl.cc ***** created 2001/June/11th *****
124 // Correct photon energy according to pi0 mass peak
125 // analysis.
126 // Version1.2.1 2002/Feb./20th Kenkichi Miya.
127 // Exp.15 proper correction curve and version control sequence
128 // were installed. option=0 and version=2 are recommended.
129 // Version1.2 2001/Nov./29th Kenkichi Miya.
130 // About exp.15, the correction curve of exp.13 will be applied
131 // for a while. 2001/Nov./29th.
132 // version control option added 2001/Nov./29th.
133 // The correction curve for exp.13 was implemented 2001/Nov./20th.
134 // Mdst_gamma and Mdst_ecl's energy and errors are corrected.
135 // input option was added 2001/July/23rd.
136 //******************************************************
137 
138 //==================================
139 // The function giving correction factor.
140 // Correcponding Data/MC so that energy in data should be divided by this.
141 //==================================
142 //Original definition.
143 //static double ecl_adhoc_corr( int Exp, int Run, double Energy, double)
144 //Modified 20081222
148  static double ecl_adhoc_corr(int Exp, int Run, int iflag05th,
149  double Energy, double)
150  {
151  double return_value;
152 
153  // Default.
154  return_value = 1.0;
155 
156  double x = std::log10(Energy * 1000.0);
157  switch (Exp) {
158  case 7:
159  //if( x < 2.5 )
160  //{
161  // return_value = 0.85445 + 0.89162e-1*x
162  // -0.96202e-2*x*x -0.10944e-2*x*x*x;
163  //}
164  // New curve for b20020416 processed mdst.
165  if (x < 2.6) {
166  return_value = 0.85878 + 0.81947e-1 * x - 0.99708e-2 * x * x
167  - 0.28161e-3 * x * x * x;
168  }
169  break;
170  case 9:
171  //if( x < 3.0 )
172  //{
173  // return_value = 0.83416 + 0.93534e-1*x
174  // -0.71632e-2*x*x -0.19072e-2*x*x*x;
175  //}
176  // New curve for b20020416 processed mdst.
177  if (x < 3.2) {
178  return_value = 0.86582 + 0.63194e-1 * x - 0.59391e-2 * x * x
179  + 0.17727e-2 * x * x * x - 0.62325e-3 * x * x * x * x;
180  }
181  break;
182  case 11:
183  // Proper curve for b20010507 processed mdst.
184  //if( x < 3.0 )
185  //{
186  // return_value = 0.84927 + 0.86308e-1*x
187  // -0.93898e-2*x*x -0.10557e-2*x*x*x;
188  //}
189  //
190  // Proper curve for b20020416 processed mdst.
191  if (x < 3.7) {
192  if (x < 3.0) {
193  return_value = 0.88989 + 0.41082e-1 * x
194  - 0.21919e-2 * x * x + 0.27116e-2 * x * x * x - 0.89113e-3 * x * x * x * x;
195  } else { // Linear extrapolation to unity.
196  return_value = 0.994441 + 4.90176e-3 * (x - 3.0);
197  }
198  }
199  break;
200  case 13:
201  // Proper curve for b20010507 processed mdst.
202  //if( x < 3.0 )
203  //{
204  //return_value = 0.78138 + 0.17899*x
205  // -0.48821e-1*x*x +0.43692e-2*x*x*x;
206  //}
207  //
208  // Proper curve for b20020416 processed mdst.
209  if (x < 3.2) {
210  return_value = 0.87434 + 0.62145e-1 * x
211  - 0.29691e-2 * x * x - 0.15843e-2 * x * x * x + 0.86858e-4 * x * x * x * x;
212  }
213  break;
214  case 15:
215  // Proper curve for b20020125 processed mdst.
216  //if( x < 3.0 )
217  //{
218  // return_value = 0.86858 + 0.75587e-1*x
219  // -0.11150e-1*x*x -0.21454e-4*x*x*x;
220  //}
221  //else if( x < 3.67 ) // Lenear extrapolation to be unity.
222  //{
223  // return_value = 0.994412 + 8.11e-3*( x - 3.0 );
224  //}
225  // Proper curve for b20020416 processed mdst.
226  if (x < 3.2) {
227  return_value = 0.83073 + 0.10137 * x
228  - 0.59946e-2 * x * x - 0.56516e-2 * x * x * x + 0.87225e-3 * x * x * x * x;
229  }
230 
231  break;
232  case 17:
233  // Proper curve for b20020416 processed mdst.
234  if (x < 3.0) {
235  return_value = 0.89260 + 0.54731e-1 * x
236  - 0.25736e-2 * x * x - 0.16493e-2 * x * x * x + 0.1032e-3 * x * x * x * x;
237  } else if (x < 3.5) {
238  return_value = 0.997459 + 0.0059039 * (x - 3.0);
239  }
240  break;
241  case 19:
242  // Correction curve obtained by Run001-299 HadronB
243  // processed by b20020416.
244  if (Run < 300) {
245  if (x < 3.0) {
246  return_value = 0.86432 + 0.87554e-1 * x
247  - 0.84182e-2 * x * x - 0.39880e-2 * x * x * x + 0.69435e-3 * x * x * x * x;
248  }
249  }
250  if (299 < Run && Run < 900) { // For Run300-899 or later by b20020416.
251  if (x < 2.9) {
252  return_value = 0.92082 + 0.45896e-1 * x
253  - 0.68067e-2 * x * x + 0.94055e-3 * x * x * x - 0.27717e-3 * x * x * x * x;
254  }
255  }
256  if (899 < Run) { // Till the end.
257  if (x < 3.1) {
258  return_value = 0.85154 + 0.97812e-1 * x
259  - 0.85774e-2 * x * x - 0.59092e-2 * x * x * x + 0.1121e-2 * x * x * x * x;
260  }
261  }
262  break;
263  case 21:
264  // Proper curve for b20021205 processed mdst.
265  if (x < 3.0) {
266  return_value = 0.84940 + 0.86266e-1 * x
267  - 0.82204e-2 * x * x - 0.12912e-2 * x * x * x;
268  }
269  break;
270  case 23:
271  // Proper curve for b20021205 processed mdst.
272  if (x < 3.1) {
273  return_value = 0.85049 + 0.85418e-1 * x
274  - 0.45747e-2 * x * x - 0.40081e-2 * x * x * x + 0.52113e-3 * x * x * x * x;
275  }
276  break;
277  case 25:
278  if (x < 3.0) {
279  return_value = 0.87001 + 0.73693e-1 * x
280  - 0.80094e-2 * x * x - 0.78527e-3 * x * x * x + 0.25888e-4 * x * x * x * x;
281  }
282  break;
283  case 27: // It turned out no signif. change through all the runs.
284  if (x < 2.85) {
285  return_value = 0.90051 + 0.56094e-1 * x
286  - 0.14842e-1 * x * x + 0.55555e-2 * x * x * x - 0.10378e-2 * x * x * x * x;
287  }
288  break;
289  case 31:
290  if (iflag05th == 1) { // For 0.5MeV threshold.
291  // For b20040429.
292  if (x < 2.95) {
293  return_value = 0.80295 + 0.13395 * x
294  - 0.10773e-1 * x * x - 0.85861e-2 * x * x * x + 0.15331e-2 * x * x * x * x;
295  }
296  } else { // For theta-dep. threshold.
297  if (x < 3.25) {
298  return_value = 0.57169 + 0.32548 * x
299  - 0.41157e-1 * x * x - 0.21971e-1 * x * x * x + 0.50114e-2 * x * x * x * x;
300  }
301  }
302  break;
303  case 33:
304  if (iflag05th == 1) { // For 0.5MeV threshold.
305  // For b20040429.
306  if (x < 2.64) {
307  return_value = 0.81153 + 0.10847 * x
308  - 0.24652e-2 * x * x - 0.54738e-2 * x * x * x + 0.41243e-3 * x * x * x * x;
309  }
310  } else { // For theta-dep. threshold, b20090127.
311  if (x < 3.15) {
312  return_value = 0.59869 + 0.29276 * x
313  - 0.15849e-1 * x * x - 0.31322e-1 * x * x * x + 0.62491e-2 * x * x * x * x;
314  }
315  }
316  break;
317  case 35:
318  if (iflag05th == 1) { // For 0.5MeV threshold.
319  // For b20040429.
320  if (x < 2.54) {
321  return_value = 0.83528 + 0.10402 * x
322  - 0.62047e-2 * x * x - 0.62411e-2 * x * x * x + 0.10312e-2 * x * x * x * x;
323  }
324  } else { // For theta-dep. threshold, b20090127.
325  if (x < 3.2) {
326  return_value = 0.58155 + 0.30642 * x
327  - 0.16981e-1 * x * x - 0.32609e-1 * x * x * x + 0.64874e-2 * x * x * x * x;
328  }
329  }
330  break;
331  case 37:
332  if (iflag05th == 1) { // For 0.5MeV threshold.
333  if (x < 2.5) {
334  return_value = 0.83706 + 0.10726 * x
335  - 0.62423e-2 * x * x - 0.42425e-2 * x * x * x;
336  }
337  } else { // For theta-dep. threshold.
338  if (x < 3.15) {
339  return_value = 0.58801 + 0.30569 * x
340  - 0.18832e-1 * x * x - 0.32116e-1 * x * x * x + 0.64899e-2 * x * x * x * x;
341  }
342  }
343  break;
344  case 39:
345  if (iflag05th == 1) { // For 0.5MeV threshold.
346  if (x < 3.0) {
347  if (x < 2.8) {
348  return_value = 0.80827 + 0.095353 * x
349  - 0.14818e-2 * x * x - 0.23854e-2 * x * x * x - 0.22454e-3 * x * x * x * x;
350  } else {
351  return_value = 0.0112468 * (x - 2.8) + 0.997847;
352  }
353  }
354  } else { // For theta-dep. threshold.
355  if (x < 3.25) {
356  return_value = 0.61133 + 0.27115 * x
357  - 0.12724e-1 * x * x - 0.28167e-1 * x * x * x + 0.54699e-2 * x * x * x * x;
358  }
359  }
360  break;
361  case 41:
362  if (iflag05th == 1) { // For 0.5MeV threshold.
363  if (Run < 800) {
364  if (x < 2.5) {
365  return_value = 0.90188 + 0.76563e-1 * x
366  - 0.16328e-1 * x * x + 0.56816e-3 * x * x * x;
367  }
368  }
369  if (799 < Run) { // Run1200-1261 no signif. change from Run800-1199.
370  if (x < 2.95) {
371  return_value = 0.88077 + 0.71720e-1 * x
372  - 0.92197e-2 * x * x - 0.15932e-2 * x * x * x + 0.38133e-3 * x * x * x * x;
373  }
374  }
375  } else { // For theta-dep. threshold.
376  if (x < 3.2) {
377  return_value = 0.57808 + 0.31083 * x
378  - 0.20247e-1 * x * x - 0.31658e-1 * x * x * x + 0.64087e-2 * x * x * x * x;
379  }
380  }
381  break;
382  case 43: // modified 20090829.
383  if (iflag05th == 1) { // For 0.5MeV threshold.
384  if (Run < 300) {
385  if (x < 3.3) {
386  return_value = 0.85592 + 0.93352e-1 * x
387  - 0.93144e-2 * x * x - 0.43681e-2 * x * x * x + 0.7971e-3 * x * x * x * x;
388  }
389  }
390  if (299 < Run && Run < 700) { // added 20050630.
391  if (x < 3.0) {
392  return_value = 0.89169 + 0.73301e-1 * x
393  - 0.13856e-1 * x * x + 0.47303e-3 * x * x * x;
394  }
395  }
396  if (699 < Run) { // added 20050727
397  if (x < 3.1) {
398  return_value = 0.90799 + 0.69815e-1 * x
399  - 0.17490e-1 * x * x + 0.14651e-2 * x * x * x;
400  }
401  }
402  } else { // For theta-dep. threshold b20090127.
403  if (x < 3.25) {
404  return_value = 0.58176 + 0.30245 * x - 0.16390e-1 * x * x
405  - 0.32258e-1 * x * x * x + 0.64121e-2 * x * x * x * x;
406  }
407  //if( 999 < Run ){// For 5S runs, 2/fb, obsolete 20091209.
408  //if( x < 3.3 ){
409  //return_value = 0.58463 + 0.29780*x -0.14441e-1*x*x
410  //-0.32323e-1*x*x*x + 0.63607e-2*x*x*x*x;
411  //}
412  //}
413  }
414  break;
415  case 45:
416  if (iflag05th == 1) { // For 0.5MeV threshold.
417  if (x < 3.1) { // added 20060413
418  return_value = 0.84823 + 0.10394 * x
419  - 0.92233e-2 * x * x - 0.72700e-2 * x * x * x + 0.14552e-2 * x * x * x * x;
420  }
421  } else { // For theta-dep. threshold b20090127.
422  if (x < 3.25) {
423  return_value = 0.62797 + 0.23897 * x
424  + 0.10261e-1 * x * x - 0.35700e-1 * x * x * x + 0.63846e-2 * x * x * x * x;
425  }
426  }
427  break;
428  case 47:
429  if (iflag05th == 1) { // For 0.5MeV threshold.
430  if (x < 3.1) { // added 20060420
431  return_value = 0.86321 + 0.82987e-1 * x
432  - 0.12139e-1 * x * x - 0.12920e-3 * x * x * x;
433  }
434  } else { // For theta-dep. threshold b20090127.
435  if (x < 3.25) {
436  return_value = 0.58061 + 0.30399 * x
437  - 0.17520e-1 * x * x - 0.31559e-1 * x * x * x + 0.62681e-2 * x * x * x * x;
438  }
439  }
440  break;
441  case 49:
442  if (iflag05th == 1) { // For 0.5MeV threshold.
443  if (x < 3.13) { // added 20060511
444  if (x < 2.80) {
445  return_value = 0.85616 + 0.78041e-1 * x
446  - 0.38987e-2 * x * x - 0.31224e-2 * x * x * x + 0.33374e-3 * x * x * x * x;
447  } else {
448  return_value = ((1.0 - 0.99608) / (3.13 - 2.80)) * (x - 2.80) + 0.99608;
449  }
450  }
451  } else { // For theta-dep. threshold b20090127.
452  if (x < 3.3) {
453  return_value = 0.57687 + 0.29948 * x
454  - 0.10594e-1 * x * x - 0.34561e-1 * x * x * x + 0.67064e-2 * x * x * x * x;
455  }
456  }
457  break;
458  case 51:
459  if (iflag05th == 1) { // For 0.5MeV threshold, added 20070323.
460  if (x < 3.07) {
461  return_value = 0.89063 + 0.72152e-1 * x
462  - 0.14143e-1 * x * x + 0.73116e-3 * x * x * x;
463  }
464  } else { // For theta-dep. threshold b20090127.
465  if (x < 3.3) {
466  return_value = 0.59310 + 0.28618 * x
467  - 0.13858e-1 * x * x - 0.30230e-1 * x * x * x + 0.59173e-2 * x * x * x * x;
468  }
469  }
470  break;
471  case 53: // modified 20090807.
472  if (iflag05th == 1) { // For 0.5MeV threshold.
473  if (x < 2.63) {
474  //std::cout<<"Exp.53 0.5MeV threshold."<<std::endl;
475  return_value = 0.83425 + 0.10468 * x
476  - 0.53641e-2 * x * x - 0.60276e-2 * x * x * x + 0.77763e-3 * x * x * x * x;
477  }
478  } else { // For theta-dep. threshold b20090127.
479  if (x < 3.1) {
480  //std::cout<<"Exp.53 theta-dep. threshold."<<std::endl;
481  return_value = 0.66100 + 0.25192 * x
482  - 0.16419e-1 * x * x - 0.25720e-1 * x * x * x + 0.52189e-2 * x * x * x * x;
483  }
484  }
485  break;
486  case 55:
487  if (iflag05th == 1) { // For 0.5MeV threshold.
488  if (x < 3.3) {
489  return_value = 0.86202 + 0.79575e-1 * x
490  - 0.66721e-2 * x * x - 0.28609e-2 * x * x * x + 0.42784e-3 * x * x * x * x;
491  }
492  } else { // For theta-dep. threshold b20090127.
493  if (x < 3.25) {
494  return_value = 0.58789 + 0.29310 * x
495  - 0.15784e-1 * x * x - 0.30619e-1 * x * x * x + 0.60648e-2 * x * x * x * x;
496  }
497  }
498  break;
499  case 61:
500  if (Run < 1210) { // On 4S by Th.-dep. threshold.
501  if (x < 3.5) {
502  return_value = 0.68839 + 0.18218 * x
503  - 0.23140e-2 * x * x - 0.17439e-1 * x * x * x + 0.29960e-2 * x * x * x * x;
504  }
505  }
506  if (Run > 1209) { // 5S scan.
507  if (iflag05th == 1) { // processed by 0.5MeV threshold.
508  if (x < 3.3) {
509  if (x < 2.8) {
510  return_value = 0.94294 - 0.77497e-2 * x
511  - 0.43464e-3 * x * x + 0.99837e-2 * x * x * x - 0.23726e-2 * x * x * x * x;
512  } else {
513  return_value = 0.0162997 * (x - 2.80) + 0.991162;
514  }
515  }
516  } else { // 5S scan. by Th.-dep. threshold.
517  if (x < 3.5) {
518  return_value = 0.64816 + 0.22492 * x
519  - 0.91745e-2 * x * x - 0.21736e-1 * x * x * x + 0.41333e-2 * x * x * x * x;
520  }
521  }
522  }
523  break;
524  case 63:
525  if (x < 3.4) {
526  return_value = 0.69302 + 0.18393 * x
527  - 0.10983e-1 * x * x - 0.14148e-1 * x * x * x + 0.27298e-2 * x * x * x * x;
528  }
529  break;
530  case 65:
531  if (Run > 1000) { // 1S run turned out to be same as 5S scan.
532  if (iflag05th == 1) { // For 0.5MeV threshold.
533  if (x < 3.3) {
534  if (x < 2.8) {
535  return_value = 0.94294 - 0.77497e-2 * x
536  - 0.43464e-3 * x * x + 0.99837e-2 * x * x * x - 0.23726e-2 * x * x * x * x;
537  } else {
538  return_value = 0.0162997 * (x - 2.80) + 0.991162;
539  }
540  }
541  } else { // For theta-dep. threshold, b20090127.
542  if (x < 3.6) {
543  return_value = 0.70987 + 0.16230 * x
544  - 0.64867e-2 * x * x - 0.12021e-1 * x * x * x + 0.20874e-2 * x * x * x * x;
545  }
546  }
547  } else { // 4S runs processed with new lib.
548  if (x < 3.4) {
549  return_value = 0.66833 + 0.20602 * x
550  - 0.14322e-1 * x * x - 0.15712e-1 * x * x * x + 0.31114e-2 * x * x * x * x;
551  }
552  }
553  break;
554  case 67:
555  if (x < 3.4) {
556  return_value = 0.64196 + 0.23069 * x
557  - 0.20303e-1 * x * x - 0.15680e-1 * x * x * x + 0.31611e-2 * x * x * x * x;
558  }
559  break;
560  case 69:
561  if (x < 3.35) {
562  //return_value = 0.64196 + 0.23752*x
563  //-0.85197e-2*x*x -0.24366e-1*x*x*x +0.46397e-2*x*x*x*x;
564  return_value = 0.99 * (0.64196 + 0.23752 * x
565  - 0.85197e-2 * x * x - 0.24366e-1 * x * x * x
566  + 0.46397e-2 * x * x * x * x);
567  }
568  break;
569  case 71:
570  if (x < 3.35) {
571  return_value = 0.66541 + 0.12579 * x
572  + 0.89999e-1 * x * x - 0.58305e-1 * x * x * x + 0.86969e-2 * x * x * x * x;
573  }
574  break;
575  case 73:
576  if (x < 3.45) {
577  return_value = 0.62368 + 0.24142 * x
578  - 0.11677e-1 * x * x - 0.22477e-1 * x * x * x + 0.42765e-2 * x * x * x * x;
579  }
580  break;
581  }
582  // This return value has to be used as division.
583  // i.e. E_corrected = E_original/return_vale.
584  return return_value;
585  }
586 
587 //=================================================================
588 // Fix up wrong calib. in Exp.45 by Isamu Nakamura, added 20060413.
589 //=================================================================
591  static double ecl_adhoc_corr_45(int exp, int /*run*/, int cid)
592  {
593 
594  int bcid[] = {
595  243, 244, 316, 317, 318, 412, 413, 414, 508, 509, 510, 604, 605, 606
596  };
597 
598  double bcfact[] = {
599  1.06215 , 1.06589 , 0.953827, 1.06607 , 0.943696, 0.910855, 1.01201 ,
600  1.01704 , 0.954612, 1.02761 , 1.00716 , 1.0319 , 1.03135 , 1.05782
601  };
602 
603  if (exp != 45) return 1.;
604 
605  for (int i = 0 ; i < 14 ; i++)
606  if (bcid[i] == cid)
607  return bcfact[i];
608 
609  return 1.;
610  }
611 
612 //=====================
613 // Correct energy scale (MC) to make pi0 peak nominal. 20080530
614 //=====================
616  static double ecl_mcx3_corr(int /*exp*/, int /*run*/, double energy, double /*theta*/)
617  {
618  //KM std::cout<<"ecl_mcx3_corr called."<<std::endl;
619  double return_value;
620  double x = std::log10(energy * 1000.0);
621  return_value = 1.0;
622  if (x < 3.2) {
623  return_value = 0.72708 + 0.22735 * x
624  - 0.20603e-1 * x * x - 0.24644e-1 * x * x * x + 0.53480e-2 * x * x * x * x;
625  }
626  //std::cout<<"energy="<<energy<<" ecl_mcx3="<<return_value<<std::endl;
627  return return_value;
628  }
629 
630 //===== mpi0pdg.cc ===== cerated 2001/07/17 =====
631 // Make MC mass peak to PDG's one.
632 // Version0.1 2001/07/17 trial version.
633 //===============================================
635  static double mpi0pdg(double Energy) // Energy in GeV.
636  {
637  // return value.
638  double return_value;
639  // pi0 mass in PDG.
640  const double m_pdg = 134.9766;
641 
642  // default = do nothing.
643  return_value = 1.0;
644 
645  if (0.1 < Energy && Energy < 1.0) {
646  double x = log10(Energy * 1000.0);
647  // Curve given by S.Uchida at 2001/07/17.
648  return_value
649  = m_pdg / (135.21 + 5.3812 * x - 4.2160 * x * x + 0.74650 * x * x * x);
650  }
651  // This return value has to be used as product.
652  // i.e. E_corrected = E_original*return_value.
653  return return_value;
654  }
655 
656 //**** pi0resol ***** implemented 2003/Feb./17th *****
657 // Standalone func. of pi0 mass resolution function
658 // coded by Hayashii-san.
659 //*******************************
660 //*** Original comments ********
661 // pi0resol: pi0 resolution function
662 //
663 //******************************
664 // input p : pi0 lab momentum in GeV
665 // theta : pi0 lab polar angle in the unit of Deg.
666 // side :
667 // ="lower" ; return the sigma of the lower(left) side.
668 // ="higher" ; return the sigma of the higher(right) side.
669 // mcdata:
670 // =0 accessing Exp. data
671 // =1 accessing MC
672 // exp : experiment number
673 //
674 // output : return pi0 mass resolution in the unit of GeV.
675 //**********************************
676 // Original definition.
677 //double pi0resol(double p, double theta, char* side, bool mcdata, int exp )
678  double B2BIIFixMdstModule::pi0resol(double p, double theta, const char* side, bool mcdata,
679  int /*exp*/, int /*option*/)
680  {
681 // option is added to become compatible with future gsim modification.
682  int iside = 0;
683 
684  if (!strcmp(side, "lower")) iside = 1;
685  if (!strcmp(side, "higher")) iside = 2;
686  if (iside == 0) {
687  B2ERROR("Error pi0resol. Wrong input parameter=" << side);
688  return -1;
689  }
690 // dout(Debugout::INFO,"B2BIIFixMdst_ecl") <<" iside="<<iside<<std::endl;
691 //----------------------------------------------------
692 // 20 <= theta <= 45
693 //----------------------------------------------------
694 // 0.15 < p < 0.8
695 // -- data L: lower side H: higher side
696  const double data_thf_pl_L[] = { 9.8233, -13.3489, 18.3271, -7.6668};
697  const double derr_thf_pl_L[] = { 0.3274, 0.755, 0.8611, 0.521};
698  const double data_thf_pl_H[] = { 6.1436, -7.9668, 12.4766, -5.3201};
699  const double derr_thf_pl_H[] = { 0.2903, 0.6754, 0.7724, 0.466};
700 // 0.8 < p < 5.3
701 // -- data L: lower side H: higher side
702  const double data_thf_ph_L[] = { 7.1936, -0.6378, 0.5912, -0.075};
703  const double derr_thf_ph_L[] = { 0.1975, 0.2036, 0.0793, 0.0121};
704  const double data_thf_ph_H[] = { 4.4923, 0.5532, 0.2658, -0.0343};
705  const double derr_thf_ph_H[] = { 0.1883, 0.2019, 0.0803, 0.0122};
706 //----------------------------------------------------
707 // 20 <= theta <= 45
708 //----------------------------------------------------
709 // 0.15 < p < 0.8
710 // -- MC L: lower side H: higher side
711  const double mc_thf_pl_L[] = { 4.8093, 3.4567, -3.7898, 1.7553};
712  const double merr_thf_pl_L[] = { 0.2145, 0.527, 0.6153, 0.3712};
713  const double mc_thf_pl_H[] = { 4.6176, -2.9049, 4.2994, -1.4776};
714  const double merr_thf_pl_H[] = { 0.1969, 0.4826, 0.555, 0.3346};
715 // 0.8 < p < 5.3
716 // -- MC L: lower side H: higher side
717  const double mc_thf_ph_L[] = { 6.3166, -0.5993, 0.5845, -0.0695};
718  const double merr_thf_ph_L[] = { 0.1444, 0.1468, 0.0571, 0.0087};
719  const double mc_thf_ph_H[] = { 2.9719, 1.7999, -0.2418, 0.028};
720  const double merr_thf_ph_H[] = { 0.1318, 0.1362, 0.0538, 0.0082};
721 
722 //----------------------------------------------------
723 // 45 <= theta <= 100
724 //----------------------------------------------------
725 // 0.15 < p < 0.8
726 // -- data L: lower side H: higher side
727  const double data_thm_pl_L[] = { 7.7573, -4.8855, 6.5561, -2.4788};
728  const double derr_thm_pl_L[] = { 0.1621, 0.4894, 0.6882, 0.4985};
729  const double data_thm_pl_H[] = { 6.9075, -10.5036, 18.5196, -9.224};
730  const double derr_thm_pl_H[] = { 0.1458, 0.4383, 0.639, 0.4726};
731 // 0.8 < p < 5.3
732 // -- data L: lower side H: higher side
733  const double data_thm_ph_L[] = { 5.2347, 2.1827, -0.6563, 0.1545};
734  const double derr_thm_ph_L[] = { 0.0986, 0.1281, 0.0627, 0.0117};
735  const double data_thm_ph_H[] = { 3.2114, 3.3806, -0.8635, 0.1371};
736  const double derr_thm_ph_H[] = { 0.0927, 0.1205, 0.058, 0.0106};
737 //
738 //----------------------------------------------------
739 // 45 <= theta <= 100
740 //----------------------------------------------------
741 // 0.15 < p < 0.8
742 // -- MC L: lower side H: higher side
743  const double mc_thm_pl_L[] = { 6.1774, -2.1831, 3.6615, -1.1813};
744  const double merr_thm_pl_L[] = { 0.11, 0.327, 0.476, 0.3655};
745  const double mc_thm_pl_H[] = { 4.0239, -0.7485, 3.6203, -1.4823};
746  const double merr_thm_pl_H[] = { 0.0991, 0.3034, 0.4223, 0.3069};
747 // 0.8 < p < 5.3
748 // -- MC L: lower side H: higher side
749  const double mc_thm_ph_L[] = { 4.5966, 2.261, -0.4938, 0.0984};
750  const double merr_thm_ph_L[] = { 0.0711, 0.0917, 0.0448, 0.0081};
751  const double mc_thm_ph_H[] = { 3.4609, 2.0069, -0.0498, -0.0018};
752  const double merr_thm_ph_H[] = { 0.0966, 0.1314, 0.0574, 0.0086};
753 //
754 //----------------------------------------------------
755 // 100 <= theta <= 150
756 //----------------------------------------------------
757 // 0.15 < p < 0.8
758 // -- data L: lower side H: higher side
759  const double data_thb_pl_L[] = {11.5829, -21.6715, 30.2368, -13.0389};
760  const double derr_thb_pl_L[] = { 0.2742, 0.7256, 0.9139, 0.6006};
761  const double data_thb_pl_H[] = { 8.0227, -14.7387, 23.1042, -10.4233};
762  const double derr_thb_pl_H[] = { 0.2466, 0.6512, 0.8239, 0.5419};
763 // 0.8 < p < 5.3
764 // -- data L: lower side H: higher side
765  const double data_thb_ph_L[] = { 7.5872, -1.8994, 1.6526, -0.2755};
766  const double derr_thb_ph_L[] = { 0.1999, 0.2638, 0.1422, 0.0335};
767  const double data_thb_ph_H[] = { 6.3542, -2.5164, 2.5763, -0.4803};
768  const double derr_thb_ph_H[] = { 0.1885, 0.2527, 0.136, 0.0318};
769 //----------------------------------------------------
770 // 100 <= theta <= 150
771 //----------------------------------------------------
772 // 0.15 < p < 0.8
773 // -- MC L: lower side H: higher side
774  const double mc_thb_pl_L[] = { 5.2707, 2.5607, -3.1377, 1.8434};
775  const double merr_thb_pl_L[] = { 0.1801, 0.5048, 0.6741, 0.4343};
776 
777  const double mc_thb_pl_H[] = { 2.5867, 7.6982, -10.0677, 5.312};
778  const double merr_thb_pl_H[] = { 0.1658, 0.4651, 0.6063, 0.3925};
779 // 0.8 < p < 5.3
780 // -- MC L: lower side H: higher side
781  const double mc_thb_ph_L[] = { 6.5206, -1.5103, 1.9054, -0.3609};
782  const double merr_thb_ph_L[] = { 0.1521, 0.2048, 0.1108, 0.0255};
783  const double mc_thb_ph_H[] = { 3.4397, 2.2372, -0.1214, -0.0004};
784  const double merr_thb_ph_H[] = { 0.1324, 0.1822, 0.1065, 0.0252};
785 
786 //
787  double resol;
788  double para[4] = { };
789  double error[4] = { };
790 
791  double pbuf = p;
792 //--
793 // theta< 45
794 //---
795  if (theta <= 45.) {
796 //--
797 // p-low
798 //---
799  if (pbuf <= 0.8) {
800  if (!mcdata) {
801  //-- data
802  for (int i = 0; i <= 3; i++) {
803  if (iside == 1) {
804  para[i] = data_thf_pl_L[i];
805  error[i] = derr_thf_pl_L[i];
806  } else {
807  para[i] = data_thf_pl_H[i];
808  error[i] = derr_thf_pl_H[i];
809  }
810  } // for
811  } else {
812  //-- mc
813  for (int i = 0; i <= 3; i++) {
814  if (iside == 1) {
815  para[i] = mc_thf_pl_L[i];
816  error[i] = merr_thf_pl_L[i];
817  } else {
818  para[i] = mc_thf_pl_H[i];
819  error[i] = merr_thf_pl_H[i];
820  }
821  } //for
822  }
823  } else {
824 //--
825 // p-high
826 //---
827 // use p= 5.2 value if p >= 5.2 for forward region
828  if (pbuf >= 5.2) {pbuf = 5.2;}
829  if (!mcdata) {
830  //-- data
831  for (int i = 0; i <= 3; i++) {
832  if (iside == 1) {
833  para[i] = data_thf_ph_L[i];
834  error[i] = derr_thf_ph_L[i];
835  } else {
836  para[i] = data_thf_ph_H[i];
837  error[i] = derr_thf_ph_H[i];
838  }
839  }// for
840  } else {
841  //-- mc
842  for (int i = 0; i <= 3; i++) {
843  if (iside == 1) {
844  para[i] = mc_thf_ph_L[i];
845  error[i] = merr_thf_ph_L[i];
846  } else {
847  para[i] = mc_thf_ph_H[i];
848  error[i] = merr_thf_ph_H[i];
849  }
850  }// for
851  }
852  } // p-range
853 //--
854 // 45< theta< 100
855 //---
856  } else if (theta > 45. && theta <= 100.) {
857 //--
858 // p-low
859 //---
860  if (pbuf <= 0.8) {
861  if (!mcdata) {
862  //-- data
863  for (int i = 0; i <= 3; i++) {
864  if (iside == 1) {
865  para[i] = data_thm_pl_L[i];
866  error[i] = derr_thm_pl_L[i];
867  } else {
868  para[i] = data_thm_pl_H[i];
869  error[i] = derr_thm_pl_H[i];
870  }
871  } // for
872  } else {
873  //-- mc
874  for (int i = 0; i <= 3; i++) {
875  if (iside == 1) {
876  para[i] = mc_thm_pl_L[i];
877  error[i] = merr_thm_pl_L[i];
878  } else {
879  para[i] = mc_thm_pl_H[i];
880  error[i] = merr_thm_pl_H[i];
881  }
882  } //for
883  }
884  } else {
885 //--
886 // p-high
887 //---
888 // use p= 5.2 value if p >= 5.2 for middle region
889  if (pbuf >= 5.2) {pbuf = 5.2;}
890  if (!mcdata) {
891  //-- data
892  for (int i = 0; i <= 3; i++) {
893  if (iside == 1) {
894  para[i] = data_thm_ph_L[i];
895  error[i] = derr_thm_ph_L[i];
896  } else {
897  para[i] = data_thm_ph_H[i];
898  error[i] = derr_thm_ph_H[i];
899  }
900  } //for
901  } else {
902  //-- mc
903  for (int i = 0; i <= 3; i++) {
904  if (iside == 1) {
905  para[i] = mc_thm_ph_L[i];
906  error[i] = merr_thm_ph_L[i];
907  } else {
908  para[i] = mc_thm_ph_H[i];
909  error[i] = merr_thm_ph_H[i];
910  }
911  } //for
912  }
913  } // p-range
914 //--
915 // theta> 100
916 //---
917  } else if (theta > 100.) {
918 //--
919 // p-low
920 //---
921  if (pbuf <= 0.8) {
922  if (!mcdata) {
923  //-- data
924  for (int i = 0; i <= 3; i++) {
925  if (iside == 1) {
926  para[i] = data_thb_pl_L[i];
927  error[i] = derr_thb_pl_L[i];
928  } else {
929  para[i] = data_thb_pl_H[i];
930  error[i] = derr_thb_pl_H[i];
931  }
932  } //for
933  } else {
934  //-- mc
935  for (int i = 0; i <= 3; i++) {
936  if (iside == 1) {
937  para[i] = mc_thb_pl_L[i];
938  error[i] = merr_thb_pl_L[i];
939  } else {
940  para[i] = mc_thb_pl_H[i];
941  error[i] = merr_thb_pl_H[i];
942  }
943  } //for
944  }
945  } else {
946 //--
947 // p-high
948 //---
949 // use p= 3.0 value if p >= 3.0 for backward region
950  if (pbuf >= 3.0) {pbuf = 3.0;}
951  if (!mcdata) {
952  //-- data
953  for (int i = 0; i <= 3; i++) {
954  if (iside == 1) {
955  para[i] = data_thb_ph_L[i];
956  error[i] = derr_thb_ph_L[i];
957  } else {
958  para[i] = data_thb_ph_H[i];
959  error[i] = derr_thb_ph_H[i];
960  }
961  } //for
962  } else {
963  //-- mc
964  for (int i = 0; i <= 3; i++) {
965  if (iside == 1) {
966  para[i] = mc_thb_ph_L[i];
967  error[i] = merr_thb_ph_L[i];
968  } else {
969  para[i] = mc_thb_ph_H[i];
970  error[i] = merr_thb_ph_H[i];
971  }
972  } //for
973  }
974  } // p-range
975 
976 
977  } //theta range
978 //--
979 // evaluate resolution in the unit of GeV.
980 //--
981  resol = para[0] + para[1] * pbuf + para[2] * pbuf * pbuf +
982  para[3] * pbuf * pbuf * pbuf;
983  resol = resol / 1000.;
984 
985 //--
986 // Evaluate the error od sigma using diagonal errors.
987 //--
988  double eresol = error[0] * error[0] + (pbuf * error[1]) * (pbuf * error[1])
989  + (pbuf * pbuf * error[2]) * (pbuf * pbuf * error[2])
990  + (pbuf * pbuf * pbuf * error[3]) * (pbuf * pbuf * pbuf * error[3]);
991  eresol = sqrt(eresol) / 1000.;
992 
993  // dout(Debugout::INFO,"B2BIIFixMdst_ecl") <<" theta="<<theta<<" p="<<pbuf
994  // <<" para="<< para[0]<<", "<< para[1]<<", "
995  // << para[2]<<", "<<para[3]
996  // <<" resol="<<resol<<" er="<<eresol<<std::endl;
997  return resol;
998  }
999 // Following lines were commented out, because asymmetric line shape is
1000 // not taken into account. 2003/Feb./17th KM.
1001 //==================================================
1002 // pi0 mass resolution function.
1003 // Original was coded by H.Hayashii
1004 // Implemented into B2BIIFixMdst_ecl in 2002/11/30 by KM
1005 //==================================================
1006 //*******************************
1007 // pi0resol: pi0 resolution function
1008 //
1009 //******************************
1010 // input p : pi0 lab momentum in GeV
1011 // theta : pi0 lab polar angle in the unit of Deg.
1012 // mcdata: =1/0 MC/Data
1013 // exp : experiment number
1014 // The unit of returned value is in GeV/c^2. (Note by KM).
1015 //static double pi0resol(double p, double theta, bool mcdata,int exp )
1016 //{
1017 //-----------------------------------------------------------
1018 // 20<=theta<=45
1019 //-----------------------------------------------------------
1020 // 0.15<= p <=0.8
1021 // -- data
1022 // const double data_thf_pl[]={9.4962, -21.906, 31.229, -13.424};
1023 // const double derr_thf_pl[]={0.4014, 0.978, 1.174, 0.786};
1024 // -- mc
1025 // const double mc_thf_pl[] ={4.3911, -2.121, 4.278, -1.696};
1026 // const double merr_thf_pl[]={0.258, 0.584, 0.617, 0.294};
1027 //---------------------------
1028 // 0.8<= p <=5.3
1029 // -- data
1030 // const double data_thf_ph[]={4.6545, 0.5212, 0.1370, -0.0214};
1031 // const double derr_thf_ph[]={0.189, 0.1866, 0.069, 0.099 };
1032 // -- mc
1033 // const double mc_thf_ph[] ={4.034, 0.8925, -0.040, -0.0034};
1034 // const double merr_thf_ph[]={0.205, 0.227, 0.088, 0.013 };
1035 //
1036 //-----------------------------------------------------------
1037 // 45<=theta<=100
1038 //-----------------------------------------------------------
1039 // 0.15<= p <=0.8
1040 // -- data
1041 // const double data_thm_pl[]={6.6403, -9.281, 16.179, -8.2410 };
1042 // const double derr_thm_pl[]={0.1353, 0.406, 0.579, 0.434 };
1043 // -- mc
1044 // const double mc_thm_pl[] ={3.7939, 1.877, -0.564, 0.098};
1045 // const double merr_thm_pl[]={0.145, 0.484, 0.705, 0.510};
1046 //---------------------------
1047 // 0.8<= p <=5.3
1048 // -- data
1049 // const double data_thm_ph[]={3.8323, 2.0226, -0.4476, 0.0833};
1050 // const double derr_thm_ph[]={0.1492, 0.2093, 0.0932, 0.0137};
1051 // -- mc
1052 // const double mc_thm_ph[] ={3.3595, 2.2174, -0.341, 0.0469};
1053 // const double merr_thm_ph[]={0.1311, 0.168, 0.077, 0.0131};
1054 //
1055 //-----------------------------------------------------------
1056 // 100<=theta<=150
1057 //-----------------------------------------------------------
1058 // 0.15<= p <=0.8
1059 // -- data
1060 // const double data_thb_pl[]={8.1320, -14.985, 22.567, -10.150 };
1061 // const double derr_thb_pl[]={0.2071, 0.537, 0.659, 0.458 };
1062 // -- mc
1063 // const double mc_thb_pl[] ={2.8217, 7.257, -7.936, 3.201};
1064 // const double merr_thb_pl[]={0.2085, 0.529, 0.642, 0.387};
1065 //---------------------------
1066 // 0.8<= p <=5.3
1067 // -- data
1068 // const double data_thb_ph[]={4.969 , 0.1843, 0.5477, -0.1011};
1069 // const double derr_thb_ph[]={0.164 , 0.2430, 0.1289, 0.0288};
1070 // -- mc
1071 // const double mc_thb_ph[] ={3.3595, 2.998 , -1.081, 0.1781};
1072 // const double merr_thb_ph[]={0.2625, 0.341, 0.180, 0.0411};
1073 //
1074 //
1075 // double resol;
1076 // double para[4];
1077 // double error[4];
1078 //
1079 // double pbuf = p;
1080 //--
1081 // theta< 45
1082 //---
1083 // if( theta <=45.) {
1084 //--
1085 // p-low
1086 //---
1087 // if( pbuf <= 0.8){
1088 // if(!mcdata){
1089 // //-- data
1090 // for( int i=0; i <= 3; i++){
1091 // para[i] = data_thf_pl[i];
1092 // error[i]= derr_thf_pl[i];}
1093 // } else{
1094 // //-- mc
1095 // for( int i=0; i <= 3; i++){
1096 // para[i] = mc_thf_pl[i];
1097 // error[i]= merr_thf_pl[i];}
1098 // }
1099 // } else {
1100 //--
1101 // p-high
1102 //---
1103 // use p= 5.2 value if p >= 5.2 for forward region
1104 // if(pbuf >= 5.2) {pbuf = 5.2;}
1105 // if(!mcdata){
1106 // //-- data
1107 // for( int i=0; i <= 3; i++){
1108 // para[i] = data_thf_ph[i];
1109 // error[i]= derr_thf_ph[i];}
1110 // } else{
1111 // //-- mc
1112 // for( int i=0; i <= 3; i++){
1113 // para[i] = mc_thf_ph[i];
1114 // error[i]= merr_thf_ph[i];}
1115 // }
1116 // } // p-range
1117 //--
1118 // 45< theta< 100
1119 //---
1120 // } else if( theta>45.&& theta<=100.){
1121 //--
1122 // p-low
1123 //---
1124 // if( pbuf <= 0.8){
1125 // if(!mcdata){
1126 // //-- data
1127 // for( int i=0; i <= 3; i++){
1128 // para[i] = data_thm_pl[i];
1129 // error[i]= derr_thm_pl[i];}
1130 // } else{
1131 // //-- mc
1132 // for( int i=0; i <= 3; i++){
1133 // para[i] = mc_thm_pl[i];
1134 // error[i]= merr_thm_pl[i];}
1135 // }
1136 // } else {
1137 //--
1138 // p-high
1139 //---
1140 // use p= 5.2 value if p >= 5.2 for middle region
1141 // if(pbuf >= 5.2) {pbuf = 5.2;}
1142 // if(!mcdata){
1143 // //-- data
1144 // for( int i=0; i <= 3; i++){
1145 // para[i] = data_thm_ph[i];
1146 // error[i]= derr_thm_ph[i];}
1147 // } else{
1148 // //-- mc
1149 // for( int i=0; i <= 3; i++){
1150 // para[i] = mc_thm_ph[i];
1151 // error[i]= merr_thm_ph[i];}
1152 // }
1153 // } // p-range
1154 //--
1155 // theta> 100
1156 //---
1157 // } else if( theta>100.){
1158 //--
1159 // p-low
1160 //---
1161 // if( pbuf <= 0.8){
1162 // if(!mcdata){
1163 // //-- data
1164 // for( int i=0; i <= 3; i++){
1165 // para[i] = data_thb_pl[i];
1166 // error[i]= derr_thb_pl[i];}
1167 // } else{
1168 // //-- mc
1169 // for( int i=0; i <= 3; i++){
1170 // para[i] = mc_thb_pl[i];
1171 // error[i]= merr_thb_pl[i];}
1172 // }
1173 // } else {
1174 //--
1175 // p-high
1176 //---
1177 // use p= 4.0 value if p >= 4.0 for backward region
1178 // if(pbuf >= 4.0) {pbuf = 4.0;}
1179 // if(!mcdata){
1180 // //-- data
1181 // for( int i=0; i <= 3; i++){
1182 // para[i] = data_thb_ph[i];
1183 // error[i]= derr_thb_ph[i];}
1184 // } else{
1185 // //-- mc
1186 // for( int i=0; i <= 3; i++){
1187 // para[i] = mc_thb_ph[i];
1188 // error[i]= merr_thb_ph[i];}
1189 // }
1190 // } // p-range
1191 //
1192 //
1193 // } //theta range
1194 //--
1195 // evaluate resolution in the unit of GeV.
1196 //--
1197 // resol = para[0] + para[1]*pbuf + para[2]*pbuf*pbuf +
1198 // para[3]*pbuf*pbuf*pbuf;
1199 // resol = resol/1000.;
1200 //
1201 //--
1202 // This error evaluation is not correct, since the off-diagonal part
1203 // is not took into account.
1204 //--
1205 // double eresol = error[0]*error[0] + (pbuf* error[1])*(pbuf*error[1])
1206 // + (pbuf*pbuf* error[2])* (pbuf*pbuf* error[2])
1207 // + (pbuf*pbuf*pbuf* error[3])*(pbuf*pbuf*pbuf* error[3]);
1208 // eresol = sqrt(eresol)/1000.;
1209 //
1210 // // dout(Debugout::INFO,"B2BIIFixMdst_ecl") <<" theta="<<theta<<" p="<<pbuf
1211 // // <<" para="<< para[0]<<", "<< para[1]<<", "
1212 // // << para[2]<<", "<<para[3]
1213 // // <<" resol="<<resol<<" er="<<eresol<<std::endl;
1214 // return resol;
1215 //}
1216 
1217 //======================================================
1218 // Correct photon's momenta and error matrix.
1219 //======================================================
1220  void B2BIIFixMdstModule::correct_ecl(int option, int version)
1221  {
1222 
1223  int ecl_threshold(0);
1224  if (m_reprocess_version == 0) {
1225  ecl_threshold = 1;
1226  } else if (m_reprocess_version == 1) {
1227  ecl_threshold = 0;
1228  }
1229 
1230  //dout(Debugout::INFO,"B2BIIFixMdst_ecl") << "This is version1.2.1" << std::endl;
1231  // dout(Debugout::INFO,"B2BIIFixMdst_ecl") << "entering correct_ecl()" << std::endl;
1232  // Check whether "version" is supported.
1233  if (version < 1 || version > 2) {
1234  B2WARNING("correct_ecl :: Warning! ");
1235  B2WARNING("version=" << version << " is not supported! ");
1236  B2WARNING("Exit doing nothing!");
1237  }
1238 
1239  //---------------------
1240  // Control sequence based on belle_event table.
1241  //---------------------
1242  Belle::Belle_event_Manager& evtmgr = Belle::Belle_event_Manager::get_manager();
1243  if (0 == evtmgr.count()) {
1244  // Do nothing if not exist, because we can not know whether
1245  // this event is Exp. data or MC.
1246  return;
1247  }
1248  //Check data or mc. Data:expmc=1, MC:expmc=2.
1249  int expmc = evtmgr[0].ExpMC();
1250  // Exp, Run, number to be given to correction func.
1251  int Eno = evtmgr[0].ExpNo();
1252  int Rno = evtmgr[0].RunNo();
1253 
1254  // If option is 0 and this event is MC, do nothing.
1255  // Modify : Exp.61 processed by b20080502, MC is also corrected.
1256  //org if( option==0 && expmc==2 )
1257  if (option == 0 && expmc == 2) {
1258  //if( Eno<55) //test
1259  //std::cout<<"ecl_threshold="<<ecl_threshold<<std::endl;
1260  // if( Eno<60 || ( Eno==61 && Rno>1209 ) /* 5S scan runs. */
1261  // || ( Eno==61 && m_correct_ecl_5s==1 ) /* Exp.61 old lib case.*/
1262  // || ( Eno==65 && Rno>1000 ) /* 1S runs.*/
1263  // || ( Eno==65 && m_correct_ecl_5s==1 ) /* Exp.65 old lib case.*/)
1264  if (ecl_threshold == 1) {
1265  return;
1266  }
1267  }
1268 
1269  // Correction curve version control.
1270  int i_previous(0);
1271  // Mdst_event_add table is there?
1272  Belle::Mdst_event_add_Manager& mevtmgr = Belle::Mdst_event_add_Manager::get_manager();
1273  // When Mdst_event_add exists,
1274  if (mevtmgr.count() > 0) {
1275  // Up to Exp.13, same treatment as V1.0 2001/Nov./30th
1276  //if( Eno <= 13 )
1277  // Modify; Exp.13 moved to same treatment as Exp.15 or later.
1278  // 20020524 Kenkichi Miyabayashi.
1279  if (Eno <= 11) {
1280  // If already the energy is corrected, exit.
1281  if (mevtmgr[0].flag(3) == 1) {
1282  return;
1283  }
1284  // Otherwise, set the flag at proper entity, flag(3).
1285  // Note that frag(1) and frag(2) have already been used
1286  // in scale_momenta().
1287  mevtmgr[0].flag(3, 1);
1288  }
1289  // Exp.13 and later...last update 2002/May./24 K.Miya.
1290  else { // Exp.15 and later...last update 2002/Feb./20
1291  // If well-established version was already applied, exit.
1292  if (mevtmgr[0].flag(3) >= version) {
1293  return;
1294  }
1295  // Otherwise, check the previous version.
1296  i_previous = mevtmgr[0].flag(3);
1297  if (i_previous == 0 || i_previous == 1) {
1298  mevtmgr[0].flag(3, version);
1299  } else { // Previous version is unsupported one.
1300  // Make Warning and exit.
1301  B2WARNING("correct_ecl :: Warning! ");
1302  B2WARNING("Previously, uncorrect version was used. ");
1303  B2WARNING(" Exit doing nothing");
1304  return;
1305  }
1306  }
1307  } else { // Create Mdst_event_add and set flag.
1308  Belle::Mdst_event_add& meadd = mevtmgr.add();
1309  // Up to Exp.13, same treatment as before.
1310  if (Eno <= 13) {
1311  meadd.flag(3, 1);
1312  } else { // For Exp.15, set version.
1313  meadd.flag(3, version);
1314  }
1315  }
1316 
1317  //--------------------------
1318  // If no ad_hoc correction has been applied so far or
1319  // correction curve is new version,
1320  // overwrite proper panther tables.
1321  //--------------------------
1322  // Get Mdst_ecl and Mdst_gamma.
1323  Belle::Mdst_ecl_Manager& eclmgr = Belle::Mdst_ecl_Manager::get_manager();
1324  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1325 
1326  // Get Mdst_ecl_aux. Added 20060413
1327  Belle::Mdst_ecl_aux_Manager& eclauxmgr = Belle::Mdst_ecl_aux_Manager::get_manager();
1328  std::vector<Belle::Mdst_ecl_aux>::iterator itaux = eclauxmgr.begin();
1329 
1330  // Correct energy and error matrix in Mdst_ecl.
1331  double factor;
1332  double factor13;
1333 
1334  for (std::vector<Belle::Mdst_ecl>::iterator itecl = eclmgr.begin();
1335  itecl != eclmgr.end(); ++itecl) {
1336  Belle::Mdst_ecl& shower = *itecl;
1337  // Shower energy and polar angle
1338  double shower_energy = shower.energy();
1339  double shower_theta = shower.theta();
1340 
1341  // Fix wrong calib. for Exp. 45.
1342  // if( Eno==45 )
1343  if (expmc == 1 && m_reprocess_version == 0 && Eno == 45) {
1344  int cellID = (*itaux).cId();
1345  double factor45 = ecl_adhoc_corr_45(Eno, Rno, cellID);
1346  //if( factor45!=1.0 )
1347  //{
1348  //int idecl=shower.get_ID();
1349  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<"Exp45 fix idecl="<<idecl<<" cellID="<<cellID
1350  //<<" factor45="<<factor45<<std::endl;
1351  //}
1352 
1353  shower.energy(shower.energy() / factor45);
1354  double factor452 = factor45 * factor45;
1355  shower.error(0, shower.error(0) / factor452); // Energy diagonal.
1356  shower.error(1, shower.error(1) / factor45); // Energy-phi.
1357  shower.error(3, shower.error(3) / factor45); // Energy-theta.
1358 
1359  // Here, also take care of Belle::Mdst_gamma.
1360  for (std::vector<Belle::Mdst_gamma>::iterator itgam45 = gammamgr.begin();
1361  itgam45 != gammamgr.end(); ++itgam45) {
1362  Belle::Mdst_gamma& gamma45 = *itgam45;
1363  if (gamma45.ecl().get_ID() == shower.get_ID()) {
1364  gamma45.px(gamma45.px() / factor45);
1365  gamma45.py(gamma45.py() / factor45);
1366  gamma45.pz(gamma45.pz() / factor45);
1367  //int idgam=gamma45.get_ID();
1368  //if( factor45!=1.0 )
1369  //{
1370  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<< "Exp45 fix idgamma="<<idgam<<std::endl;
1371  //}
1372  }
1373  }
1374  }
1375 
1376  // control sequence by option and expmc.
1377  // option=0
1378  switch (option) {
1379  case 0: // Option set as default.
1380  if (expmc == 2) { // This event record is MC data.
1381  // processed b20080331 or older, already skip.
1382  //KM std::cout<<"mdst_ecl"<<std::endl;
1383  factor
1384  = ecl_mcx3_corr(Eno, Rno, shower_energy, cos(shower_theta));
1385  } else { // Exp. data.
1386  //Original definition.
1387  //factor
1388  // = ecl_adhoc_corr( Eno, Rno,
1389  // shower_energy, cos(shower_theta));
1390  //Modified 20081222
1391  factor
1392  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1393  shower_energy, cos(shower_theta));
1394  // Special treatment of Exp.15 processed by b20020125.
1395  if (Eno == 15 && i_previous == 1) {
1396  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1397  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "ecl factor=" << factor << " " ;
1398  // Original scaling is done by Exp.13 curve.
1399  //Original definition.
1400  //factor13
1401  //= ecl_adhoc_corr( 13, 1,
1402  //shower_energy, cos(shower_theta));
1403  //Modified 20081222
1404  factor13
1405  = ecl_adhoc_corr(13, 1, ecl_threshold,
1406  shower_energy, cos(shower_theta));
1407  factor = factor / factor13;
1408  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor(rel)=" << factor << " ";
1409  }
1410  }
1411  break;
1412  case 1:
1413  if (expmc == 2) { // This event record is MC data.
1414  factor = 1.0 / mpi0pdg(shower.energy());
1415 
1416  } else { // This event record is real data.
1417  //Original definition.
1418  //factor
1419  // = ecl_adhoc_corr( Eno, Rno,
1420  // shower.energy(), cos(shower.theta()));
1421  //Modified 20081222
1422  factor
1423  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1424  shower.energy(), cos(shower.theta()));
1425  // Special treatment of Exp.15 processed by b20020125.
1426  if (Eno == 15 && i_previous == 1) {
1427  // Original scaling is done by Exp.13 curve.
1428  //Original definition.
1429  //factor13
1430  // = ecl_adhoc_corr( 13, 1,
1431  // shower_energy, cos(shower_theta));
1432  //Modified 20081222
1433  factor13
1434  = ecl_adhoc_corr(13, 1, ecl_threshold,
1435  shower_energy, cos(shower_theta));
1436  factor = factor / factor13;
1437  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor(rel)=" << factor << " ";
1438  }
1439 
1440  factor = factor / mpi0pdg(shower.energy());
1441  }
1442  break;
1443  default:
1444  factor = 1.0;
1445  }
1446 
1447  // energy correction.
1448  shower.energy(shower.energy() / factor);
1449 
1450  // error matrix correction.
1451  double factor2 = factor * factor;
1452  shower.error(0, shower.error(0) / factor2); // Energy diagonal.
1453  shower.error(1, shower.error(1) / factor); // Energy-phi.
1454  shower.error(3, shower.error(3) / factor); // Energy-theta.
1455 
1456  // Incriment Belle::Mdst_ecl_aux pointer.
1457  ++itaux;
1458  }
1459 
1460  // Correct energy in Belle::Mdst_gamma.
1461  for (std::vector<Belle::Mdst_gamma>::iterator itgam = gammamgr.begin();
1462  itgam != gammamgr.end(); ++itgam) {
1463  Belle::Mdst_gamma& gamma = *itgam;
1464  // Create the gamma's 3vector
1465  CLHEP::Hep3Vector gamma_3v(gamma.px(), gamma.py(), gamma.pz());
1466  double gamma_energy = gamma_3v.mag();
1467  double gamma_cos = gamma_3v.cosTheta();
1468 
1469  // control sequence by option and expmc.
1470  switch (option) {
1471  case 0: // Option as default.
1472  if (expmc == 2) { // This event record is MC data.
1473  // processed b20080331 or older, already skip.
1474  //KM std::cout<<"mdst_gamma"<<std::endl;
1475  factor
1476  = ecl_mcx3_corr(Eno, Rno, gamma_energy, gamma_cos);
1477  } else { // Exp. data.
1478  //Original definition
1479  //factor
1480  // = ecl_adhoc_corr( Eno, Rno,
1481  // gamma_energy, gamma_cos);
1482  //Modified 20081222
1483  factor
1484  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1485  gamma_energy, gamma_cos);
1486  // Special treatment of Exp.15 processed by b20020125.
1487  if (Eno == 15 && i_previous == 1) {
1488  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1489  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor=" << factor << " " ;
1490  // Original scaling is done by Exp.13 curve.
1491  //Original definition.
1492  //factor13
1493  // = ecl_adhoc_corr( 13, 1,
1494  // gamma_energy, gamma_cos);
1495  //Modified 20081222
1496  factor13
1497  = ecl_adhoc_corr(13, 1, ecl_threshold,
1498  gamma_energy, gamma_cos);
1499  factor = factor / factor13;
1500  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "gamma factor(rel)=" << factor << " " << std::endl;
1501  }
1502  }
1503  break;
1504  case 1:
1505  if (expmc == 2) {
1506  factor = 1.0 / mpi0pdg(gamma_3v.mag());
1507 
1508  } else {
1509  //Original definition.
1510  //factor
1511  // = ecl_adhoc_corr( Eno, Rno,
1512  // gamma_3v.mag(), gamma_3v.cosTheta());
1513  //Modified 20081222
1514  factor
1515  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1516  gamma_3v.mag(), gamma_3v.cosTheta());
1517  // Special treatment of Exp.15 processed by b20020125.
1518  if (Eno == 15 && i_previous == 1) {
1519  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1520  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor=" << factor << " " ;
1521  // Original scaling is done by Exp.13 curve.
1522  //Original definition.
1523  //factor13
1524  // = ecl_adhoc_corr( 13, 1,
1525  // gamma_energy, gamma_cos);
1526  //Modified 20081222
1527  factor13
1528  = ecl_adhoc_corr(13, 1, ecl_threshold,
1529  gamma_energy, gamma_cos);
1530  factor = factor / factor13;
1531  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "gamma factor(rel)=" << factor << " " << std::endl;
1532  }
1533  factor = factor / mpi0pdg(gamma_3v.mag());
1534  }
1535  break;
1536  default:
1537  factor = 1.0;
1538  }
1539 
1540  // factor should be used as a division.
1541  gamma.px(gamma.px() / factor);
1542  gamma.py(gamma.py() / factor);
1543  gamma.pz(gamma.pz() / factor);
1544  }
1545  }
1546 
1547 
1548 //=====================================================================
1549 //***** make_pi0.cc ********** created 2001/Jul./21st *****
1550 // Create Belle::Mdst_pi0 from Belle::Mdst_gamma and Belle::Mdst_ecl to let
1551 // people get mass-constraint fitted momentum of pi0
1552 // after ad_hoc correction.
1553 // input : option = 0; same as the existent Rececl_pi0
1554 // low_limit and up_limit are ignored.
1555 // option = 1; users can modify mass window as they
1556 // like. The boundary of window is defined
1557 // by low_limit ans up_limit (in GeV).
1558 // option = 2; users can modify mass window in the manner of
1559 // -Xsigma < Mgg - Mpi0 < +Xsigma. The value of
1560 // sigma(mass resolution) is calculated by pi0resol
1561 // function which is originally coded by Hayashii-san.
1562 //*********************************************************
1563  void B2BIIFixMdstModule::make_pi0(int option, double low_limit, double up_limit)
1564  {
1565  //---------------------
1566  // Check Exp. number and Data/MC.
1567  //---------------------
1568  Belle::Belle_event_Manager& evtmgr = Belle::Belle_event_Manager::get_manager();
1569  if (0 == evtmgr.count()) {
1570  // Do nothing if not exist, because we can not know whether
1571  // this event is Exp. data or MC.
1572  return;
1573  }
1574  //Check data or mc. Data:expmc=1, MC:expmc=2.
1575  int expmc = evtmgr[0].ExpMC();
1576  // Check Exp. number.
1577  int Eno = evtmgr[0].ExpNo();
1578  // Set mcdata frag to be compatible with pi0resol func.
1579  bool mcdata = 0;
1580  if (expmc == 2) {
1581  mcdata = 1;
1582  }
1583 
1584  // pi0 mass of PDG.
1585  const double mpi0_pdg = 0.1349739;
1586 
1587  // Default mass window.
1588  // const double low_default = 0.1178;
1589  // const double up_default = 0.1502;
1590  const double low_default = 0.080; // modified 20040606
1591  const double up_default = 0.180; // modified 20040606
1592 
1593  // Maximum iteration of fit.
1594  const int iter_max = 5;
1595 
1596  // Check whether proper option and mass window are given.
1597  switch (option) {
1598  case 0: // option=0 case, set default mass window.
1599  low_limit = low_default;
1600  up_limit = up_default;
1601  break;
1602  case 1: // option=1 case, check given mass window.
1603  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
1604  // If mass window is not correct, do nothing.
1605  B2ERROR("Invalid mass window between ");
1606  B2ERROR(" and " << up_limit);
1607  return;
1608  }
1609  break;
1610  case 2: // pi0 cand. are selected by -Xsigma < Mgg-Mpi0 < +Xsigma.
1611  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl") << "option=2 was selected." << std::endl;
1612  //dbg dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "mass window is " << low_limit << " & ";
1613  //dbg dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << up_limit << std::endl;
1614  if (0.0 < low_limit || up_limit < 0.0) {
1615  B2ERROR("option=2 was selected. ");
1616  B2ERROR("Invalid mass window! " << low_limit);
1617  B2ERROR(" sould be negative, or " << up_limit);
1618  B2ERROR(" should be positive.");
1619  return;
1620  }
1621  break;
1622  default: // Otherwise, invalid option, do nothing.
1623  B2ERROR("Invalid option=" << option);
1624  return;
1625  }
1626 
1627  // At first, clear already existing Belle::Mdst_pi0.
1628  Belle::Mdst_pi0_Manager::get_manager().remove();
1629 
1630  // Get Belle::Mdst_gamma for photon's momentum
1631  // and get Belle::Mdst_ecl for error matrix.
1632  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1633 
1634  // Re-allocate Belle::Mdst_pi0 table.
1635  Belle::Mdst_pi0_Manager& pi0mgr = Belle::Mdst_pi0_Manager::get_manager();
1636 
1637  // If Only one photon in the event, no need to do anything.
1638  if (gammamgr.count() < 2) {
1639  return;
1640  }
1641 
1642  // Make combination of two Belle::Mdst_gamma.
1643  for (std::vector<Belle::Mdst_gamma>::iterator itgamma = gammamgr.begin();
1644  itgamma != gammamgr.end(); ++itgamma) {
1645  Belle::Mdst_gamma& gamma1 = *itgamma;
1646  CLHEP::Hep3Vector gamma1_3v(gamma1.px(), gamma1.py(), gamma1.pz());
1647  CLHEP::HepLorentzVector gamma1_lv(gamma1_3v, gamma1_3v.mag());
1648 
1649  Belle::Mdst_ecl& ecl1 = gamma1.ecl();
1650 
1651  for (std::vector<Belle::Mdst_gamma>::iterator jtgamma = itgamma + 1;
1652  jtgamma != gammamgr.end(); ++jtgamma) {
1653  Belle::Mdst_gamma& gamma2 = *jtgamma;
1654  CLHEP::Hep3Vector gamma2_3v(gamma2.px(), gamma2.py(), gamma2.pz());
1655  CLHEP::HepLorentzVector gamma2_lv(gamma2_3v, gamma2_3v.mag());
1656 
1657  Belle::Mdst_ecl& ecl2 = gamma2.ecl();
1658 
1659  // Invariant mass before const. fit.
1660  CLHEP::HepLorentzVector gamgam_lv = gamma1_lv + gamma2_lv;
1661  const double mass_before = gamgam_lv.mag();
1662 
1663  // In the case of option=0 or 1, criteria is controlled
1664  // by the inv. mass.
1665  double mass_ctrl{0};
1666  if (option == 0 || option == 1) {
1667  mass_ctrl = mass_before;
1668  }
1669  if (option == 2) {
1670  // Asymmetric lineshape is taken into account.
1671  if (mass_before * mpi0pdg(mass_before) < mpi0_pdg) {
1672  mass_ctrl = (mass_before * mpi0pdg(mass_before) - mpi0_pdg) /
1673  pi0resol(gamgam_lv.vect().mag(),
1674  gamgam_lv.theta() * 180. / M_PI, "lower",
1675  mcdata, Eno, option);
1676  } else {
1677  mass_ctrl = (mass_before * mpi0pdg(mass_before) - mpi0_pdg) /
1678  pi0resol(gamgam_lv.vect().mag(),
1679  gamgam_lv.theta() * 180. / M_PI, "higher",
1680  mcdata, Eno, option);
1681  }
1682  }
1683 
1684  // If invariant mass is inside the window,
1685  // if( low_limit < mass_before && mass_before < up_limit ) // old.
1686  if (low_limit < mass_ctrl && mass_ctrl < up_limit) {
1687  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<"mass="<<mass_before;
1688  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" p="<<gamgam_lv.vect().mag()<<" theta=";
1689  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<gamgam_lv.theta()<<" mcdata="<<mcdata;
1690  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" higher="<<
1691  // pi0resol(gamgam_lv.vect().mag(),
1692  // gamgam_lv.theta()*180./M_PI,
1693  // "higher",mcdata, Eno, option )<<" or ";
1694  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" lower="<<
1695  // pi0resol(gamgam_lv.vect().mag(),
1696  // gamgam_lv.theta()*180./M_PI,
1697  // "lower",mcdata, Eno, option )<<std::endl;
1698  // Error matrix(covariant matrix)
1699  CLHEP::HepMatrix V(6, 6, 0);
1700 
1701  V[0][0] = ecl1.error(0);
1702  V[0][1] = V[1][0] = ecl1.error(1);
1703  V[1][1] = ecl1.error(2);
1704  V[0][2] = V[2][0] = ecl1.error(3);
1705  V[1][2] = V[2][1] = ecl1.error(4);
1706  V[2][2] = ecl1.error(5);
1707  V[3][3] = ecl2.error(0);
1708  V[3][4] = V[4][3] = ecl2.error(1);
1709  V[4][4] = ecl2.error(2);
1710  V[3][5] = V[5][3] = ecl2.error(3);
1711  V[4][5] = V[5][4] = ecl2.error(4);
1712  V[5][5] = ecl2.error(5);
1713 
1714  // Measurements; i.e. initial parameters.
1715  CLHEP::HepMatrix y0(6, 1);
1716  y0[0][0] = ecl1.energy();
1717  y0[1][0] = ecl1.phi();
1718  y0[2][0] = ecl1.theta();
1719  y0[3][0] = ecl2.energy();
1720  y0[4][0] = ecl2.phi();
1721  y0[5][0] = ecl2.theta();
1722 
1723  // Copy them to proper matrix which is given to fit.
1724  CLHEP::HepMatrix y(y0);
1725  // Delivative.
1726  CLHEP::HepMatrix Dy(6, 1, 0);
1727 
1728  int iter = 0;
1729  double f_old = DBL_MAX;
1730  double chi2_old = DBL_MAX;
1731  double /*mass_gg,*/ chi2 = DBL_MAX;
1732  bool exit_flag = false;
1733 
1734  // Set parameters to decide whether converged.
1735  const double Df_limit = 0.1;
1736  const double Dchi2_limit = 0.1;
1737  // Do mass constraint fit; iterate until convergence.
1738  while (1) {
1739  const double& E1 = y[0][0];
1740  const double& E2 = y[3][0];
1741  const double sin_theta1 = std::sin(y[2][0]);
1742  const double cos_theta1 = std::cos(y[2][0]);
1743  const double sin_theta2 = std::sin(y[5][0]);
1744  const double cos_theta2 = std::cos(y[5][0]);
1745  const double Dphi = y[1][0] - y[4][0];
1746  const double cos_Dphi = std::cos(Dphi);
1747  const double sin_Dphi = std::sin(Dphi);
1748  const double open_angle = sin_theta1 * sin_theta2 * cos_Dphi
1749  + cos_theta1 * cos_theta2;
1750  const double mass2_gg = 2 * E1 * E2 * (1 - open_angle);
1751  //mass_gg = (mass2_gg > 0) ? std::sqrt(mass2_gg) : -std::sqrt(-mass2_gg);
1752 
1753  // No more iteration. Break to return.
1754  if (exit_flag || ++iter > iter_max)
1755  break;
1756 
1757  // constraint
1758  CLHEP::HepMatrix f(1, 1);
1759  f[0][0] = mass2_gg - (mpi0_pdg * mpi0_pdg);
1760 
1761  // dG/dq_i, G = (M_gg - M_pi0) = 0 is the constraint.
1762  CLHEP::HepMatrix B(1, 6);
1763  B[0][0] = mass2_gg / E1;
1764  B[0][1] = 2 * E1 * E2 * sin_theta1 * sin_theta2 * sin_Dphi;
1765  B[0][2] = 2 * E1 * E2 * (-cos_theta1 * sin_theta2 * cos_Dphi
1766  + sin_theta1 * cos_theta2);
1767  B[0][3] = mass2_gg / E2;
1768  B[0][4] = -B[0][1];
1769  B[0][5] = 2 * E1 * E2 * (-sin_theta1 * cos_theta2 * cos_Dphi
1770  + cos_theta1 * sin_theta2);
1771 
1772  const double sigma2_mass2_gg = (B * V * B.T())[0][0];
1773 
1774  // Proceed one step to get newer value y.
1775  Dy = V * B.T() * (B * Dy - f) / sigma2_mass2_gg;
1776  y = y0 + Dy;
1777  int ierr;
1778  chi2 = (Dy.T() * V.inverse(ierr) * Dy)[0][0];
1779  double Dchi2 = fabs(chi2 - chi2_old);
1780  chi2_old = chi2;
1781  double Df = fabs(f[0][0] - f_old);
1782  f_old = f[0][0];
1783 
1784  // When chi-sq. change is small enough and mass is
1785  if (Dchi2 < Dchi2_limit && Df < Df_limit)
1786  exit_flag = true;;
1787 
1788  }
1789  // Sort fitted result into proper variables.
1790  double pi0_E = y[0][0] + y[3][0];
1791  double pi0_px = y[0][0] * std::cos(y[1][0]) * std::sin(y[2][0])
1792  + y[3][0] * std::cos(y[4][0]) * std::sin(y[5][0]);
1793  double pi0_py = y[0][0] * std::sin(y[1][0]) * std::sin(y[2][0])
1794  + y[3][0] * std::sin(y[4][0]) * std::sin(y[5][0]);
1795  double pi0_pz = y[0][0] * std::cos(y[2][0]) + y[3][0] * std::cos(y[5][0]);
1796  double pi0_mass = mass_before;
1797 
1798  // m_chi2 = chi2; // protect...
1799  double pi0_chi2 = (chi2 > FLT_MAX) ? FLT_MAX : chi2;
1800 
1801  // Fill Belle::Mdst_pi0 based on the fit result.
1802  Belle::Mdst_pi0& pi0 = pi0mgr.add();
1803  pi0.gamma(0, gamma1);
1804  pi0.gamma(1, gamma2);
1805  pi0.px(pi0_px);
1806  pi0.py(pi0_py);
1807  pi0.pz(pi0_pz);
1808  pi0.energy(pi0_E);
1809  pi0.mass(pi0_mass);
1810  pi0.chisq(pi0_chi2);
1811  }
1812  }
1813  }
1814  }
1815 
1816  void B2BIIFixMdstModule::make_pi0_primary_vertex(int option, double low_limit, double up_limit,
1817  const HepPoint3D& /*epvtx*/,
1818  const CLHEP::HepSymMatrix& epvtx_err)
1819  {
1820 #if 0
1821  // pi0 mass of PDG.
1822  const double mpi0_pdg = 0.1349739;
1823 
1824  // Default mass window.
1825  const double low_default = 0.1178;
1826  const double up_default = 0.1502;
1827 
1828  // Maximum iteration of fit.
1829  const int iter_max = 5;
1830 
1831  // Check whether proper option and mass window are given.
1832  switch (option) {
1833  case 0: // option=0 case, set default mass window.
1834  low_limit = low_default;
1835  up_limit = up_default;
1836  break;
1837  case 1: // option=1 case, check given mass window.
1838  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
1839  // If mass window is not correct, do nothing.
1840  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << "Invalid mass window between " << low_limit;
1841  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << " and " << up_limit << std::endl;
1842  return;
1843  }
1844  break;
1845  default: // Otherwise, invalid option, do nothing.
1846  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << "Invalid option=" << option << std::endl;
1847  return;
1848  }
1849 
1850  // At first, clear already existing Belle::Mdst_pi0.
1851  Belle::Mdst_pi0_Manager::get_manager().remove();
1852 
1853  // Get Belle::Mdst_gamma for photon's momentum
1854  // and get Belle::Mdst_ecl for error matrix.
1855  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1856 
1857  // Re-allocate Belle::Mdst_pi0 table.
1858  Belle::Mdst_pi0_Manager& pi0mgr = Belle::Mdst_pi0_Manager::get_manager();
1859 
1860  // If Only one photon in the event, no need to do anything.
1861  if (gammamgr.count() < 2) {
1862  return;
1863  }
1864 
1865  // Make combination of two Belle::Mdst_gamma.
1866  for (std::vector<Belle::Mdst_gamma>::iterator itgamma = gammamgr.begin();
1867  itgamma != gammamgr.end(); ++itgamma) {
1868  Belle::Mdst_gamma& gamma1 = *itgamma;
1869  CLHEP::Hep3Vector gamma1_3v(gamma1.px(), gamma1.py(), gamma1.pz());
1870  CLHEP::HepLorentzVector gamma1_lv(gamma1_3v, gamma1_3v.mag());
1871 
1872  Belle::Mdst_ecl& ecl1 = gamma1.ecl();
1873 
1874  const double r_i = ecl1.r();
1875  const double dzr_i = std::sqrt(epvtx_err[2][2]) / r_i;
1876 
1877  const double theta_i0 = ecl1.theta();
1878  const double sin_th_i0 = std::sin(theta_i0);
1879 
1880  for (std::vector<Belle::Mdst_gamma>::iterator jtgamma = itgamma + 1;
1881  jtgamma != gammamgr.end(); ++jtgamma) {
1882  Belle::Mdst_gamma& gamma2 = *jtgamma;
1883  CLHEP::Hep3Vector gamma2_3v(gamma2.px(), gamma2.py(), gamma2.pz());
1884  CLHEP::HepLorentzVector gamma2_lv(gamma2_3v, gamma2_3v.mag());
1885 
1886  Belle::Mdst_ecl& ecl2 = gamma2.ecl();
1887 
1888  // Invariant mass before const. fit.
1889  CLHEP::HepLorentzVector gamgam_lv = gamma1_lv + gamma2_lv;
1890  const double mass_before = gamgam_lv.mag();
1891 
1892  // If invariant mass is inside the window,
1893  if (low_limit < mass_before && mass_before < up_limit) {
1894  // Error matrix(covariant matrix)
1895  CLHEP::HepMatrix V(6, 6, 0);
1896 
1897  V[0][0] = ecl1.error(0);
1898  V[0][1] = V[1][0] = ecl1.error(1);
1899  V[1][1] = ecl1.error(2);
1900  V[0][2] = V[2][0] = ecl1.error(3);
1901  V[1][2] = V[2][1] = ecl1.error(4);
1902  V[2][2] = ecl1.error(5);
1903  V[3][3] = ecl2.error(0);
1904  V[3][4] = V[4][3] = ecl2.error(1);
1905  V[4][4] = ecl2.error(2);
1906  V[3][5] = V[5][3] = ecl2.error(3);
1907  V[4][5] = V[5][4] = ecl2.error(4);
1908  V[5][5] = ecl2.error(5);
1909 
1910  // Correlation term
1911  const double r_j = ecl2.r();
1912  const double dzr_j = sqrt(epvtx_err[2][2]) / r_j;
1913 
1914  const double theta_j0 = ecl2.theta();
1915  const double sin_th_j0 = std::sin(theta_j0);
1916 
1917  V[2][5] = V[5][2] = dzr_i * sin_th_i0 * dzr_j * sin_th_j0;
1918  // Measurements; i.e. initial parameters.
1919  CLHEP::HepMatrix y0(6, 1);
1920  y0[0][0] = ecl1.energy();
1921  y0[1][0] = ecl1.phi();
1922  y0[2][0] = ecl1.theta();
1923  y0[3][0] = ecl2.energy();
1924  y0[4][0] = ecl2.phi();
1925  y0[5][0] = ecl2.theta();
1926 
1927  // Copy them to proper matrix which is given to fit.
1928  CLHEP::HepMatrix y(y0);
1929  // Delivative.
1930  CLHEP::HepMatrix Dy(6, 1, 0);
1931 
1932  int iter = 0;
1933  double Df, f_old = DBL_MAX;
1934  double Dchi2, chi2_old = DBL_MAX;
1935  double mass_gg, chi2 = DBL_MAX;
1936  bool exit_flag = false;
1937 
1938  // Set parameters to decide whether converged.
1939  const double Df_limit = 0.1;
1940  const double Dchi2_limit = 0.1;
1941  // Do mass constraint fit; iterate until convergence.
1942  while (1) {
1943  const double& E1 = y[0][0];
1944  const double& E2 = y[3][0];
1945  const double sin_theta1 = std::sin(y[2][0]);
1946  const double cos_theta1 = std::cos(y[2][0]);
1947  const double sin_theta2 = std::sin(y[5][0]);
1948  const double cos_theta2 = std::cos(y[5][0]);
1949  const double Dphi = y[1][0] - y[4][0];
1950  const double cos_Dphi = std::cos(Dphi);
1951  const double sin_Dphi = std::sin(Dphi);
1952  const double open_angle = sin_theta1 * sin_theta2 * cos_Dphi
1953  + cos_theta1 * cos_theta2;
1954  const double mass2_gg = 2 * E1 * E2 * (1 - open_angle);
1955  mass_gg = (mass2_gg > 0) ? std::sqrt(mass2_gg) : -std::sqrt(-mass2_gg);
1956 
1957  // No more iteration. Break to return.
1958  if (exit_flag || ++iter > iter_max)
1959  break;
1960 
1961  // constraint
1962  CLHEP::HepMatrix f(1, 1);
1963  f[0][0] = mass2_gg - (mpi0_pdg * mpi0_pdg);
1964 
1965  // dG/dq_i, G = (M_gg - M_pi0) = 0 is the constraint.
1966  CLHEP::HepMatrix B(1, 6);
1967  B[0][0] = mass2_gg / E1;
1968  B[0][1] = 2 * E1 * E2 * sin_theta1 * sin_theta2 * sin_Dphi;
1969  B[0][2] = 2 * E1 * E2 * (-cos_theta1 * sin_theta2 * cos_Dphi
1970  + sin_theta1 * cos_theta2);
1971  B[0][3] = mass2_gg / E2;
1972  B[0][4] = -B[0][1];
1973  B[0][5] = 2 * E1 * E2 * (-sin_theta1 * cos_theta2 * cos_Dphi
1974  + cos_theta1 * sin_theta2);
1975 
1976  const double sigma2_mass2_gg = (B * V * B.T())[0][0];
1977 
1978  // Proceed one step to get newer value y.
1979  Dy = V * B.T() * (B * Dy - f) / sigma2_mass2_gg;
1980  y = y0 + Dy;
1981  int ierr;
1982  chi2 = (Dy.T() * V.inverse(ierr) * Dy)[0][0];
1983  Dchi2 = fabs(chi2 - chi2_old);
1984  chi2_old = chi2;
1985  Df = fabs(f[0][0] - f_old);
1986  f_old = f[0][0];
1987 
1988  // When chi-sq. change is small enough and mass is
1989  if (Dchi2 < Dchi2_limit && Df < Df_limit)
1990  exit_flag = true;;
1991 
1992  }
1993  // Sort fitted result into proper variables.
1994  double pi0_E = y[0][0] + y[3][0];
1995  double pi0_px = y[0][0] * std::cos(y[1][0]) * std::sin(y[2][0])
1996  + y[3][0] * std::cos(y[4][0]) * std::sin(y[5][0]);
1997  double pi0_py = y[0][0] * std::sin(y[1][0]) * std::sin(y[2][0])
1998  + y[3][0] * std::sin(y[4][0]) * std::sin(y[5][0]);
1999  double pi0_pz = y[0][0] * std::cos(y[2][0]) + y[3][0] * std::cos(y[5][0]);
2000  double pi0_mass = mass_before;
2001 
2002  // m_chi2 = chi2; // protect...
2003  double pi0_chi2 = (chi2 > FLT_MAX) ? FLT_MAX : chi2;
2004 
2005  // Fill Belle::Mdst_pi0 based on the fit result.
2006  Belle::Mdst_pi0& pi0 = pi0mgr.add();
2007  pi0.gamma(0, gamma1);
2008  pi0.gamma(1, gamma2);
2009  pi0.px(pi0_px);
2010  pi0.py(pi0_py);
2011  pi0.pz(pi0_pz);
2012  pi0.energy(pi0_E);
2013  pi0.mass(pi0_mass);
2014  pi0.chisq(pi0_chi2);
2015  }
2016  }
2017  }
2018 #else
2019  // pi0 mass of PDG. ;
2020  const double mpi0_pdg = 0.1349739;
2021  const double m2pi0_pdg = mpi0_pdg * mpi0_pdg;
2022  // Default mass window.;
2023  const double low_default = 0.1178;
2024  const double up_default = 0.1502;
2025  // Maximum iteration of fit.;
2026  const int iter_max = 5;
2027 
2028  // Check whether proper option and mass window are given.;
2029  switch (option) {
2030  case 0: // option=0 case, set default mass window.;
2031  low_limit = low_default;
2032  up_limit = up_default;
2033  break;
2034  case 1: // option=1 case, check given mass window.;
2035  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
2036  // If mass window is not correct, do nothing.;
2037  B2ERROR("Invalid mass window between " << low_limit);
2038  B2ERROR(" and " << up_limit);
2039  return;
2040  }
2041  break;
2042  default: // Otherwise, invalid option, do nothing.;
2043  B2ERROR("Invalid option=" << option);
2044  return;
2045  }
2046  Belle::Mdst_pi0_Manager& pi0_mgr = Belle::Mdst_pi0_Manager::get_manager();
2047  Belle::Mdst_gamma_Manager& gamma_mgr = Belle::Mdst_gamma_Manager::get_manager();
2048 
2049  // At first, clear already existing Belle::Mdst_pi0.;
2050  Belle::Mdst_pi0_Manager::get_manager().remove();
2051  // If Only one photon in the event, no need to do anything.;
2052  if (gamma_mgr.count() < 2) {
2053  return;
2054  }
2055 
2056  for (std::vector<Belle::Mdst_gamma>::iterator i = gamma_mgr.begin();
2057  i != gamma_mgr.end(); ++i) {
2058  const Belle::Mdst_gamma& gamma_i = *i;
2059  if (!gamma_i.ecl()) {
2060  continue;
2061  }
2062  const Belle::Mdst_ecl& ecl_i = gamma_i.ecl();
2063  const double r_i = ecl_i.r();
2064  const double e_i0 = ecl_i.energy();
2065  const double phi_i0 = ecl_i.phi();
2066  const double theta_i0 = ecl_i.theta();
2067 
2068  const double sin_th_i0 = std::sin(theta_i0);
2069  const double cos_th_i0 = std::cos(theta_i0);
2070 
2071  CLHEP::HepSymMatrix err_i(3, 0);
2072  err_i[0][0] = ecl_i.error(0);
2073  err_i[1][0] = ecl_i.error(1); err_i[1][1] = ecl_i.error(2);
2074  err_i[2][0] = ecl_i.error(3); err_i[2][1] = ecl_i.error(4); err_i[2][2] = ecl_i.error(5);
2075 
2076  const double dzr_i = std::sqrt(epvtx_err[2][2]) / r_i;
2077 
2078  for (std::vector<Belle::Mdst_gamma>::iterator j = i + 1; j != gamma_mgr.end(); ++j) {
2079  const Belle::Mdst_gamma& gamma_j = *j;
2080  if (!gamma_j.ecl()) {
2081  continue;
2082  }
2083  const Belle::Mdst_ecl& ecl_j = gamma_j.ecl();
2084  const double r_j = ecl_j.r();
2085  const double e_j0 = ecl_j.energy();
2086  const double phi_j0 = ecl_j.phi();
2087  const double theta_j0 = ecl_j.theta();
2088 
2089  const double sin_th_j0 = std::sin(theta_j0);
2090  const double cos_th_j0 = std::cos(theta_j0);
2091 
2092  CLHEP::HepSymMatrix err_j(3, 0);
2093  err_j[0][0] = ecl_j.error(0);
2094  err_j[1][0] = ecl_j.error(1); err_j[1][1] = ecl_j.error(2);
2095  err_j[2][0] = ecl_j.error(3); err_j[2][1] = ecl_j.error(4); err_j[2][2] = ecl_j.error(5);
2096 
2097  const double dzr_j = sqrt(epvtx_err[2][2]) / r_j;
2098 
2099  // -------------------------------------------------------------------- ;
2100 
2101  double dth_i_x_dth_j = dzr_i * sin_th_i0 * dzr_j * sin_th_j0;
2102 
2103  double det = err_i[2][2] * err_j[2][2] - dth_i_x_dth_j * dth_i_x_dth_j;
2104 
2105 
2106  double e_i = e_i0, phi_i = phi_i0, theta_i = theta_i0, sin_th_i = sin_th_i0, cos_th_i = cos_th_i0;
2107  double e_j = e_j0, phi_j = phi_j0, theta_j = theta_j0, sin_th_j = sin_th_j0, cos_th_j = cos_th_j0;
2108 
2109 
2110  double dphi = phi_i - phi_j;
2111  double sin_dphi = std::sin(dphi);
2112  double cos_dphi = std::cos(dphi);
2113 
2114  // Cos of opening angle ;
2115  double cos_th_ij = sin_th_i * sin_th_j * cos_dphi + cos_th_i * cos_th_j;
2116  double m2_ij = 2.0 * e_i * e_j * (1.0 - cos_th_ij);
2117  double mass = std::sqrt(m2_ij);
2118 
2119  if (mass < low_limit || mass > up_limit)
2120  continue;
2121 
2122  double de_i = 0.0, dphi_i = 0.0, dtheta_i = 0.0;
2123  double de_j = 0.0, dphi_j = 0.0, dtheta_j = 0.0;
2124 
2125  int it = 0;
2126 
2127  double f = 1.0e+6;
2128  double chisq = 1.0e+6;
2129  double df = 0.0, dchisq = 0.0;
2130 
2131  const double del_chisq = 0.1;
2132  const double del_f = 0.1;
2133  //const int icovar = 1;
2134 
2135  double mass0 = 0.0;
2136  double dmass0 = 0.0;
2137  //double dmass = 0.0;
2138  //double sdev = 0.0;
2139 
2140  do {
2141 
2142  const double dcos_th_ij_dphi_i = -sin_th_i * sin_th_j * sin_dphi;
2143  const double dcos_th_ij_dphi_j = sin_th_i * sin_th_j * sin_dphi;
2144  const double dcos_th_ij_dth_i = cos_th_i * sin_th_j * cos_dphi - sin_th_i * cos_th_j;
2145  const double dcos_th_ij_dth_j = cos_th_j * sin_th_i * cos_dphi - sin_th_j * cos_th_i;
2146 
2147  const double dm2_de_i = m2_ij / e_i;
2148  const double dm2_de_j = m2_ij / e_j;
2149 
2150  const double dm2_dcos_th_ij = -2.0 * e_i * e_j;
2151 
2152  const double dm2_dphi_i = dm2_dcos_th_ij * dcos_th_ij_dphi_i;
2153  const double dm2_dphi_j = dm2_dcos_th_ij * dcos_th_ij_dphi_j;
2154  const double dm2_dth_i = dm2_dcos_th_ij * dcos_th_ij_dth_i;
2155  const double dm2_dth_j = dm2_dcos_th_ij * dcos_th_ij_dth_j;
2156 
2157  double dm2_2
2158  = dm2_de_i * dm2_de_i * err_i[0][0]
2159  + dm2_de_j * dm2_de_j * err_j[0][0]
2160  + dm2_dphi_i * dm2_dphi_i * err_i[1][0]
2161  + dm2_dphi_j * dm2_dphi_j * err_j[1][0]
2162  + dm2_dth_i * dm2_dth_i * err_i[2][0]
2163  + dm2_dth_j * dm2_dth_j * err_j[2][0]
2164  + 2.0 * dm2_dth_i * dm2_dth_j * dth_i_x_dth_j;
2165 
2166  if (dm2_2 < 0.0) {
2167  dm2_2 = 0.0;
2168  }
2169  if (it++ == 0) {
2170  mass0 = mass;
2171  dmass0 = std::sqrt(dm2_2);
2172  /*if (mass0 > 0.0) {
2173  dmass = 0.5 * dmass0 / mass;
2174  }*/
2175  }
2176  //const double residual = mass - mpi0_pdg;
2177  //const double pull = residual / dmass;
2178  if (it >= iter_max ||
2179  mass < low_default || mass > up_default) {
2180  it = -it;
2181  //sdev = (dmass > 0.0) ? pull : -100;
2182  break;
2183  }
2184 
2185  const double del_m2
2186  = (m2pi0_pdg - m2_ij
2187  + dm2_de_i * de_i
2188  + dm2_de_j * de_j
2189  + dm2_dphi_i * dphi_i
2190  + dm2_dphi_j * dphi_j
2191  + dm2_dth_i * dtheta_i
2192  + dm2_dth_j * dtheta_j) / dm2_2;
2193 
2194  de_i = del_m2 * dm2_de_i * err_i[0][0];
2195  de_j = del_m2 * dm2_de_j * err_j[0][0];
2196  dphi_i = del_m2 * dm2_dphi_i * err_i[1][1];
2197  dphi_j = del_m2 * dm2_dphi_j * err_j[1][1];
2198  dtheta_i = del_m2 * (dm2_dth_i * err_i[2][2] + dm2_dth_j * dth_i_x_dth_j);
2199  dtheta_j = del_m2 * (dm2_dth_j * err_j[2][2] + dm2_dth_i * dth_i_x_dth_j);
2200 
2201 
2202  dchisq = chisq;
2203  chisq
2204  = de_i * de_i / err_i[0][0]
2205  + de_j * de_j / err_j[0][0]
2206  + dphi_i * dphi_i / err_i[1][1]
2207  + dphi_j * dphi_j / err_j[1][1]
2208  + (dtheta_i * dtheta_i * err_i[2][2]
2209  + dtheta_j * dtheta_j * err_j[2][2]
2210  + 2.0 * dtheta_i * dtheta_j * dth_i_x_dth_j) / det;
2211  dchisq -= chisq;
2212 
2213 
2214  e_i = e_i0 + de_i;
2215  e_j = e_j0 + de_j;
2216  phi_i = phi_i0 + dphi_i;
2217  phi_j = phi_j0 + dphi_j;
2218  theta_i = theta_i0 + dtheta_i;
2219  theta_j = theta_j0 + dtheta_j;
2220 
2221  sin_th_i = std::sin(theta_i);
2222  cos_th_i = std::cos(theta_i);
2223  sin_th_j = std::sin(theta_j);
2224  cos_th_j = std::cos(theta_j);
2225 
2226  dth_i_x_dth_j = dzr_i * sin_th_i * dzr_j * sin_th_j;
2227  det = err_i[2][2] * err_j[2][2] - dth_i_x_dth_j * dth_i_x_dth_j;
2228 
2229  dphi = phi_i - phi_j;
2230  sin_dphi = std::sin(dphi);
2231  cos_dphi = std::cos(dphi);
2232 
2233  cos_th_ij = sin_th_i * sin_th_j * cos_dphi + cos_th_i * cos_th_j;
2234  m2_ij = 2.0 * e_i * e_j * (1.0 - cos_th_ij);
2235  mass = std::sqrt(m2_ij);
2236 
2237  df = f;
2238  f = fabs(m2pi0_pdg - m2_ij) / dmass0;
2239  df -= f;
2240  ++it;
2241  } while (std::fabs(df) > del_f || std::fabs(dchisq) > del_chisq);
2242  const double cos_phi_i = std::cos(phi_i);
2243  const double cos_phi_j = std::cos(phi_j);
2244  const double sin_phi_i = std::sin(phi_i);
2245  const double sin_phi_j = std::sin(phi_j);
2246 
2247  const CLHEP::HepLorentzVector p4_i(e_i * sin_th_i * cos_phi_i,
2248  e_i * sin_th_i * sin_phi_i, e_i * cos_th_i, e_i);
2249  const CLHEP::HepLorentzVector p4_j(e_j * sin_th_j * cos_phi_j,
2250  e_j * sin_th_j * sin_phi_j, e_j * cos_th_j, e_j);
2251 
2252  const CLHEP::HepLorentzVector p4_pi0(p4_i + p4_j);
2253 
2254  // Fill Belle::Mdst_pi0 based on the fit result.;
2255  Belle::Mdst_pi0& pi0 = pi0_mgr.add();
2256  pi0.gamma(0, gamma_i);
2257  pi0.gamma(1, gamma_j);
2258  pi0.px(p4_pi0.x());
2259  pi0.py(p4_pi0.y());
2260  pi0.pz(p4_pi0.z());
2261  pi0.energy(p4_pi0.e());
2262  pi0.mass(mass0);
2263  pi0.chisq(chisq);
2264  }
2265  }
2266 #endif
2267  return;
2268  }
2269 
2270  void B2BIIFixMdstModule::correct_ecl_primary_vertex(const HepPoint3D& epvtx,
2271  const CLHEP::HepSymMatrix& epvtx_err)
2272  {
2273 
2274  Belle::Mdst_gamma_Manager& Gamma = Belle::Mdst_gamma_Manager::get_manager();
2275 
2276  for (std::vector<Belle::Mdst_gamma>::iterator
2277  it = Gamma.begin(); it != Gamma.end(); ++it) {
2278  double r(it->ecl().r());
2279  double theta(it->ecl().theta());
2280  double phi(it->ecl().phi());
2281  // double stheta(std::sqrt(it->ecl().error(2)));
2282  // double sphi(std::sqrt(it->ecl().error(5)));
2283  double st(std::sin(theta));
2284  double ct(std::cos(theta));
2285  double sp(std::sin(phi));
2286  double cp(std::cos(phi));
2287  HepPoint3D gamma_pos(r * st * cp, r * st * sp, r * ct);
2288  CLHEP::Hep3Vector gamma_vec(gamma_pos - epvtx);
2289  double hsq(gamma_vec.perp2());
2290  double rsq(gamma_vec.mag2());
2291  double stheta_sq_new = it->ecl().error(5) + epvtx_err(3, 3) * (hsq / (rsq * rsq));
2292  CLHEP::Hep3Vector gamma_dir(gamma_vec.unit());
2293  double e(it->ecl().energy());
2294  it->px(e * gamma_dir.x());
2295  it->py(e * gamma_dir.y());
2296  it->pz(e * gamma_dir.z());
2297  it->ecl().theta(gamma_dir.theta());
2298  it->ecl().r(gamma_vec.mag());
2299  it->ecl().error(5, stheta_sq_new);
2300  }
2301  }
2302 
2303 
2305 } // namespace Belle
Belle2::ecl_mcx3_corr
static double ecl_mcx3_corr(int, int, double energy, double)
Correct energy scale (MC) to make pi0 peak nominal.
Definition: B2BIIFixMdstModule_ecl.cc:616
Belle2::ecl_adhoc_corr_45
static double ecl_adhoc_corr_45(int exp, int, int cid)
The function giving correction factor in Exp.45.
Definition: B2BIIFixMdstModule_ecl.cc:591
Belle2::B2BIIFixMdstModule::m_reprocess_version
int m_reprocess_version
Reprocess version (=0:old; =1:new)
Definition: B2BIIFixMdstModule.h:168
Belle2::B2BIIFixMdstModule::pi0resol
static double pi0resol(double, double, const char *, bool, int, int)
Treat pi0 mass width as a func.
Definition: B2BIIFixMdstModule_ecl.cc:678
Belle2::B2BIIFixMdstModule::correct_ecl
void correct_ecl(int, int)
Correct photon's momenta and error matrix.
Definition: B2BIIFixMdstModule_ecl.cc:1220
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::B2BIIFixMdstModule::make_pi0
void make_pi0(int, double, double)
Create Mdst_pi0 from Mdst_gamma and Mdst_ecl to let people get mass-constraint fitted momentum of pi0...
Definition: B2BIIFixMdstModule_ecl.cc:1563
Belle2::mpi0pdg
static double mpi0pdg(double Energy)
Make MC mass peak to PDG value.
Definition: B2BIIFixMdstModule_ecl.cc:635
HepGeom::Point3D< double >
Belle2::B2BIIFixMdstModule::correct_ecl_primary_vertex
void correct_ecl_primary_vertex(const HepPoint3D &, const CLHEP::HepSymMatrix &)
Correct ecl using primary vertex.
Definition: B2BIIFixMdstModule_ecl.cc:2270
Belle2::B2BIIFixMdstModule::make_pi0_primary_vertex
void make_pi0_primary_vertex(int, double, double, const HepPoint3D &, const CLHEP::HepSymMatrix &)
Fill Mdst_pi0 based on the fit result.
Definition: B2BIIFixMdstModule_ecl.cc:1816
Belle2::ecl_adhoc_corr
static double ecl_adhoc_corr(int Exp, int Run, int iflag05th, double Energy, double)
The function giving correction factor.
Definition: B2BIIFixMdstModule_ecl.cc:148