Belle II Software development
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// update 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 corrections 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
122namespace 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// Corresponding 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 ===== created 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, incorrect 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 and 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 // Derivative.
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 // Derivative.
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:178
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
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::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::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::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::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
void correct_ecl(int, int)
Correct photon's momenta and error matrix.
Definition: B2BIIFixMdstModule_ecl.cc:1226
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::mpi0pdg
static double mpi0pdg(double Energy)
Make MC mass peak to PDG value.
Definition: B2BIIFixMdstModule_ecl.cc:642
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