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