Belle II Software  release-08-01-10
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  const int bcid[] = {
602  243, 244, 316, 317, 318, 412, 413, 414, 508, 509, 510, 604, 605, 606
603  };
604 
605  const 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 <= 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 {
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  B2DEBUG(19, "B2BIIFixMdst_ecl" << LogVar("theta", theta) << LogVar("p", pbuf)
1001  << LogVar("para0", para[0]) << LogVar("para1", para[1]) << LogVar("para2", para[2]) << LogVar("para3", para[3])
1002  << LogVar("resol", resol) << LogVar("er", eresol));
1003  return resol;
1004  }
1005 // Following lines were commented out, because asymmetric line shape is
1006 // not taken into account. 2003/Feb./17th KM.
1007 //==================================================
1008 // pi0 mass resolution function.
1009 // Original was coded by H.Hayashii
1010 // Implemented into B2BIIFixMdst_ecl in 2002/11/30 by KM
1011 //==================================================
1012 //*******************************
1013 // pi0resol: pi0 resolution function
1014 //
1015 //******************************
1016 // input p : pi0 lab momentum in GeV
1017 // theta : pi0 lab polar angle in the unit of Deg.
1018 // mcdata: =1/0 MC/Data
1019 // exp : experiment number
1020 // The unit of returned value is in GeV/c^2. (Note by KM).
1021 //static double pi0resol(double p, double theta, bool mcdata,int exp )
1022 //{
1023 //-----------------------------------------------------------
1024 // 20<=theta<=45
1025 //-----------------------------------------------------------
1026 // 0.15<= p <=0.8
1027 // -- data
1028 // const double data_thf_pl[]={9.4962, -21.906, 31.229, -13.424};
1029 // const double derr_thf_pl[]={0.4014, 0.978, 1.174, 0.786};
1030 // -- mc
1031 // const double mc_thf_pl[] ={4.3911, -2.121, 4.278, -1.696};
1032 // const double merr_thf_pl[]={0.258, 0.584, 0.617, 0.294};
1033 //---------------------------
1034 // 0.8<= p <=5.3
1035 // -- data
1036 // const double data_thf_ph[]={4.6545, 0.5212, 0.1370, -0.0214};
1037 // const double derr_thf_ph[]={0.189, 0.1866, 0.069, 0.099 };
1038 // -- mc
1039 // const double mc_thf_ph[] ={4.034, 0.8925, -0.040, -0.0034};
1040 // const double merr_thf_ph[]={0.205, 0.227, 0.088, 0.013 };
1041 //
1042 //-----------------------------------------------------------
1043 // 45<=theta<=100
1044 //-----------------------------------------------------------
1045 // 0.15<= p <=0.8
1046 // -- data
1047 // const double data_thm_pl[]={6.6403, -9.281, 16.179, -8.2410 };
1048 // const double derr_thm_pl[]={0.1353, 0.406, 0.579, 0.434 };
1049 // -- mc
1050 // const double mc_thm_pl[] ={3.7939, 1.877, -0.564, 0.098};
1051 // const double merr_thm_pl[]={0.145, 0.484, 0.705, 0.510};
1052 //---------------------------
1053 // 0.8<= p <=5.3
1054 // -- data
1055 // const double data_thm_ph[]={3.8323, 2.0226, -0.4476, 0.0833};
1056 // const double derr_thm_ph[]={0.1492, 0.2093, 0.0932, 0.0137};
1057 // -- mc
1058 // const double mc_thm_ph[] ={3.3595, 2.2174, -0.341, 0.0469};
1059 // const double merr_thm_ph[]={0.1311, 0.168, 0.077, 0.0131};
1060 //
1061 //-----------------------------------------------------------
1062 // 100<=theta<=150
1063 //-----------------------------------------------------------
1064 // 0.15<= p <=0.8
1065 // -- data
1066 // const double data_thb_pl[]={8.1320, -14.985, 22.567, -10.150 };
1067 // const double derr_thb_pl[]={0.2071, 0.537, 0.659, 0.458 };
1068 // -- mc
1069 // const double mc_thb_pl[] ={2.8217, 7.257, -7.936, 3.201};
1070 // const double merr_thb_pl[]={0.2085, 0.529, 0.642, 0.387};
1071 //---------------------------
1072 // 0.8<= p <=5.3
1073 // -- data
1074 // const double data_thb_ph[]={4.969 , 0.1843, 0.5477, -0.1011};
1075 // const double derr_thb_ph[]={0.164 , 0.2430, 0.1289, 0.0288};
1076 // -- mc
1077 // const double mc_thb_ph[] ={3.3595, 2.998 , -1.081, 0.1781};
1078 // const double merr_thb_ph[]={0.2625, 0.341, 0.180, 0.0411};
1079 //
1080 //
1081 // double resol;
1082 // double para[4];
1083 // double error[4];
1084 //
1085 // double pbuf = p;
1086 //--
1087 // theta< 45
1088 //---
1089 // if( theta <=45.) {
1090 //--
1091 // p-low
1092 //---
1093 // if( pbuf <= 0.8){
1094 // if(!mcdata){
1095 // //-- data
1096 // for( int i=0; i <= 3; i++){
1097 // para[i] = data_thf_pl[i];
1098 // error[i]= derr_thf_pl[i];}
1099 // } else{
1100 // //-- mc
1101 // for( int i=0; i <= 3; i++){
1102 // para[i] = mc_thf_pl[i];
1103 // error[i]= merr_thf_pl[i];}
1104 // }
1105 // } else {
1106 //--
1107 // p-high
1108 //---
1109 // use p= 5.2 value if p >= 5.2 for forward region
1110 // if(pbuf >= 5.2) {pbuf = 5.2;}
1111 // if(!mcdata){
1112 // //-- data
1113 // for( int i=0; i <= 3; i++){
1114 // para[i] = data_thf_ph[i];
1115 // error[i]= derr_thf_ph[i];}
1116 // } else{
1117 // //-- mc
1118 // for( int i=0; i <= 3; i++){
1119 // para[i] = mc_thf_ph[i];
1120 // error[i]= merr_thf_ph[i];}
1121 // }
1122 // } // p-range
1123 //--
1124 // 45< theta< 100
1125 //---
1126 // } else if( theta>45.&& theta<=100.){
1127 //--
1128 // p-low
1129 //---
1130 // if( pbuf <= 0.8){
1131 // if(!mcdata){
1132 // //-- data
1133 // for( int i=0; i <= 3; i++){
1134 // para[i] = data_thm_pl[i];
1135 // error[i]= derr_thm_pl[i];}
1136 // } else{
1137 // //-- mc
1138 // for( int i=0; i <= 3; i++){
1139 // para[i] = mc_thm_pl[i];
1140 // error[i]= merr_thm_pl[i];}
1141 // }
1142 // } else {
1143 //--
1144 // p-high
1145 //---
1146 // use p= 5.2 value if p >= 5.2 for middle region
1147 // if(pbuf >= 5.2) {pbuf = 5.2;}
1148 // if(!mcdata){
1149 // //-- data
1150 // for( int i=0; i <= 3; i++){
1151 // para[i] = data_thm_ph[i];
1152 // error[i]= derr_thm_ph[i];}
1153 // } else{
1154 // //-- mc
1155 // for( int i=0; i <= 3; i++){
1156 // para[i] = mc_thm_ph[i];
1157 // error[i]= merr_thm_ph[i];}
1158 // }
1159 // } // p-range
1160 //--
1161 // theta> 100
1162 //---
1163 // } else if( theta>100.){
1164 //--
1165 // p-low
1166 //---
1167 // if( pbuf <= 0.8){
1168 // if(!mcdata){
1169 // //-- data
1170 // for( int i=0; i <= 3; i++){
1171 // para[i] = data_thb_pl[i];
1172 // error[i]= derr_thb_pl[i];}
1173 // } else{
1174 // //-- mc
1175 // for( int i=0; i <= 3; i++){
1176 // para[i] = mc_thb_pl[i];
1177 // error[i]= merr_thb_pl[i];}
1178 // }
1179 // } else {
1180 //--
1181 // p-high
1182 //---
1183 // use p= 4.0 value if p >= 4.0 for backward region
1184 // if(pbuf >= 4.0) {pbuf = 4.0;}
1185 // if(!mcdata){
1186 // //-- data
1187 // for( int i=0; i <= 3; i++){
1188 // para[i] = data_thb_ph[i];
1189 // error[i]= derr_thb_ph[i];}
1190 // } else{
1191 // //-- mc
1192 // for( int i=0; i <= 3; i++){
1193 // para[i] = mc_thb_ph[i];
1194 // error[i]= merr_thb_ph[i];}
1195 // }
1196 // } // p-range
1197 //
1198 //
1199 // } //theta range
1200 //--
1201 // evaluate resolution in the unit of GeV.
1202 //--
1203 // resol = para[0] + para[1]*pbuf + para[2]*pbuf*pbuf +
1204 // para[3]*pbuf*pbuf*pbuf;
1205 // resol = resol/1000.;
1206 //
1207 //--
1208 // This error evaluation is not correct, since the off-diagonal part
1209 // is not took into account.
1210 //--
1211 // double eresol = error[0]*error[0] + (pbuf* error[1])*(pbuf*error[1])
1212 // + (pbuf*pbuf* error[2])* (pbuf*pbuf* error[2])
1213 // + (pbuf*pbuf*pbuf* error[3])*(pbuf*pbuf*pbuf* error[3]);
1214 // eresol = sqrt(eresol)/1000.;
1215 //
1216 // // dout(Debugout::INFO,"B2BIIFixMdst_ecl") <<" theta="<<theta<<" p="<<pbuf
1217 // // <<" para="<< para[0]<<", "<< para[1]<<", "
1218 // // << para[2]<<", "<<para[3]
1219 // // <<" resol="<<resol<<" er="<<eresol<<std::endl;
1220 // return resol;
1221 //}
1222 
1223 //======================================================
1224 // Correct photon's momenta and error matrix.
1225 //======================================================
1226  void B2BIIFixMdstModule::correct_ecl(int option, int version)
1227  {
1228 
1229  int ecl_threshold(0);
1230  if (m_reprocess_version == 0) {
1231  ecl_threshold = 1;
1232  } else if (m_reprocess_version == 1) {
1233  ecl_threshold = 0;
1234  }
1235 
1236  //dout(Debugout::INFO,"B2BIIFixMdst_ecl") << "This is version1.2.1" << std::endl;
1237  // dout(Debugout::INFO,"B2BIIFixMdst_ecl") << "entering correct_ecl()" << std::endl;
1238  // Check whether "version" is supported.
1239  if (version < 1 || version > 2) {
1240  B2WARNING("correct_ecl :: Warning! ");
1241  B2WARNING("version=" << version << " is not supported! ");
1242  B2WARNING("Exit doing nothing!");
1243  }
1244 
1245  //---------------------
1246  // Control sequence based on belle_event table.
1247  //---------------------
1248  Belle::Belle_event_Manager& evtmgr = Belle::Belle_event_Manager::get_manager();
1249  if (0 == evtmgr.count()) {
1250  // Do nothing if not exist, because we can not know whether
1251  // this event is Exp. data or MC.
1252  return;
1253  }
1254  //Check data or mc. Data:expmc=1, MC:expmc=2.
1255  int expmc = evtmgr[0].ExpMC();
1256  // Exp, Run, number to be given to correction func.
1257  int Eno = evtmgr[0].ExpNo();
1258  int Rno = evtmgr[0].RunNo();
1259 
1260  // If option is 0 and this event is MC, do nothing.
1261  // Modify : Exp.61 processed by b20080502, MC is also corrected.
1262  //org if( option==0 && expmc==2 )
1263  if (option == 0 && expmc == 2) {
1264  //if( Eno<55) //test
1265  //std::cout<<"ecl_threshold="<<ecl_threshold<<std::endl;
1266  // if( Eno<60 || ( Eno==61 && Rno>1209 ) /* 5S scan runs. */
1267  // || ( Eno==61 && m_correct_ecl_5s==1 ) /* Exp.61 old lib case.*/
1268  // || ( Eno==65 && Rno>1000 ) /* 1S runs.*/
1269  // || ( Eno==65 && m_correct_ecl_5s==1 ) /* Exp.65 old lib case.*/)
1270  if (ecl_threshold == 1) {
1271  return;
1272  }
1273  }
1274 
1275  // Correction curve version control.
1276  int i_previous(0);
1277  // Mdst_event_add table is there?
1278  Belle::Mdst_event_add_Manager& mevtmgr = Belle::Mdst_event_add_Manager::get_manager();
1279  // When Mdst_event_add exists,
1280  if (mevtmgr.count() > 0) {
1281  // Up to Exp.13, same treatment as V1.0 2001/Nov./30th
1282  //if( Eno <= 13 )
1283  // Modify; Exp.13 moved to same treatment as Exp.15 or later.
1284  // 20020524 Kenkichi Miyabayashi.
1285  if (Eno <= 11) {
1286  // If already the energy is corrected, exit.
1287  if (mevtmgr[0].flag(3) == 1) {
1288  return;
1289  }
1290  // Otherwise, set the flag at proper entity, flag(3).
1291  // Note that frag(1) and frag(2) have already been used
1292  // in scale_momenta().
1293  mevtmgr[0].flag(3, 1);
1294  }
1295  // Exp.13 and later...last update 2002/May./24 K.Miya.
1296  else { // Exp.15 and later...last update 2002/Feb./20
1297  // If well-established version was already applied, exit.
1298  if (mevtmgr[0].flag(3) >= version) {
1299  return;
1300  }
1301  // Otherwise, check the previous version.
1302  i_previous = mevtmgr[0].flag(3);
1303  if (i_previous == 0 || i_previous == 1) {
1304  mevtmgr[0].flag(3, version);
1305  } else { // Previous version is unsupported one.
1306  // Make Warning and exit.
1307  B2WARNING("correct_ecl :: Warning! ");
1308  B2WARNING("Previously, uncorrect version was used. ");
1309  B2WARNING(" Exit doing nothing");
1310  return;
1311  }
1312  }
1313  } else { // Create Mdst_event_add and set flag.
1314  Belle::Mdst_event_add& meadd = mevtmgr.add();
1315  // Up to Exp.13, same treatment as before.
1316  if (Eno <= 13) {
1317  meadd.flag(3, 1);
1318  } else { // For Exp.15, set version.
1319  meadd.flag(3, version);
1320  }
1321  }
1322 
1323  //--------------------------
1324  // If no ad_hoc correction has been applied so far or
1325  // correction curve is new version,
1326  // overwrite proper panther tables.
1327  //--------------------------
1328  // Get Mdst_ecl and Mdst_gamma.
1329  Belle::Mdst_ecl_Manager& eclmgr = Belle::Mdst_ecl_Manager::get_manager();
1330  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1331 
1332  // Get Mdst_ecl_aux. Added 20060413
1333  Belle::Mdst_ecl_aux_Manager& eclauxmgr = Belle::Mdst_ecl_aux_Manager::get_manager();
1334  std::vector<Belle::Mdst_ecl_aux>::iterator itaux = eclauxmgr.begin();
1335 
1336  // Correct energy and error matrix in Mdst_ecl.
1337  double factor;
1338  double factor13;
1339 
1340  for (std::vector<Belle::Mdst_ecl>::iterator itecl = eclmgr.begin();
1341  itecl != eclmgr.end(); ++itecl) {
1342  Belle::Mdst_ecl& shower = *itecl;
1343  // Shower energy and polar angle
1344  double shower_energy = shower.energy();
1345  double shower_theta = shower.theta();
1346 
1347  // Fix wrong calib. for Exp. 45.
1348  // if( Eno==45 )
1349  if (expmc == 1 && m_reprocess_version == 0 && Eno == 45) {
1350  int cellID = (*itaux).cId();
1351  double factor45 = ecl_adhoc_corr_45(Eno, Rno, cellID);
1352  //if( factor45!=1.0 )
1353  //{
1354  //int idecl=shower.get_ID();
1355  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<"Exp45 fix idecl="<<idecl<<" cellID="<<cellID
1356  //<<" factor45="<<factor45<<std::endl;
1357  //}
1358 
1359  shower.energy(shower.energy() / factor45);
1360  double factor452 = factor45 * factor45;
1361  shower.error(0, shower.error(0) / factor452); // Energy diagonal.
1362  shower.error(1, shower.error(1) / factor45); // Energy-phi.
1363  shower.error(3, shower.error(3) / factor45); // Energy-theta.
1364 
1365  // Here, also take care of Belle::Mdst_gamma.
1366  for (std::vector<Belle::Mdst_gamma>::iterator itgam45 = gammamgr.begin();
1367  itgam45 != gammamgr.end(); ++itgam45) {
1368  Belle::Mdst_gamma& gamma45 = *itgam45;
1369  if (gamma45.ecl().get_ID() == shower.get_ID()) {
1370  gamma45.px(gamma45.px() / factor45);
1371  gamma45.py(gamma45.py() / factor45);
1372  gamma45.pz(gamma45.pz() / factor45);
1373  //int idgam=gamma45.get_ID();
1374  //if( factor45!=1.0 )
1375  //{
1376  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<< "Exp45 fix idgamma="<<idgam<<std::endl;
1377  //}
1378  }
1379  }
1380  }
1381 
1382  // control sequence by option and expmc.
1383  // option=0
1384  switch (option) {
1385  case 0: // Option set as default.
1386  if (expmc == 2) { // This event record is MC data.
1387  // processed b20080331 or older, already skip.
1388  //KM std::cout<<"mdst_ecl"<<std::endl;
1389  factor
1390  = ecl_mcx3_corr(Eno, Rno, shower_energy, cos(shower_theta));
1391  } else { // Exp. data.
1392  //Original definition.
1393  //factor
1394  // = ecl_adhoc_corr( Eno, Rno,
1395  // shower_energy, cos(shower_theta));
1396  //Modified 20081222
1397  factor
1398  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1399  shower_energy, cos(shower_theta));
1400  // Special treatment of Exp.15 processed by b20020125.
1401  if (Eno == 15 && i_previous == 1) {
1402  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1403  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "ecl factor=" << factor << " " ;
1404  // Original scaling is done by Exp.13 curve.
1405  //Original definition.
1406  //factor13
1407  //= ecl_adhoc_corr( 13, 1,
1408  //shower_energy, cos(shower_theta));
1409  //Modified 20081222
1410  factor13
1411  = ecl_adhoc_corr(13, 1, ecl_threshold,
1412  shower_energy, cos(shower_theta));
1413  factor = factor / factor13;
1414  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor(rel)=" << factor << " ";
1415  }
1416  }
1417  break;
1418  case 1:
1419  if (expmc == 2) { // This event record is MC data.
1420  factor = 1.0 / mpi0pdg(shower.energy());
1421 
1422  } else { // This event record is real data.
1423  //Original definition.
1424  //factor
1425  // = ecl_adhoc_corr( Eno, Rno,
1426  // shower.energy(), cos(shower.theta()));
1427  //Modified 20081222
1428  factor
1429  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1430  shower.energy(), cos(shower.theta()));
1431  // Special treatment of Exp.15 processed by b20020125.
1432  if (Eno == 15 && i_previous == 1) {
1433  // Original scaling is done by Exp.13 curve.
1434  //Original definition.
1435  //factor13
1436  // = ecl_adhoc_corr( 13, 1,
1437  // shower_energy, cos(shower_theta));
1438  //Modified 20081222
1439  factor13
1440  = ecl_adhoc_corr(13, 1, ecl_threshold,
1441  shower_energy, cos(shower_theta));
1442  factor = factor / factor13;
1443  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor(rel)=" << factor << " ";
1444  }
1445 
1446  factor = factor / mpi0pdg(shower.energy());
1447  }
1448  break;
1449  default:
1450  factor = 1.0;
1451  }
1452 
1453  // energy correction.
1454  shower.energy(shower.energy() / factor);
1455 
1456  // error matrix correction.
1457  double factor2 = factor * factor;
1458  shower.error(0, shower.error(0) / factor2); // Energy diagonal.
1459  shower.error(1, shower.error(1) / factor); // Energy-phi.
1460  shower.error(3, shower.error(3) / factor); // Energy-theta.
1461 
1462  // Incriment Belle::Mdst_ecl_aux pointer.
1463  ++itaux;
1464  }
1465 
1466  // Correct energy in Belle::Mdst_gamma.
1467  for (std::vector<Belle::Mdst_gamma>::iterator itgam = gammamgr.begin();
1468  itgam != gammamgr.end(); ++itgam) {
1469  Belle::Mdst_gamma& gamma = *itgam;
1470  // Create the gamma's 3vector
1471  CLHEP::Hep3Vector gamma_3v(gamma.px(), gamma.py(), gamma.pz());
1472  double gamma_energy = gamma_3v.mag();
1473  double gamma_cos = gamma_3v.cosTheta();
1474 
1475  // control sequence by option and expmc.
1476  switch (option) {
1477  case 0: // Option as default.
1478  if (expmc == 2) { // This event record is MC data.
1479  // processed b20080331 or older, already skip.
1480  //KM std::cout<<"mdst_gamma"<<std::endl;
1481  factor
1482  = ecl_mcx3_corr(Eno, Rno, gamma_energy, gamma_cos);
1483  } else { // Exp. data.
1484  //Original definition
1485  //factor
1486  // = ecl_adhoc_corr( Eno, Rno,
1487  // gamma_energy, gamma_cos);
1488  //Modified 20081222
1489  factor
1490  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1491  gamma_energy, gamma_cos);
1492  // Special treatment of Exp.15 processed by b20020125.
1493  if (Eno == 15 && i_previous == 1) {
1494  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1495  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor=" << factor << " " ;
1496  // Original scaling is done by Exp.13 curve.
1497  //Original definition.
1498  //factor13
1499  // = ecl_adhoc_corr( 13, 1,
1500  // gamma_energy, gamma_cos);
1501  //Modified 20081222
1502  factor13
1503  = ecl_adhoc_corr(13, 1, ecl_threshold,
1504  gamma_energy, gamma_cos);
1505  factor = factor / factor13;
1506  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "gamma factor(rel)=" << factor << " " << std::endl;
1507  }
1508  }
1509  break;
1510  case 1:
1511  if (expmc == 2) {
1512  factor = 1.0 / mpi0pdg(gamma_3v.mag());
1513 
1514  } else {
1515  //Original definition.
1516  //factor
1517  // = ecl_adhoc_corr( Eno, Rno,
1518  // gamma_3v.mag(), gamma_3v.cosTheta());
1519  //Modified 20081222
1520  factor
1521  = ecl_adhoc_corr(Eno, Rno, ecl_threshold,
1522  gamma_3v.mag(), gamma_3v.cosTheta());
1523  // Special treatment of Exp.15 processed by b20020125.
1524  if (Eno == 15 && i_previous == 1) {
1525  // dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "Exp.15 version=1 applied case!" << std::endl;
1526  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "factor=" << factor << " " ;
1527  // Original scaling is done by Exp.13 curve.
1528  //Original definition.
1529  //factor13
1530  // = ecl_adhoc_corr( 13, 1,
1531  // gamma_energy, gamma_cos);
1532  //Modified 20081222
1533  factor13
1534  = ecl_adhoc_corr(13, 1, ecl_threshold,
1535  gamma_energy, gamma_cos);
1536  factor = factor / factor13;
1537  //dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "gamma factor(rel)=" << factor << " " << std::endl;
1538  }
1539  factor = factor / mpi0pdg(gamma_3v.mag());
1540  }
1541  break;
1542  default:
1543  factor = 1.0;
1544  }
1545 
1546  // factor should be used as a division.
1547  gamma.px(gamma.px() / factor);
1548  gamma.py(gamma.py() / factor);
1549  gamma.pz(gamma.pz() / factor);
1550  }
1551  }
1552 
1553 
1554 //=====================================================================
1555 //***** make_pi0.cc ********** created 2001/Jul./21st *****
1556 // Create Belle::Mdst_pi0 from Belle::Mdst_gamma and Belle::Mdst_ecl to let
1557 // people get mass-constraint fitted momentum of pi0
1558 // after ad_hoc correction.
1559 // input : option = 0; same as the existent Rececl_pi0
1560 // low_limit and up_limit are ignored.
1561 // option = 1; users can modify mass window as they
1562 // like. The boundary of window is defined
1563 // by low_limit ans up_limit (in GeV).
1564 // option = 2; users can modify mass window in the manner of
1565 // -Xsigma < Mgg - Mpi0 < +Xsigma. The value of
1566 // sigma(mass resolution) is calculated by pi0resol
1567 // function which is originally coded by Hayashii-san.
1568 //*********************************************************
1569  void B2BIIFixMdstModule::make_pi0(int option, double low_limit, double up_limit)
1570  {
1571  //---------------------
1572  // Check Exp. number and Data/MC.
1573  //---------------------
1574  Belle::Belle_event_Manager& evtmgr = Belle::Belle_event_Manager::get_manager();
1575  if (0 == evtmgr.count()) {
1576  // Do nothing if not exist, because we can not know whether
1577  // this event is Exp. data or MC.
1578  return;
1579  }
1580  //Check data or mc. Data:expmc=1, MC:expmc=2.
1581  int expmc = evtmgr[0].ExpMC();
1582  // Check Exp. number.
1583  int Eno = evtmgr[0].ExpNo();
1584  // Set mcdata frag to be compatible with pi0resol func.
1585  bool mcdata = 0;
1586  if (expmc == 2) {
1587  mcdata = 1;
1588  }
1589 
1590  // pi0 mass of PDG.
1591  const double mpi0_pdg = 0.1349739;
1592 
1593  // Default mass window.
1594  // const double low_default = 0.1178;
1595  // const double up_default = 0.1502;
1596  const double low_default = 0.080; // modified 20040606
1597  const double up_default = 0.180; // modified 20040606
1598 
1599  // Maximum iteration of fit.
1600  const int iter_max = 5;
1601 
1602  // Check whether proper option and mass window are given.
1603  switch (option) {
1604  case 0: // option=0 case, set default mass window.
1605  low_limit = low_default;
1606  up_limit = up_default;
1607  break;
1608  case 1: // option=1 case, check given mass window.
1609  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
1610  // If mass window is not correct, do nothing.
1611  B2ERROR("Invalid mass window between ");
1612  B2ERROR(" and " << up_limit);
1613  return;
1614  }
1615  break;
1616  case 2: // pi0 cand. are selected by -Xsigma < Mgg-Mpi0 < +Xsigma.
1617  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl") << "option=2 was selected." << std::endl;
1618  //dbg dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << "mass window is " << low_limit << " & ";
1619  //dbg dout(Debugout::INFO,"B2BIIFixBelle::Mdst_ecl") << up_limit << std::endl;
1620  if (0.0 < low_limit || up_limit < 0.0) {
1621  B2ERROR("option=2 was selected. ");
1622  B2ERROR("Invalid mass window! " << low_limit);
1623  B2ERROR(" should be negative, and " << up_limit);
1624  B2ERROR(" should be positive.");
1625  return;
1626  }
1627  break;
1628  default: // Otherwise, invalid option, do nothing.
1629  B2ERROR("Invalid option=" << option);
1630  return;
1631  }
1632 
1633  // At first, clear already existing Belle::Mdst_pi0.
1634  Belle::Mdst_pi0_Manager::get_manager().remove();
1635 
1636  // Get Belle::Mdst_gamma for photon's momentum
1637  // and get Belle::Mdst_ecl for error matrix.
1638  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1639 
1640  // Re-allocate Belle::Mdst_pi0 table.
1641  Belle::Mdst_pi0_Manager& pi0mgr = Belle::Mdst_pi0_Manager::get_manager();
1642 
1643  // If Only one photon in the event, no need to do anything.
1644  if (gammamgr.count() < 2) {
1645  return;
1646  }
1647 
1648  // Make combination of two Belle::Mdst_gamma.
1649  for (std::vector<Belle::Mdst_gamma>::iterator itgamma = gammamgr.begin();
1650  itgamma != gammamgr.end(); ++itgamma) {
1651  Belle::Mdst_gamma& gamma1 = *itgamma;
1652  CLHEP::Hep3Vector gamma1_3v(gamma1.px(), gamma1.py(), gamma1.pz());
1653  CLHEP::HepLorentzVector gamma1_lv(gamma1_3v, gamma1_3v.mag());
1654 
1655  Belle::Mdst_ecl& ecl1 = gamma1.ecl();
1656 
1657  for (std::vector<Belle::Mdst_gamma>::iterator jtgamma = itgamma + 1;
1658  jtgamma != gammamgr.end(); ++jtgamma) {
1659  Belle::Mdst_gamma& gamma2 = *jtgamma;
1660  CLHEP::Hep3Vector gamma2_3v(gamma2.px(), gamma2.py(), gamma2.pz());
1661  CLHEP::HepLorentzVector gamma2_lv(gamma2_3v, gamma2_3v.mag());
1662 
1663  Belle::Mdst_ecl& ecl2 = gamma2.ecl();
1664 
1665  // Invariant mass before const. fit.
1666  CLHEP::HepLorentzVector gamgam_lv = gamma1_lv + gamma2_lv;
1667  const double mass_before = gamgam_lv.mag();
1668 
1669  // In the case of option=0 or 1, criteria is controlled
1670  // by the inv. mass.
1671  double mass_ctrl{0};
1672  if (option == 0 || option == 1) {
1673  mass_ctrl = mass_before;
1674  }
1675  if (option == 2) {
1676  // Asymmetric lineshape is taken into account.
1677  if (mass_before * mpi0pdg(mass_before) < mpi0_pdg) {
1678  mass_ctrl = (mass_before * mpi0pdg(mass_before) - mpi0_pdg) /
1679  pi0resol(gamgam_lv.vect().mag(),
1680  gamgam_lv.theta() * 180. / M_PI, "lower",
1681  mcdata, Eno, option);
1682  } else {
1683  mass_ctrl = (mass_before * mpi0pdg(mass_before) - mpi0_pdg) /
1684  pi0resol(gamgam_lv.vect().mag(),
1685  gamgam_lv.theta() * 180. / M_PI, "higher",
1686  mcdata, Eno, option);
1687  }
1688  }
1689 
1690  // If invariant mass is inside the window,
1691  // if( low_limit < mass_before && mass_before < up_limit ) // old.
1692  if (low_limit < mass_ctrl && mass_ctrl < up_limit) {
1693  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<"mass="<<mass_before;
1694  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" p="<<gamgam_lv.vect().mag()<<" theta=";
1695  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<gamgam_lv.theta()<<" mcdata="<<mcdata;
1696  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" higher="<<
1697  // pi0resol(gamgam_lv.vect().mag(),
1698  // gamgam_lv.theta()*180./M_PI,
1699  // "higher",mcdata, Eno, option )<<" or ";
1700  // dout(Debugout::DDEBUG,"B2BIIFixBelle::Mdst_ecl")<<" lower="<<
1701  // pi0resol(gamgam_lv.vect().mag(),
1702  // gamgam_lv.theta()*180./M_PI,
1703  // "lower",mcdata, Eno, option )<<std::endl;
1704  // Error matrix(covariant matrix)
1705  CLHEP::HepMatrix V(6, 6, 0);
1706 
1707  V[0][0] = ecl1.error(0);
1708  V[0][1] = V[1][0] = ecl1.error(1);
1709  V[1][1] = ecl1.error(2);
1710  V[0][2] = V[2][0] = ecl1.error(3);
1711  V[1][2] = V[2][1] = ecl1.error(4);
1712  V[2][2] = ecl1.error(5);
1713  V[3][3] = ecl2.error(0);
1714  V[3][4] = V[4][3] = ecl2.error(1);
1715  V[4][4] = ecl2.error(2);
1716  V[3][5] = V[5][3] = ecl2.error(3);
1717  V[4][5] = V[5][4] = ecl2.error(4);
1718  V[5][5] = ecl2.error(5);
1719 
1720  // Measurements; i.e. initial parameters.
1721  CLHEP::HepMatrix y0(6, 1);
1722  y0[0][0] = ecl1.energy();
1723  y0[1][0] = ecl1.phi();
1724  y0[2][0] = ecl1.theta();
1725  y0[3][0] = ecl2.energy();
1726  y0[4][0] = ecl2.phi();
1727  y0[5][0] = ecl2.theta();
1728 
1729  // Copy them to proper matrix which is given to fit.
1730  CLHEP::HepMatrix y(y0);
1731  // Delivative.
1732  CLHEP::HepMatrix Dy(6, 1, 0);
1733 
1734  int iter = 0;
1735  double f_old = DBL_MAX;
1736  double chi2_old = DBL_MAX;
1737  double /*mass_gg,*/ chi2 = DBL_MAX;
1738  bool exit_flag = false;
1739 
1740  // Set parameters to decide whether converged.
1741  const double Df_limit = 0.1;
1742  const double Dchi2_limit = 0.1;
1743  // Do mass constraint fit; iterate until convergence.
1744  while (1) {
1745  const double& E1 = y[0][0];
1746  const double& E2 = y[3][0];
1747  const double sin_theta1 = std::sin(y[2][0]);
1748  const double cos_theta1 = std::cos(y[2][0]);
1749  const double sin_theta2 = std::sin(y[5][0]);
1750  const double cos_theta2 = std::cos(y[5][0]);
1751  const double Dphi = y[1][0] - y[4][0];
1752  const double cos_Dphi = std::cos(Dphi);
1753  const double sin_Dphi = std::sin(Dphi);
1754  const double open_angle = sin_theta1 * sin_theta2 * cos_Dphi
1755  + cos_theta1 * cos_theta2;
1756  const double mass2_gg = 2 * E1 * E2 * (1 - open_angle);
1757  //mass_gg = (mass2_gg > 0) ? std::sqrt(mass2_gg) : -std::sqrt(-mass2_gg);
1758 
1759  // No more iteration. Break to return.
1760  if (exit_flag || ++iter > iter_max)
1761  break;
1762 
1763  // constraint
1764  CLHEP::HepMatrix f(1, 1);
1765  f[0][0] = mass2_gg - (mpi0_pdg * mpi0_pdg);
1766 
1767  // dG/dq_i, G = (M_gg - M_pi0) = 0 is the constraint.
1768  CLHEP::HepMatrix B(1, 6);
1769  B[0][0] = mass2_gg / E1;
1770  B[0][1] = 2 * E1 * E2 * sin_theta1 * sin_theta2 * sin_Dphi;
1771  B[0][2] = 2 * E1 * E2 * (-cos_theta1 * sin_theta2 * cos_Dphi
1772  + sin_theta1 * cos_theta2);
1773  B[0][3] = mass2_gg / E2;
1774  B[0][4] = -B[0][1];
1775  B[0][5] = 2 * E1 * E2 * (-sin_theta1 * cos_theta2 * cos_Dphi
1776  + cos_theta1 * sin_theta2);
1777 
1778  const double sigma2_mass2_gg = (B * V * B.T())[0][0];
1779 
1780  // Proceed one step to get newer value y.
1781  Dy = V * B.T() * (B * Dy - f) / sigma2_mass2_gg;
1782  y = y0 + Dy;
1783  int ierr;
1784  chi2 = (Dy.T() * V.inverse(ierr) * Dy)[0][0];
1785  double Dchi2 = fabs(chi2 - chi2_old);
1786  chi2_old = chi2;
1787  double Df = fabs(f[0][0] - f_old);
1788  f_old = f[0][0];
1789 
1790  // When chi-sq. change is small enough and mass is
1791  if (Dchi2 < Dchi2_limit && Df < Df_limit)
1792  exit_flag = true;;
1793 
1794  }
1795  // Sort fitted result into proper variables.
1796  double pi0_E = y[0][0] + y[3][0];
1797  double pi0_px = y[0][0] * std::cos(y[1][0]) * std::sin(y[2][0])
1798  + y[3][0] * std::cos(y[4][0]) * std::sin(y[5][0]);
1799  double pi0_py = y[0][0] * std::sin(y[1][0]) * std::sin(y[2][0])
1800  + y[3][0] * std::sin(y[4][0]) * std::sin(y[5][0]);
1801  double pi0_pz = y[0][0] * std::cos(y[2][0]) + y[3][0] * std::cos(y[5][0]);
1802  double pi0_mass = mass_before;
1803 
1804  // m_chi2 = chi2; // protect...
1805  double pi0_chi2 = (chi2 > FLT_MAX) ? FLT_MAX : chi2;
1806 
1807  // Fill Belle::Mdst_pi0 based on the fit result.
1808  Belle::Mdst_pi0& pi0 = pi0mgr.add();
1809  pi0.gamma(0, gamma1);
1810  pi0.gamma(1, gamma2);
1811  pi0.px(pi0_px);
1812  pi0.py(pi0_py);
1813  pi0.pz(pi0_pz);
1814  pi0.energy(pi0_E);
1815  pi0.mass(pi0_mass);
1816  pi0.chisq(pi0_chi2);
1817  }
1818  }
1819  }
1820  }
1821 
1822  void B2BIIFixMdstModule::make_pi0_primary_vertex(int option, double low_limit, double up_limit,
1823  const HepPoint3D& /*epvtx*/,
1824  const CLHEP::HepSymMatrix& epvtx_err)
1825  {
1826 #if 0
1827  // pi0 mass of PDG.
1828  const double mpi0_pdg = 0.1349739;
1829 
1830  // Default mass window.
1831  const double low_default = 0.1178;
1832  const double up_default = 0.1502;
1833 
1834  // Maximum iteration of fit.
1835  const int iter_max = 5;
1836 
1837  // Check whether proper option and mass window are given.
1838  switch (option) {
1839  case 0: // option=0 case, set default mass window.
1840  low_limit = low_default;
1841  up_limit = up_default;
1842  break;
1843  case 1: // option=1 case, check given mass window.
1844  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
1845  // If mass window is not correct, do nothing.
1846  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << "Invalid mass window between " << low_limit;
1847  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << " and " << up_limit << std::endl;
1848  return;
1849  }
1850  break;
1851  default: // Otherwise, invalid option, do nothing.
1852  dout(Debugout::ERR, "B2BIIFixBelle::Mdst_ecl") << "Invalid option=" << option << std::endl;
1853  return;
1854  }
1855 
1856  // At first, clear already existing Belle::Mdst_pi0.
1857  Belle::Mdst_pi0_Manager::get_manager().remove();
1858 
1859  // Get Belle::Mdst_gamma for photon's momentum
1860  // and get Belle::Mdst_ecl for error matrix.
1861  Belle::Mdst_gamma_Manager& gammamgr = Belle::Mdst_gamma_Manager::get_manager();
1862 
1863  // Re-allocate Belle::Mdst_pi0 table.
1864  Belle::Mdst_pi0_Manager& pi0mgr = Belle::Mdst_pi0_Manager::get_manager();
1865 
1866  // If Only one photon in the event, no need to do anything.
1867  if (gammamgr.count() < 2) {
1868  return;
1869  }
1870 
1871  // Make combination of two Belle::Mdst_gamma.
1872  for (std::vector<Belle::Mdst_gamma>::iterator itgamma = gammamgr.begin();
1873  itgamma != gammamgr.end(); ++itgamma) {
1874  Belle::Mdst_gamma& gamma1 = *itgamma;
1875  CLHEP::Hep3Vector gamma1_3v(gamma1.px(), gamma1.py(), gamma1.pz());
1876  CLHEP::HepLorentzVector gamma1_lv(gamma1_3v, gamma1_3v.mag());
1877 
1878  Belle::Mdst_ecl& ecl1 = gamma1.ecl();
1879 
1880  const double r_i = ecl1.r();
1881  const double dzr_i = std::sqrt(epvtx_err[2][2]) / r_i;
1882 
1883  const double theta_i0 = ecl1.theta();
1884  const double sin_th_i0 = std::sin(theta_i0);
1885 
1886  for (std::vector<Belle::Mdst_gamma>::iterator jtgamma = itgamma + 1;
1887  jtgamma != gammamgr.end(); ++jtgamma) {
1888  Belle::Mdst_gamma& gamma2 = *jtgamma;
1889  CLHEP::Hep3Vector gamma2_3v(gamma2.px(), gamma2.py(), gamma2.pz());
1890  CLHEP::HepLorentzVector gamma2_lv(gamma2_3v, gamma2_3v.mag());
1891 
1892  Belle::Mdst_ecl& ecl2 = gamma2.ecl();
1893 
1894  // Invariant mass before const. fit.
1895  CLHEP::HepLorentzVector gamgam_lv = gamma1_lv + gamma2_lv;
1896  const double mass_before = gamgam_lv.mag();
1897 
1898  // If invariant mass is inside the window,
1899  if (low_limit < mass_before && mass_before < up_limit) {
1900  // Error matrix(covariant matrix)
1901  CLHEP::HepMatrix V(6, 6, 0);
1902 
1903  V[0][0] = ecl1.error(0);
1904  V[0][1] = V[1][0] = ecl1.error(1);
1905  V[1][1] = ecl1.error(2);
1906  V[0][2] = V[2][0] = ecl1.error(3);
1907  V[1][2] = V[2][1] = ecl1.error(4);
1908  V[2][2] = ecl1.error(5);
1909  V[3][3] = ecl2.error(0);
1910  V[3][4] = V[4][3] = ecl2.error(1);
1911  V[4][4] = ecl2.error(2);
1912  V[3][5] = V[5][3] = ecl2.error(3);
1913  V[4][5] = V[5][4] = ecl2.error(4);
1914  V[5][5] = ecl2.error(5);
1915 
1916  // Correlation term
1917  const double r_j = ecl2.r();
1918  const double dzr_j = sqrt(epvtx_err[2][2]) / r_j;
1919 
1920  const double theta_j0 = ecl2.theta();
1921  const double sin_th_j0 = std::sin(theta_j0);
1922 
1923  V[2][5] = V[5][2] = dzr_i * sin_th_i0 * dzr_j * sin_th_j0;
1924  // Measurements; i.e. initial parameters.
1925  CLHEP::HepMatrix y0(6, 1);
1926  y0[0][0] = ecl1.energy();
1927  y0[1][0] = ecl1.phi();
1928  y0[2][0] = ecl1.theta();
1929  y0[3][0] = ecl2.energy();
1930  y0[4][0] = ecl2.phi();
1931  y0[5][0] = ecl2.theta();
1932 
1933  // Copy them to proper matrix which is given to fit.
1934  CLHEP::HepMatrix y(y0);
1935  // Delivative.
1936  CLHEP::HepMatrix Dy(6, 1, 0);
1937 
1938  int iter = 0;
1939  double Df, f_old = DBL_MAX;
1940  double Dchi2, chi2_old = DBL_MAX;
1941  double mass_gg, chi2 = DBL_MAX;
1942  bool exit_flag = false;
1943 
1944  // Set parameters to decide whether converged.
1945  const double Df_limit = 0.1;
1946  const double Dchi2_limit = 0.1;
1947  // Do mass constraint fit; iterate until convergence.
1948  while (1) {
1949  const double& E1 = y[0][0];
1950  const double& E2 = y[3][0];
1951  const double sin_theta1 = std::sin(y[2][0]);
1952  const double cos_theta1 = std::cos(y[2][0]);
1953  const double sin_theta2 = std::sin(y[5][0]);
1954  const double cos_theta2 = std::cos(y[5][0]);
1955  const double Dphi = y[1][0] - y[4][0];
1956  const double cos_Dphi = std::cos(Dphi);
1957  const double sin_Dphi = std::sin(Dphi);
1958  const double open_angle = sin_theta1 * sin_theta2 * cos_Dphi
1959  + cos_theta1 * cos_theta2;
1960  const double mass2_gg = 2 * E1 * E2 * (1 - open_angle);
1961  mass_gg = (mass2_gg > 0) ? std::sqrt(mass2_gg) : -std::sqrt(-mass2_gg);
1962 
1963  // No more iteration. Break to return.
1964  if (exit_flag || ++iter > iter_max)
1965  break;
1966 
1967  // constraint
1968  CLHEP::HepMatrix f(1, 1);
1969  f[0][0] = mass2_gg - (mpi0_pdg * mpi0_pdg);
1970 
1971  // dG/dq_i, G = (M_gg - M_pi0) = 0 is the constraint.
1972  CLHEP::HepMatrix B(1, 6);
1973  B[0][0] = mass2_gg / E1;
1974  B[0][1] = 2 * E1 * E2 * sin_theta1 * sin_theta2 * sin_Dphi;
1975  B[0][2] = 2 * E1 * E2 * (-cos_theta1 * sin_theta2 * cos_Dphi
1976  + sin_theta1 * cos_theta2);
1977  B[0][3] = mass2_gg / E2;
1978  B[0][4] = -B[0][1];
1979  B[0][5] = 2 * E1 * E2 * (-sin_theta1 * cos_theta2 * cos_Dphi
1980  + cos_theta1 * sin_theta2);
1981 
1982  const double sigma2_mass2_gg = (B * V * B.T())[0][0];
1983 
1984  // Proceed one step to get newer value y.
1985  Dy = V * B.T() * (B * Dy - f) / sigma2_mass2_gg;
1986  y = y0 + Dy;
1987  int ierr;
1988  chi2 = (Dy.T() * V.inverse(ierr) * Dy)[0][0];
1989  Dchi2 = fabs(chi2 - chi2_old);
1990  chi2_old = chi2;
1991  Df = fabs(f[0][0] - f_old);
1992  f_old = f[0][0];
1993 
1994  // When chi-sq. change is small enough and mass is
1995  if (Dchi2 < Dchi2_limit && Df < Df_limit)
1996  exit_flag = true;;
1997 
1998  }
1999  // Sort fitted result into proper variables.
2000  double pi0_E = y[0][0] + y[3][0];
2001  double pi0_px = y[0][0] * std::cos(y[1][0]) * std::sin(y[2][0])
2002  + y[3][0] * std::cos(y[4][0]) * std::sin(y[5][0]);
2003  double pi0_py = y[0][0] * std::sin(y[1][0]) * std::sin(y[2][0])
2004  + y[3][0] * std::sin(y[4][0]) * std::sin(y[5][0]);
2005  double pi0_pz = y[0][0] * std::cos(y[2][0]) + y[3][0] * std::cos(y[5][0]);
2006  double pi0_mass = mass_before;
2007 
2008  // m_chi2 = chi2; // protect...
2009  double pi0_chi2 = (chi2 > FLT_MAX) ? FLT_MAX : chi2;
2010 
2011  // Fill Belle::Mdst_pi0 based on the fit result.
2012  Belle::Mdst_pi0& pi0 = pi0mgr.add();
2013  pi0.gamma(0, gamma1);
2014  pi0.gamma(1, gamma2);
2015  pi0.px(pi0_px);
2016  pi0.py(pi0_py);
2017  pi0.pz(pi0_pz);
2018  pi0.energy(pi0_E);
2019  pi0.mass(pi0_mass);
2020  pi0.chisq(pi0_chi2);
2021  }
2022  }
2023  }
2024 #else
2025  // pi0 mass of PDG. ;
2026  const double mpi0_pdg = 0.1349739;
2027  const double m2pi0_pdg = mpi0_pdg * mpi0_pdg;
2028  // Default mass window.;
2029  const double low_default = 0.1178;
2030  const double up_default = 0.1502;
2031  // Maximum iteration of fit.;
2032  const int iter_max = 5;
2033 
2034  // Check whether proper option and mass window are given.;
2035  switch (option) {
2036  case 0: // option=0 case, set default mass window.;
2037  low_limit = low_default;
2038  up_limit = up_default;
2039  break;
2040  case 1: // option=1 case, check given mass window.;
2041  if (mpi0_pdg < low_limit || up_limit < mpi0_pdg) {
2042  // If mass window is not correct, do nothing.;
2043  B2ERROR("Invalid mass window between " << low_limit);
2044  B2ERROR(" and " << up_limit);
2045  return;
2046  }
2047  break;
2048  default: // Otherwise, invalid option, do nothing.;
2049  B2ERROR("Invalid option=" << option);
2050  return;
2051  }
2052  Belle::Mdst_pi0_Manager& pi0_mgr = Belle::Mdst_pi0_Manager::get_manager();
2053  Belle::Mdst_gamma_Manager& gamma_mgr = Belle::Mdst_gamma_Manager::get_manager();
2054 
2055  // At first, clear already existing Belle::Mdst_pi0.;
2056  Belle::Mdst_pi0_Manager::get_manager().remove();
2057  // If Only one photon in the event, no need to do anything.;
2058  if (gamma_mgr.count() < 2) {
2059  return;
2060  }
2061 
2062  for (std::vector<Belle::Mdst_gamma>::iterator i = gamma_mgr.begin();
2063  i != gamma_mgr.end(); ++i) {
2064  const Belle::Mdst_gamma& gamma_i = *i;
2065  if (!gamma_i.ecl()) {
2066  continue;
2067  }
2068  const Belle::Mdst_ecl& ecl_i = gamma_i.ecl();
2069  const double r_i = ecl_i.r();
2070  const double e_i0 = ecl_i.energy();
2071  const double phi_i0 = ecl_i.phi();
2072  const double theta_i0 = ecl_i.theta();
2073 
2074  const double sin_th_i0 = std::sin(theta_i0);
2075  const double cos_th_i0 = std::cos(theta_i0);
2076 
2077  CLHEP::HepSymMatrix err_i(3, 0);
2078  err_i[0][0] = ecl_i.error(0);
2079  err_i[1][0] = ecl_i.error(1); err_i[1][1] = ecl_i.error(2);
2080  err_i[2][0] = ecl_i.error(3); err_i[2][1] = ecl_i.error(4); err_i[2][2] = ecl_i.error(5);
2081 
2082  const double dzr_i = std::sqrt(epvtx_err[2][2]) / r_i;
2083 
2084  for (std::vector<Belle::Mdst_gamma>::iterator j = i + 1; j != gamma_mgr.end(); ++j) {
2085  const Belle::Mdst_gamma& gamma_j = *j;
2086  if (!gamma_j.ecl()) {
2087  continue;
2088  }
2089  const Belle::Mdst_ecl& ecl_j = gamma_j.ecl();
2090  const double r_j = ecl_j.r();
2091  const double e_j0 = ecl_j.energy();
2092  const double phi_j0 = ecl_j.phi();
2093  const double theta_j0 = ecl_j.theta();
2094 
2095  const double sin_th_j0 = std::sin(theta_j0);
2096  const double cos_th_j0 = std::cos(theta_j0);
2097 
2098  CLHEP::HepSymMatrix err_j(3, 0);
2099  err_j[0][0] = ecl_j.error(0);
2100  err_j[1][0] = ecl_j.error(1); err_j[1][1] = ecl_j.error(2);
2101  err_j[2][0] = ecl_j.error(3); err_j[2][1] = ecl_j.error(4); err_j[2][2] = ecl_j.error(5);
2102 
2103  const double dzr_j = sqrt(epvtx_err[2][2]) / r_j;
2104 
2105  // -------------------------------------------------------------------- ;
2106 
2107  double dth_i_x_dth_j = dzr_i * sin_th_i0 * dzr_j * sin_th_j0;
2108 
2109  double det = err_i[2][2] * err_j[2][2] - dth_i_x_dth_j * dth_i_x_dth_j;
2110 
2111 
2112  double e_i = e_i0, phi_i = phi_i0, sin_th_i = sin_th_i0, cos_th_i = cos_th_i0;
2113  double e_j = e_j0, phi_j = phi_j0, sin_th_j = sin_th_j0, cos_th_j = cos_th_j0;
2114 
2115 
2116  double dphi = phi_i - phi_j;
2117  double sin_dphi = std::sin(dphi);
2118  double cos_dphi = std::cos(dphi);
2119 
2120  // Cos of opening angle ;
2121  double cos_th_ij = sin_th_i * sin_th_j * cos_dphi + cos_th_i * cos_th_j;
2122  double m2_ij = 2.0 * e_i * e_j * (1.0 - cos_th_ij);
2123  double mass = std::sqrt(m2_ij);
2124 
2125  if (mass < low_limit || mass > up_limit)
2126  continue;
2127 
2128  double de_i = 0.0, dphi_i = 0.0, dtheta_i = 0.0;
2129  double de_j = 0.0, dphi_j = 0.0, dtheta_j = 0.0;
2130 
2131  int it = 0;
2132 
2133  double f = 1.0e+6;
2134  double chisq = 1.0e+6;
2135  double df = 0.0, dchisq = 0.0;
2136 
2137  const double del_chisq = 0.1;
2138  const double del_f = 0.1;
2139  //const int icovar = 1;
2140 
2141  double mass0 = 0.0;
2142  double dmass0 = 0.0;
2143  //double dmass = 0.0;
2144  //double sdev = 0.0;
2145 
2146  do {
2147 
2148  const double dcos_th_ij_dphi_i = -sin_th_i * sin_th_j * sin_dphi;
2149  const double dcos_th_ij_dphi_j = sin_th_i * sin_th_j * sin_dphi;
2150  const double dcos_th_ij_dth_i = cos_th_i * sin_th_j * cos_dphi - sin_th_i * cos_th_j;
2151  const double dcos_th_ij_dth_j = cos_th_j * sin_th_i * cos_dphi - sin_th_j * cos_th_i;
2152 
2153  const double dm2_de_i = m2_ij / e_i;
2154  const double dm2_de_j = m2_ij / e_j;
2155 
2156  const double dm2_dcos_th_ij = -2.0 * e_i * e_j;
2157 
2158  const double dm2_dphi_i = dm2_dcos_th_ij * dcos_th_ij_dphi_i;
2159  const double dm2_dphi_j = dm2_dcos_th_ij * dcos_th_ij_dphi_j;
2160  const double dm2_dth_i = dm2_dcos_th_ij * dcos_th_ij_dth_i;
2161  const double dm2_dth_j = dm2_dcos_th_ij * dcos_th_ij_dth_j;
2162 
2163  double dm2_2
2164  = dm2_de_i * dm2_de_i * err_i[0][0]
2165  + dm2_de_j * dm2_de_j * err_j[0][0]
2166  + dm2_dphi_i * dm2_dphi_i * err_i[1][0]
2167  + dm2_dphi_j * dm2_dphi_j * err_j[1][0]
2168  + dm2_dth_i * dm2_dth_i * err_i[2][0]
2169  + dm2_dth_j * dm2_dth_j * err_j[2][0]
2170  + 2.0 * dm2_dth_i * dm2_dth_j * dth_i_x_dth_j;
2171 
2172  if (dm2_2 < 0.0) {
2173  dm2_2 = 0.0;
2174  }
2175  if (it++ == 0) {
2176  mass0 = mass;
2177  dmass0 = std::sqrt(dm2_2);
2178  /*if (mass0 > 0.0) {
2179  dmass = 0.5 * dmass0 / mass;
2180  }*/
2181  }
2182  //const double residual = mass - mpi0_pdg;
2183  //const double pull = residual / dmass;
2184  if (it >= iter_max ||
2185  mass < low_default || mass > up_default) {
2186  //sdev = (dmass > 0.0) ? pull : -100;
2187  break;
2188  }
2189 
2190  const double del_m2
2191  = (m2pi0_pdg - m2_ij
2192  + dm2_de_i * de_i
2193  + dm2_de_j * de_j
2194  + dm2_dphi_i * dphi_i
2195  + dm2_dphi_j * dphi_j
2196  + dm2_dth_i * dtheta_i
2197  + dm2_dth_j * dtheta_j) / dm2_2;
2198 
2199  de_i = del_m2 * dm2_de_i * err_i[0][0];
2200  de_j = del_m2 * dm2_de_j * err_j[0][0];
2201  dphi_i = del_m2 * dm2_dphi_i * err_i[1][1];
2202  dphi_j = del_m2 * dm2_dphi_j * err_j[1][1];
2203  dtheta_i = del_m2 * (dm2_dth_i * err_i[2][2] + dm2_dth_j * dth_i_x_dth_j);
2204  dtheta_j = del_m2 * (dm2_dth_j * err_j[2][2] + dm2_dth_i * dth_i_x_dth_j);
2205 
2206 
2207  dchisq = chisq;
2208  chisq
2209  = de_i * de_i / err_i[0][0]
2210  + de_j * de_j / err_j[0][0]
2211  + dphi_i * dphi_i / err_i[1][1]
2212  + dphi_j * dphi_j / err_j[1][1]
2213  + (dtheta_i * dtheta_i * err_i[2][2]
2214  + dtheta_j * dtheta_j * err_j[2][2]
2215  + 2.0 * dtheta_i * dtheta_j * dth_i_x_dth_j) / det;
2216  dchisq -= chisq;
2217 
2218 
2219  e_i = e_i0 + de_i;
2220  e_j = e_j0 + de_j;
2221  phi_i = phi_i0 + dphi_i;
2222  phi_j = phi_j0 + dphi_j;
2223  double theta_i = theta_i0 + dtheta_i;
2224  double theta_j = theta_j0 + dtheta_j;
2225 
2226  sin_th_i = std::sin(theta_i);
2227  cos_th_i = std::cos(theta_i);
2228  sin_th_j = std::sin(theta_j);
2229  cos_th_j = std::cos(theta_j);
2230 
2231  dth_i_x_dth_j = dzr_i * sin_th_i * dzr_j * sin_th_j;
2232  det = err_i[2][2] * err_j[2][2] - dth_i_x_dth_j * dth_i_x_dth_j;
2233 
2234  dphi = phi_i - phi_j;
2235  sin_dphi = std::sin(dphi);
2236  cos_dphi = std::cos(dphi);
2237 
2238  cos_th_ij = sin_th_i * sin_th_j * cos_dphi + cos_th_i * cos_th_j;
2239  m2_ij = 2.0 * e_i * e_j * (1.0 - cos_th_ij);
2240  mass = std::sqrt(m2_ij);
2241 
2242  df = f;
2243  f = fabs(m2pi0_pdg - m2_ij) / dmass0;
2244  df -= f;
2245  ++it;
2246  } while (std::fabs(df) > del_f || std::fabs(dchisq) > del_chisq);
2247  const double cos_phi_i = std::cos(phi_i);
2248  const double cos_phi_j = std::cos(phi_j);
2249  const double sin_phi_i = std::sin(phi_i);
2250  const double sin_phi_j = std::sin(phi_j);
2251 
2252  const CLHEP::HepLorentzVector p4_i(e_i * sin_th_i * cos_phi_i,
2253  e_i * sin_th_i * sin_phi_i, e_i * cos_th_i, e_i);
2254  const CLHEP::HepLorentzVector p4_j(e_j * sin_th_j * cos_phi_j,
2255  e_j * sin_th_j * sin_phi_j, e_j * cos_th_j, e_j);
2256 
2257  const CLHEP::HepLorentzVector p4_pi0(p4_i + p4_j);
2258 
2259  // Fill Belle::Mdst_pi0 based on the fit result.;
2260  Belle::Mdst_pi0& pi0 = pi0_mgr.add();
2261  pi0.gamma(0, gamma_i);
2262  pi0.gamma(1, gamma_j);
2263  pi0.px(p4_pi0.x());
2264  pi0.py(p4_pi0.y());
2265  pi0.pz(p4_pi0.z());
2266  pi0.energy(p4_pi0.e());
2267  pi0.mass(mass0);
2268  pi0.chisq(chisq);
2269  }
2270  }
2271 #endif
2272  return;
2273  }
2274 
2275  void B2BIIFixMdstModule::correct_ecl_primary_vertex(const HepPoint3D& epvtx,
2276  const CLHEP::HepSymMatrix& epvtx_err)
2277  {
2278 
2279  Belle::Mdst_gamma_Manager& Gamma = Belle::Mdst_gamma_Manager::get_manager();
2280 
2281  for (std::vector<Belle::Mdst_gamma>::iterator
2282  it = Gamma.begin(); it != Gamma.end(); ++it) {
2283  double r(it->ecl().r());
2284  double theta(it->ecl().theta());
2285  double phi(it->ecl().phi());
2286  // double stheta(std::sqrt(it->ecl().error(2)));
2287  // double sphi(std::sqrt(it->ecl().error(5)));
2288  double st(std::sin(theta));
2289  double ct(std::cos(theta));
2290  double sp(std::sin(phi));
2291  double cp(std::cos(phi));
2292  HepPoint3D gamma_pos(r * st * cp, r * st * sp, r * ct);
2293  CLHEP::Hep3Vector gamma_vec(gamma_pos - epvtx);
2294  double hsq(gamma_vec.perp2());
2295  double rsq(gamma_vec.mag2());
2296  double stheta_sq_new = it->ecl().error(5) + epvtx_err(3, 3) * (hsq / (rsq * rsq));
2297  CLHEP::Hep3Vector gamma_dir(gamma_vec.unit());
2298  double e(it->ecl().energy());
2299  it->px(e * gamma_dir.x());
2300  it->py(e * gamma_dir.y());
2301  it->pz(e * gamma_dir.z());
2302  it->ecl().theta(gamma_dir.theta());
2303  it->ecl().r(gamma_vec.mag());
2304  it->ecl().error(5, stheta_sq_new);
2305  }
2306  }
2307 
2308 
2310 } // namespace Belle
Belle2::B2BIIFixMdstModule::m_reprocess_version
int m_reprocess_version
Reprocess version (=0:old; =1:new)
Definition: B2BIIFixMdstModule.h:181
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
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:1569
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:1822
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:2275
Belle2::B2BIIFixMdstModule::correct_ecl
void correct_ecl(int, int)
Correct photon's momenta and error matrix.
Definition: B2BIIFixMdstModule_ecl.cc:1226
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::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17