diff --git a/data/T2K/CC1pip/CH/MomentumPion.root b/data/T2K/CC1pip/CH/MomentumPion.root new file mode 100644 index 0000000..3d52fea Binary files /dev/null and b/data/T2K/CC1pip/CH/MomentumPion.root differ diff --git a/data/T2K/CC1pip/CH/MomentumPion.rootout.root b/data/T2K/CC1pip/CH/MomentumPion.rootout.root new file mode 100644 index 0000000..5b2c080 Binary files /dev/null and b/data/T2K/CC1pip/CH/MomentumPion.rootout.root differ diff --git a/data/T2K/CC1pip/CH/MomentumPion.txt b/data/T2K/CC1pip/CH/MomentumPion.txt new file mode 100644 index 0000000..af65692 --- /dev/null +++ b/data/T2K/CC1pip/CH/MomentumPion.txt @@ -0,0 +1,184 @@ + + Variable : Momentum_pion + +(Low end - Up end ) (GeV) cross-section (10^{-38}cm2/GeV/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (0.0000,0.2000) 0.0322 +Bin 1 (0.2000,0.3000) 0.1133 +Bin 2 (0.3000,0.4000) 0.0979 +Bin 3 (0.4000,0.5000) 0.0552 +Bin 4 (0.5000,0.6000) 0.0329 +Bin 5 (0.6000,0.7000) 0.0182 +Bin 6 (0.7000,0.8000) 0.0164 +Bin 7 (0.8000,0.9000) 0.0146 +Bin 8 (0.9000,1.0000) 0.0107 +Bin 9 (1.0000,1.1000) 0.0103 +Bin 10 (1.1000,1.2000) 0.0102 +Bin 11 (1.2000,1.3000) 0.0082 +Bin 12 (1.3000,1.4000) 0.0050 +Bin 13 (1.4000,1.6000) 0.0043 +Bin 14 (1.6000,2.0000) 0.0021 +Bin 15 (2.0000,3.0000) 0.0013 +Bin 16 (3.0000,15.0000) 0.0001 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.00025742 + ( 0-1 ) 0.00044177 + ( 0-2 ) 0.00021643 + ( 0-3 ) 0.00012093 + ( 0-4 ) 0.00007565 + ( 0-5 ) 0.00004334 + ( 0-6 ) 0.00003825 + ( 0-7 ) 0.00003113 + ( 0-8 ) 0.00002540 + ( 0-9 ) 0.00002340 + ( 0-10 ) 0.00002229 + ( 0-11 ) 0.00001811 + ( 0-12 ) 0.00001323 + ( 0-13 ) 0.00001236 + ( 0-14 ) 0.00000895 + ( 0-15 ) 0.00000450 + ( 0-16 ) 0.00000032 + ( 1-1 ) 0.00109701 + ( 1-2 ) 0.00050478 + ( 1-3 ) 0.00026115 + ( 1-4 ) 0.00016197 + ( 1-5 ) 0.00009815 + ( 1-6 ) 0.00008905 + ( 1-7 ) 0.00007268 + ( 1-8 ) 0.00005803 + ( 1-9 ) 0.00005465 + ( 1-10 ) 0.00005414 + ( 1-11 ) 0.00004486 + ( 1-12 ) 0.00003219 + ( 1-13 ) 0.00003090 + ( 1-14 ) 0.00002213 + ( 1-15 ) 0.00001092 + ( 1-16 ) 0.00000075 + ( 2-2 ) 0.00038075 + ( 2-3 ) 0.00020255 + ( 2-4 ) 0.00012507 + ( 2-5 ) 0.00007649 + ( 2-6 ) 0.00007052 + ( 2-7 ) 0.00005703 + ( 2-8 ) 0.00004560 + ( 2-9 ) 0.00004171 + ( 2-10 ) 0.00004124 + ( 2-11 ) 0.00003390 + ( 2-12 ) 0.00002404 + ( 2-13 ) 0.00002223 + ( 2-14 ) 0.00001484 + ( 2-15 ) 0.00000746 + ( 2-16 ) 0.00000050 + ( 3-3 ) 0.00014639 + ( 3-4 ) 0.00008776 + ( 3-5 ) 0.00005136 + ( 3-6 ) 0.00004617 + ( 3-7 ) 0.00003615 + ( 3-8 ) 0.00002870 + ( 3-9 ) 0.00002560 + ( 3-10 ) 0.00002514 + ( 3-11 ) 0.00002085 + ( 3-12 ) 0.00001431 + ( 3-13 ) 0.00001275 + ( 3-14 ) 0.00000800 + ( 3-15 ) 0.00000409 + ( 3-16 ) 0.00000028 + ( 4-4 ) 0.00007383 + ( 4-5 ) 0.00004221 + ( 4-6 ) 0.00003450 + ( 4-7 ) 0.00002653 + ( 4-8 ) 0.00002096 + ( 4-9 ) 0.00001829 + ( 4-10 ) 0.00001726 + ( 4-11 ) 0.00001435 + ( 4-12 ) 0.00000981 + ( 4-13 ) 0.00000845 + ( 4-14 ) 0.00000478 + ( 4-15 ) 0.00000257 + ( 4-16 ) 0.00000018 + ( 5-5 ) 0.00003302 + ( 5-6 ) 0.00002597 + ( 5-7 ) 0.00001857 + ( 5-8 ) 0.00001449 + ( 5-9 ) 0.00001213 + ( 5-10 ) 0.00001114 + ( 5-11 ) 0.00000905 + ( 5-12 ) 0.00000613 + ( 5-13 ) 0.00000496 + ( 5-14 ) 0.00000258 + ( 5-15 ) 0.00000155 + ( 5-16 ) 0.00000011 + ( 6-6 ) 0.00002747 + ( 6-7 ) 0.00001951 + ( 6-8 ) 0.00001433 + ( 6-9 ) 0.00001219 + ( 6-10 ) 0.00001137 + ( 6-11 ) 0.00000929 + ( 6-12 ) 0.00000628 + ( 6-13 ) 0.00000538 + ( 6-14 ) 0.00000305 + ( 6-15 ) 0.00000166 + ( 6-16 ) 0.00000011 + ( 7-7 ) 0.00001952 + ( 7-8 ) 0.00001352 + ( 7-9 ) 0.00001079 + ( 7-10 ) 0.00000979 + ( 7-11 ) 0.00000801 + ( 7-12 ) 0.00000541 + ( 7-13 ) 0.00000476 + ( 7-14 ) 0.00000269 + ( 7-15 ) 0.00000143 + ( 7-16 ) 0.00000009 + ( 8-8 ) 0.00001162 + ( 8-9 ) 0.00000970 + ( 8-10 ) 0.00000851 + ( 8-11 ) 0.00000682 + ( 8-12 ) 0.00000471 + ( 8-13 ) 0.00000429 + ( 8-14 ) 0.00000252 + ( 8-15 ) 0.00000125 + ( 8-16 ) 0.00000008 + ( 9-9 ) 0.00000999 + ( 9-10 ) 0.00000892 + ( 9-11 ) 0.00000712 + ( 9-12 ) 0.00000495 + ( 9-13 ) 0.00000482 + ( 9-14 ) 0.00000313 + ( 9-15 ) 0.00000140 + ( 9-16 ) 0.00000008 + ( 10-10 ) 0.00000988 + ( 10-11 ) 0.00000818 + ( 10-12 ) 0.00000554 + ( 10-13 ) 0.00000520 + ( 10-14 ) 0.00000341 + ( 10-15 ) 0.00000143 + ( 10-16 ) 0.00000008 + ( 11-11 ) 0.00000771 + ( 11-12 ) 0.00000525 + ( 11-13 ) 0.00000492 + ( 11-14 ) 0.00000328 + ( 11-15 ) 0.00000131 + ( 11-16 ) 0.00000007 + ( 12-12 ) 0.00000428 + ( 12-13 ) 0.00000421 + ( 12-14 ) 0.00000286 + ( 12-15 ) 0.00000108 + ( 12-16 ) 0.00000006 + ( 13-13 ) 0.00000488 + ( 13-14 ) 0.00000351 + ( 13-15 ) 0.00000125 + ( 13-16 ) 0.00000006 + ( 14-14 ) 0.00000305 + ( 14-15 ) 0.00000108 + ( 14-16 ) 0.00000005 + ( 15-15 ) 0.00000051 + ( 15-16 ) 0.00000003 + ( 16-16 ) 0.00000000 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/PmuThetamu.root b/data/T2K/CC1pip/CH/PmuThetamu.root new file mode 100644 index 0000000..aba8a63 Binary files /dev/null and b/data/T2K/CC1pip/CH/PmuThetamu.root differ diff --git a/data/T2K/CC1pip/CH/PmuThetamu.rootout.root b/data/T2K/CC1pip/CH/PmuThetamu.rootout.root new file mode 100644 index 0000000..f60fbd2 Binary files /dev/null and b/data/T2K/CC1pip/CH/PmuThetamu.rootout.root differ diff --git a/data/T2K/CC1pip/CH/PmuThetamu.txt b/data/T2K/CC1pip/CH/PmuThetamu.txt new file mode 100644 index 0000000..f0d18d0 --- /dev/null +++ b/data/T2K/CC1pip/CH/PmuThetamu.txt @@ -0,0 +1,244 @@ + + Variable : p_mu theta_mu + +Cos (Low end-Up end) Mom (Low end-Up end) GeV cross-section 10^{-38}cm2/GeV/Nucleon +--------------------------------------------------------------------------------------- +Bin 0 (0.0000,0.8000) (0.0000,0.4000) 0.1004 +Bin 1 (0.0000,0.8000) (0.4000,1.2000) 0.0233 +Bin 2 (0.0000,0.8000) (1.2000,1.6000) 0.0017 +Bin 3 (0.0000,0.8000) (1.6000,2.0000) 0.0004 +Bin 4 (0.0000,0.8000) (2.0000,15.0000) 0.0000 +Bin 5 (0.8000,0.8500) (0.0000,0.4000) 0.1175 +Bin 6 (0.8000,0.8500) (0.4000,1.2000) 0.1062 +Bin 7 (0.8000,0.8500) (1.2000,1.6000) 0.0182 +Bin 8 (0.8000,0.8500) (1.6000,2.0000) 0.0043 +Bin 9 (0.8000,0.8500) (2.0000,15.0000) 0.0000 +Bin 10 (0.8500,0.9000) (0.0000,0.4000) 0.0998 +Bin 11 (0.8500,0.9000) (0.4000,1.2000) 0.1343 +Bin 12 (0.8500,0.9000) (1.2000,1.6000) 0.0415 +Bin 13 (0.8500,0.9000) (1.6000,2.0000) 0.0114 +Bin 14 (0.8500,0.9000) (2.0000,15.0000) 0.0000 +Bin 15 (0.9000,1.0000) (0.0000,0.4000) 0.1145 +Bin 16 (0.9000,1.0000) (0.4000,1.2000) 0.1413 +Bin 17 (0.9000,1.0000) (1.2000,1.6000) 0.1202 +Bin 18 (0.9000,1.0000) (1.6000,2.0000) 0.1007 +Bin 19 (0.9000,1.0000) (2.0000,15.0000) 0.0140 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.0005666501 + ( 0-1 ) 0.0005686983 + ( 0-2 ) 0.0005103095 + ( 0-3 ) 0.0005338387 + ( 0-4 ) 0.0000902192 + ( 0-5 ) 0.0003452594 + ( 0-6 ) 0.0004299629 + ( 0-7 ) 0.0004411195 + ( 0-8 ) 0.0000077690 + ( 0-9 ) 0.0000744048 + ( 0-10 ) 0.0001425246 + ( 0-11 ) 0.0003744206 + ( 0-12 ) 0.0000019880 + ( 0-13 ) 0.0000186993 + ( 0-14 ) 0.0000526413 + ( 0-15 ) 0.0003187756 + ( 0-16 ) 0.0000000100 + ( 0-17 ) 0.0000000422 + ( 0-18 ) 0.0000011576 + ( 0-19 ) 0.0000462853 + ( 1-1 ) 0.0009526706 + ( 1-2 ) 0.0006600224 + ( 1-3 ) 0.0006030986 + ( 1-4 ) 0.0001031060 + ( 1-5 ) 0.0004118135 + ( 1-6 ) 0.0005013412 + ( 1-7 ) 0.0005253508 + ( 1-8 ) 0.0000092039 + ( 1-9 ) 0.0000936286 + ( 1-10 ) 0.0001691692 + ( 1-11 ) 0.0004353796 + ( 1-12 ) 0.0000025182 + ( 1-13 ) 0.0000253229 + ( 1-14 ) 0.0000670219 + ( 1-15 ) 0.0003677954 + ( 1-16 ) 0.0000000123 + ( 1-17 ) 0.0000001411 + ( 1-18 ) 0.0000014346 + ( 1-19 ) 0.0000545252 + ( 2-2 ) 0.0008167382 + ( 2-3 ) 0.0006346576 + ( 2-4 ) 0.0000928845 + ( 2-5 ) 0.0003565603 + ( 2-6 ) 0.0004752307 + ( 2-7 ) 0.0004861377 + ( 2-8 ) 0.0000083624 + ( 2-9 ) 0.0000929873 + ( 2-10 ) 0.0001626229 + ( 2-11 ) 0.0004038681 + ( 2-12 ) 0.0000020398 + ( 2-13 ) 0.0000239370 + ( 2-14 ) 0.0000646321 + ( 2-15 ) 0.0003482652 + ( 2-16 ) 0.0000000094 + ( 2-17 ) 0.0000001508 + ( 2-18 ) 0.0000014976 + ( 2-19 ) 0.0000494646 + ( 3-3 ) 0.0009090564 + ( 3-4 ) 0.0001001871 + ( 3-5 ) 0.0003974689 + ( 3-6 ) 0.0005277231 + ( 3-7 ) 0.0005636140 + ( 3-8 ) 0.0000082616 + ( 3-9 ) 0.0000900370 + ( 3-10 ) 0.0001835398 + ( 3-11 ) 0.0004683568 + ( 3-12 ) 0.0000021445 + ( 3-13 ) 0.0000222478 + ( 3-14 ) 0.0000713213 + ( 3-15 ) 0.0004095885 + ( 3-16 ) 0.0000000101 + ( 3-17 ) -0.0000000100 + ( 3-18 ) 0.0000014909 + ( 3-19 ) 0.0000561390 + ( 4-4 ) 0.0000262635 + ( 4-5 ) 0.0000879958 + ( 4-6 ) 0.0001045915 + ( 4-7 ) 0.0001058880 + ( 4-8 ) 0.0000025257 + ( 4-9 ) 0.0000241173 + ( 4-10 ) 0.0000381947 + ( 4-11 ) 0.0000887736 + ( 4-12 ) 0.0000006183 + ( 4-13 ) 0.0000079510 + ( 4-14 ) 0.0000197338 + ( 4-15 ) 0.0000777699 + ( 4-16 ) 0.0000000036 + ( 4-17 ) 0.0000000706 + ( 4-18 ) 0.0000005831 + ( 4-19 ) 0.0000120574 + ( 5-5 ) 0.0005124229 + ( 5-6 ) 0.0004582615 + ( 5-7 ) 0.0004396641 + ( 5-8 ) 0.0000085306 + ( 5-9 ) 0.0001110136 + ( 5-10 ) 0.0001557185 + ( 5-11 ) 0.0003669045 + ( 5-12 ) 0.0000021566 + ( 5-13 ) 0.0000340846 + ( 5-14 ) 0.0000695727 + ( 5-15 ) 0.0003203564 + ( 5-16 ) 0.0000000138 + ( 5-17 ) 0.0000002800 + ( 5-18 ) 0.0000017090 + ( 5-19 ) 0.0000484742 + ( 6-6 ) 0.0007498285 + ( 6-7 ) 0.0005724277 + ( 6-8 ) 0.0000094604 + ( 6-9 ) 0.0001108240 + ( 6-10 ) 0.0002146363 + ( 6-11 ) 0.0004685564 + ( 6-12 ) 0.0000023533 + ( 6-13 ) 0.0000337459 + ( 6-14 ) 0.0000796033 + ( 6-15 ) 0.0004034379 + ( 6-16 ) 0.0000000131 + ( 6-17 ) 0.0000002866 + ( 6-18 ) 0.0000017176 + ( 6-19 ) 0.0000599575 + ( 7-7 ) 0.0006808361 + ( 7-8 ) 0.0000098310 + ( 7-9 ) 0.0001067350 + ( 7-10 ) 0.0001995776 + ( 7-11 ) 0.0004939694 + ( 7-12 ) 0.0000025445 + ( 7-13 ) 0.0000315475 + ( 7-14 ) 0.0000844284 + ( 7-15 ) 0.0004166798 + ( 7-16 ) 0.0000000145 + ( 7-17 ) 0.0000002355 + ( 7-18 ) 0.0000019627 + ( 7-19 ) 0.0000610722 + ( 8-8 ) 0.0000008523 + ( 8-9 ) 0.0000035144 + ( 8-10 ) 0.0000041586 + ( 8-11 ) 0.0000079261 + ( 8-12 ) 0.0000001840 + ( 8-13 ) 0.0000013560 + ( 8-14 ) 0.0000030671 + ( 8-15 ) 0.0000074010 + ( 8-16 ) 0.0000000013 + ( 8-17 ) 0.0000000208 + ( 8-18 ) 0.0000001472 + ( 8-19 ) 0.0000013361 + ( 9-9 ) 0.0001270294 + ( 9-10 ) 0.0000588923 + ( 9-11 ) 0.0000887266 + ( 9-12 ) 0.0000010689 + ( 9-13 ) 0.0000365389 + ( 9-14 ) 0.0000476800 + ( 9-15 ) 0.0000836867 + ( 9-16 ) 0.0000000066 + ( 9-17 ) 0.0000005111 + ( 9-18 ) 0.0000020389 + ( 9-19 ) 0.0000162138 + ( 10-10 ) 0.0001885256 + ( 10-11 ) 0.0001734856 + ( 10-12 ) 0.0000011919 + ( 10-13 ) 0.0000192385 + ( 10-14 ) 0.0000691699 + ( 10-15 ) 0.0001563459 + ( 10-16 ) 0.0000000067 + ( 10-17 ) 0.0000001822 + ( 10-18 ) 0.0000017787 + ( 10-19 ) 0.0000236055 + ( 11-11 ) 0.0005476729 + ( 11-12 ) 0.0000021589 + ( 11-13 ) 0.0000249483 + ( 11-14 ) 0.0000732101 + ( 11-15 ) 0.0003973613 + ( 11-16 ) 0.0000000119 + ( 11-17 ) 0.0000001797 + ( 11-18 ) 0.0000019774 + ( 11-19 ) 0.0000519585 + ( 12-12 ) 0.0000000810 + ( 12-13 ) 0.0000005620 + ( 12-14 ) 0.0000010464 + ( 12-15 ) 0.0000019127 + ( 12-16 ) 0.0000000005 + ( 12-17 ) 0.0000000104 + ( 12-18 ) 0.0000000585 + ( 12-19 ) 0.0000003503 + ( 13-13 ) 0.0000290034 + ( 13-14 ) 0.0000215253 + ( 13-15 ) 0.0000267870 + ( 13-16 ) 0.0000000032 + ( 13-17 ) 0.0000004370 + ( 13-18 ) 0.0000009273 + ( 13-19 ) 0.0000059893 + ( 14-14 ) 0.0001028494 + ( 14-15 ) 0.0000700266 + ( 14-16 ) 0.0000000068 + ( 14-17 ) 0.0000003461 + ( 14-18 ) 0.0000034968 + ( 14-19 ) 0.0000141327 + ( 15-15 ) 0.0004168317 + ( 15-16 ) 0.0000000112 + ( 15-17 ) 0.0000001791 + ( 15-18 ) 0.0000019472 + ( 15-19 ) 0.0000507654 + ( 16-16 ) 0.0000000000 + ( 16-17 ) 0.0000000000 + ( 16-18 ) 0.0000000000 + ( 16-19 ) 0.0000000023 + ( 17-17 ) 0.0000000221 + ( 17-18 ) 0.0000000103 + ( 17-19 ) 0.0000000536 + ( 18-18 ) 0.0000001294 + ( 18-19 ) 0.0000005010 + ( 19-19 ) 0.0000095819 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/Q2.root b/data/T2K/CC1pip/CH/Q2.root new file mode 100644 index 0000000..058de9e Binary files /dev/null and b/data/T2K/CC1pip/CH/Q2.root differ diff --git a/data/T2K/CC1pip/CH/Q2.rootout.root b/data/T2K/CC1pip/CH/Q2.rootout.root new file mode 100644 index 0000000..f23e9b1 Binary files /dev/null and b/data/T2K/CC1pip/CH/Q2.rootout.root differ diff --git a/data/T2K/CC1pip/CH/Q2.txt b/data/T2K/CC1pip/CH/Q2.txt new file mode 100644 index 0000000..e13c180 --- /dev/null +++ b/data/T2K/CC1pip/CH/Q2.txt @@ -0,0 +1,166 @@ + + Variable : Q2 + +(Low end - Up end ) (GeV2) cross-section (10^{-38}cm2/GeV2/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (0.0000,0.1000) 0.0654 +Bin 1 (0.1000,0.2000) 0.0693 +Bin 2 (0.2000,0.3000) 0.0643 +Bin 3 (0.3000,0.4000) 0.0511 +Bin 4 (0.4000,0.5000) 0.0380 +Bin 5 (0.5000,0.6000) 0.0331 +Bin 6 (0.6000,0.7000) 0.0257 +Bin 7 (0.7000,0.8000) 0.0169 +Bin 8 (0.8000,0.9000) 0.0118 +Bin 9 (0.9000,1.0000) 0.0091 +Bin 10 (1.0000,1.2000) 0.0067 +Bin 11 (1.2000,1.4000) 0.0053 +Bin 12 (1.4000,1.6000) 0.0038 +Bin 13 (1.6000,1.8000) 0.0030 +Bin 14 (1.8000,2.0000) 0.0022 +Bin 15 (2.0000,3.3000) 0.0008 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.00021423 + ( 0-1 ) 0.00018968 + ( 0-2 ) 0.00016073 + ( 0-3 ) 0.00012758 + ( 0-4 ) 0.00009625 + ( 0-5 ) 0.00008277 + ( 0-6 ) 0.00006519 + ( 0-7 ) 0.00004784 + ( 0-8 ) 0.00003549 + ( 0-9 ) 0.00002957 + ( 0-10 ) 0.00002330 + ( 0-11 ) 0.00001893 + ( 0-12 ) 0.00001445 + ( 0-13 ) 0.00001238 + ( 0-14 ) 0.00000976 + ( 0-15 ) 0.00000465 + ( 1-1 ) 0.00024190 + ( 1-2 ) 0.00019118 + ( 1-3 ) 0.00014710 + ( 1-4 ) 0.00010846 + ( 1-5 ) 0.00009261 + ( 1-6 ) 0.00007173 + ( 1-7 ) 0.00005271 + ( 1-8 ) 0.00003889 + ( 1-9 ) 0.00003244 + ( 1-10 ) 0.00002575 + ( 1-11 ) 0.00002133 + ( 1-12 ) 0.00001606 + ( 1-13 ) 0.00001387 + ( 1-14 ) 0.00001080 + ( 1-15 ) 0.00000525 + ( 2-2 ) 0.00018773 + ( 2-3 ) 0.00014367 + ( 2-4 ) 0.00010201 + ( 2-5 ) 0.00008628 + ( 2-6 ) 0.00006694 + ( 2-7 ) 0.00004931 + ( 2-8 ) 0.00003627 + ( 2-9 ) 0.00003046 + ( 2-10 ) 0.00002429 + ( 2-11 ) 0.00002031 + ( 2-12 ) 0.00001529 + ( 2-13 ) 0.00001340 + ( 2-14 ) 0.00001049 + ( 2-15 ) 0.00000523 + ( 3-3 ) 0.00012674 + ( 3-4 ) 0.00009028 + ( 3-5 ) 0.00007519 + ( 3-6 ) 0.00005759 + ( 3-7 ) 0.00004205 + ( 3-8 ) 0.00003097 + ( 3-9 ) 0.00002582 + ( 3-10 ) 0.00002074 + ( 3-11 ) 0.00001732 + ( 3-12 ) 0.00001310 + ( 3-13 ) 0.00001156 + ( 3-14 ) 0.00000919 + ( 3-15 ) 0.00000468 + ( 4-4 ) 0.00007242 + ( 4-5 ) 0.00006044 + ( 4-6 ) 0.00004592 + ( 4-7 ) 0.00003299 + ( 4-8 ) 0.00002420 + ( 4-9 ) 0.00002000 + ( 4-10 ) 0.00001597 + ( 4-11 ) 0.00001345 + ( 4-12 ) 0.00001020 + ( 4-13 ) 0.00000911 + ( 4-14 ) 0.00000727 + ( 4-15 ) 0.00000374 + ( 5-5 ) 0.00005521 + ( 5-6 ) 0.00004267 + ( 5-7 ) 0.00003048 + ( 5-8 ) 0.00002212 + ( 5-9 ) 0.00001801 + ( 5-10 ) 0.00001432 + ( 5-11 ) 0.00001208 + ( 5-12 ) 0.00000919 + ( 5-13 ) 0.00000823 + ( 5-14 ) 0.00000656 + ( 5-15 ) 0.00000339 + ( 6-6 ) 0.00003580 + ( 6-7 ) 0.00002554 + ( 6-8 ) 0.00001847 + ( 6-9 ) 0.00001493 + ( 6-10 ) 0.00001187 + ( 6-11 ) 0.00000999 + ( 6-12 ) 0.00000752 + ( 6-13 ) 0.00000677 + ( 6-14 ) 0.00000537 + ( 6-15 ) 0.00000275 + ( 7-7 ) 0.00002050 + ( 7-8 ) 0.00001476 + ( 7-9 ) 0.00001221 + ( 7-10 ) 0.00000962 + ( 7-11 ) 0.00000804 + ( 7-12 ) 0.00000599 + ( 7-13 ) 0.00000545 + ( 7-14 ) 0.00000428 + ( 7-15 ) 0.00000221 + ( 8-8 ) 0.00001193 + ( 8-9 ) 0.00000972 + ( 8-10 ) 0.00000757 + ( 8-11 ) 0.00000613 + ( 8-12 ) 0.00000459 + ( 8-13 ) 0.00000409 + ( 8-14 ) 0.00000323 + ( 8-15 ) 0.00000165 + ( 9-9 ) 0.00000872 + ( 9-10 ) 0.00000671 + ( 9-11 ) 0.00000549 + ( 9-12 ) 0.00000400 + ( 9-13 ) 0.00000356 + ( 9-14 ) 0.00000277 + ( 9-15 ) 0.00000142 + ( 10-10 ) 0.00000570 + ( 10-11 ) 0.00000466 + ( 10-12 ) 0.00000343 + ( 10-13 ) 0.00000307 + ( 10-14 ) 0.00000239 + ( 10-15 ) 0.00000123 + ( 11-11 ) 0.00000422 + ( 11-12 ) 0.00000311 + ( 11-13 ) 0.00000284 + ( 11-14 ) 0.00000223 + ( 11-15 ) 0.00000119 + ( 12-12 ) 0.00000255 + ( 12-13 ) 0.00000231 + ( 12-14 ) 0.00000186 + ( 12-15 ) 0.00000099 + ( 13-13 ) 0.00000237 + ( 13-14 ) 0.00000184 + ( 13-15 ) 0.00000101 + ( 14-14 ) 0.00000160 + ( 14-15 ) 0.00000090 + ( 15-15 ) 0.00000056 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/Thetapimu.root b/data/T2K/CC1pip/CH/Thetapimu.root new file mode 100644 index 0000000..a11fd49 Binary files /dev/null and b/data/T2K/CC1pip/CH/Thetapimu.root differ diff --git a/data/T2K/CC1pip/CH/Thetapimu.rootout.root b/data/T2K/CC1pip/CH/Thetapimu.rootout.root new file mode 100644 index 0000000..154bf5a Binary files /dev/null and b/data/T2K/CC1pip/CH/Thetapimu.rootout.root differ diff --git a/data/T2K/CC1pip/CH/Thetapimu.txt b/data/T2K/CC1pip/CH/Thetapimu.txt new file mode 100644 index 0000000..96fbcbe --- /dev/null +++ b/data/T2K/CC1pip/CH/Thetapimu.txt @@ -0,0 +1,166 @@ + + Variable : Theta (pi,mu) (rads) + +(Low end - Up end ) (rads) cross-section (10^{-38}cm2/rad/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (0.0000,0.1000) 0.0039 +Bin 1 (0.1000,0.2000) 0.0096 +Bin 2 (0.2000,0.3000) 0.0160 +Bin 3 (0.3000,0.4000) 0.0223 +Bin 4 (0.4000,0.5000) 0.0259 +Bin 5 (0.5000,0.6000) 0.0295 +Bin 6 (0.6000,0.7000) 0.0313 +Bin 7 (0.7000,0.8000) 0.0351 +Bin 8 (0.8000,0.9000) 0.0321 +Bin 9 (0.9000,1.0000) 0.0298 +Bin 10 (1.0000,1.2000) 0.0299 +Bin 11 (1.2000,1.4000) 0.0318 +Bin 12 (1.4000,1.6000) 0.0194 +Bin 13 (1.6000,1.8000) 0.0104 +Bin 14 (1.8000,2.0000) 0.0041 +Bin 15 (2.0000,3.1415) 0.0014 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 1-1 ) 0.00000300 + ( 1-2 ) 0.00000318 + ( 1-3 ) 0.00000354 + ( 1-4 ) 0.00000445 + ( 1-5 ) 0.00000525 + ( 1-6 ) 0.00000586 + ( 1-7 ) 0.00000624 + ( 1-8 ) 0.00000664 + ( 1-9 ) 0.00000638 + ( 1-10 ) 0.00000605 + ( 1-11 ) 0.00000618 + ( 1-12 ) 0.00000585 + ( 1-13 ) 0.00000372 + ( 1-14 ) 0.00000200 + ( 1-15 ) 0.00000076 + ( 1-16 ) 0.00000015 + ( 2-2 ) 0.00001147 + ( 2-3 ) 0.00001106 + ( 2-4 ) 0.00001296 + ( 2-5 ) 0.00001519 + ( 2-6 ) 0.00001628 + ( 2-7 ) 0.00001794 + ( 2-8 ) 0.00001923 + ( 2-9 ) 0.00001788 + ( 2-10 ) 0.00001756 + ( 2-11 ) 0.00001734 + ( 2-12 ) 0.00001527 + ( 2-13 ) 0.00000981 + ( 2-14 ) 0.00000496 + ( 2-15 ) 0.00000210 + ( 2-16 ) 0.00000059 + ( 3-3 ) 0.00002394 + ( 3-4 ) 0.00002256 + ( 3-5 ) 0.00002509 + ( 3-6 ) 0.00002699 + ( 3-7 ) 0.00003013 + ( 3-8 ) 0.00003102 + ( 3-9 ) 0.00002961 + ( 3-10 ) 0.00002758 + ( 3-11 ) 0.00002691 + ( 3-12 ) 0.00002406 + ( 3-13 ) 0.00001500 + ( 3-14 ) 0.00000804 + ( 3-15 ) 0.00000315 + ( 3-16 ) 0.00000080 + ( 4-4 ) 0.00003734 + ( 4-5 ) 0.00003513 + ( 4-6 ) 0.00003568 + ( 4-7 ) 0.00003954 + ( 4-8 ) 0.00004162 + ( 4-9 ) 0.00003923 + ( 4-10 ) 0.00003625 + ( 4-11 ) 0.00003509 + ( 4-12 ) 0.00003260 + ( 4-13 ) 0.00002040 + ( 4-14 ) 0.00001046 + ( 4-15 ) 0.00000414 + ( 4-16 ) 0.00000115 + ( 5-5 ) 0.00005177 + ( 5-6 ) 0.00004526 + ( 5-7 ) 0.00004765 + ( 5-8 ) 0.00005024 + ( 5-9 ) 0.00004699 + ( 5-10 ) 0.00004413 + ( 5-11 ) 0.00004270 + ( 5-12 ) 0.00003915 + ( 5-13 ) 0.00002479 + ( 5-14 ) 0.00001289 + ( 5-15 ) 0.00000496 + ( 5-16 ) 0.00000136 + ( 6-6 ) 0.00005906 + ( 6-7 ) 0.00005358 + ( 6-8 ) 0.00005466 + ( 6-9 ) 0.00005176 + ( 6-10 ) 0.00004726 + ( 6-11 ) 0.00004622 + ( 6-12 ) 0.00004236 + ( 6-13 ) 0.00002692 + ( 6-14 ) 0.00001421 + ( 6-15 ) 0.00000577 + ( 6-16 ) 0.00000141 + ( 7-7 ) 0.00007188 + ( 7-8 ) 0.00006461 + ( 7-9 ) 0.00005849 + ( 7-10 ) 0.00005310 + ( 7-11 ) 0.00005241 + ( 7-12 ) 0.00004771 + ( 7-13 ) 0.00003023 + ( 7-14 ) 0.00001546 + ( 7-15 ) 0.00000652 + ( 7-16 ) 0.00000157 + ( 8-8 ) 0.00008295 + ( 8-9 ) 0.00006546 + ( 8-10 ) 0.00005810 + ( 8-11 ) 0.00005815 + ( 8-12 ) 0.00005356 + ( 8-13 ) 0.00003414 + ( 8-14 ) 0.00001733 + ( 8-15 ) 0.00000717 + ( 8-16 ) 0.00000174 + ( 9-9 ) 0.00007636 + ( 9-10 ) 0.00005787 + ( 9-11 ) 0.00005499 + ( 9-12 ) 0.00005147 + ( 9-13 ) 0.00003291 + ( 9-14 ) 0.00001688 + ( 9-15 ) 0.00000733 + ( 9-16 ) 0.00000166 + ( 10-10 ) 0.00006983 + ( 10-11 ) 0.00005616 + ( 10-12 ) 0.00004833 + ( 10-13 ) 0.00003118 + ( 10-14 ) 0.00001639 + ( 10-15 ) 0.00000714 + ( 10-16 ) 0.00000165 + ( 11-11 ) 0.00006693 + ( 11-12 ) 0.00005218 + ( 11-13 ) 0.00003290 + ( 11-14 ) 0.00001702 + ( 11-15 ) 0.00000726 + ( 11-16 ) 0.00000157 + ( 12-12 ) 0.00006748 + ( 12-13 ) 0.00003475 + ( 12-14 ) 0.00001690 + ( 12-15 ) 0.00000707 + ( 12-16 ) 0.00000160 + ( 13-13 ) 0.00003325 + ( 13-14 ) 0.00001299 + ( 13-15 ) 0.00000500 + ( 13-16 ) 0.00000103 + ( 14-14 ) 0.00001326 + ( 14-15 ) 0.00000304 + ( 14-16 ) 0.00000054 + ( 15-15 ) 0.00000504 + ( 15-16 ) 0.00000039 + ( 16-16 ) 0.00000064 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/Thetapion.root b/data/T2K/CC1pip/CH/Thetapion.root new file mode 100644 index 0000000..6c9a432 Binary files /dev/null and b/data/T2K/CC1pip/CH/Thetapion.root differ diff --git a/data/T2K/CC1pip/CH/Thetapion.rootout.root b/data/T2K/CC1pip/CH/Thetapion.rootout.root new file mode 100644 index 0000000..53e2b5a Binary files /dev/null and b/data/T2K/CC1pip/CH/Thetapion.rootout.root differ diff --git a/data/T2K/CC1pip/CH/Thetapion.txt b/data/T2K/CC1pip/CH/Thetapion.txt new file mode 100644 index 0000000..2d09a2e --- /dev/null +++ b/data/T2K/CC1pip/CH/Thetapion.txt @@ -0,0 +1,118 @@ + + Variable : Theta_pion + +(Low end - Up end ) (rad) cross-section (10^{-38}cm2/rad/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (0.0000,0.1000) 0.0111 +Bin 1 (0.1000,0.2000) 0.0187 +Bin 2 (0.2000,0.3000) 0.0338 +Bin 3 (0.3000,0.4000) 0.0480 +Bin 4 (0.4000,0.5000) 0.0431 +Bin 5 (0.5000,0.6000) 0.0517 +Bin 6 (0.6000,0.7000) 0.0473 +Bin 7 (0.7000,0.8000) 0.0340 +Bin 8 (0.8000,0.9000) 0.0333 +Bin 9 (0.9000,1.0000) 0.0312 +Bin 10 (1.0000,1.2000) 0.0216 +Bin 11 (1.2000,1.4000) 0.0178 +Bin 12 (1.4000,3.1415) 0.0013 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.00002172 + ( 0-1 ) 0.00002430 + ( 0-2 ) 0.00002696 + ( 0-3 ) 0.00003337 + ( 0-4 ) 0.00003079 + ( 0-5 ) 0.00003312 + ( 0-6 ) 0.00002957 + ( 0-7 ) 0.00002266 + ( 0-8 ) 0.00002138 + ( 0-9 ) 0.00001937 + ( 0-10 ) 0.00001515 + ( 0-11 ) 0.00001401 + ( 0-12 ) 0.00000117 + ( 1-1 ) 0.00005862 + ( 1-2 ) 0.00005504 + ( 1-3 ) 0.00006442 + ( 1-4 ) 0.00005963 + ( 1-5 ) 0.00006369 + ( 1-6 ) 0.00005710 + ( 1-7 ) 0.00004586 + ( 1-8 ) 0.00004401 + ( 1-9 ) 0.00004036 + ( 1-10 ) 0.00003014 + ( 1-11 ) 0.00002675 + ( 1-12 ) 0.00000220 + ( 2-2 ) 0.00008754 + ( 2-3 ) 0.00008892 + ( 2-4 ) 0.00008083 + ( 2-5 ) 0.00008640 + ( 2-6 ) 0.00007711 + ( 2-7 ) 0.00006260 + ( 2-8 ) 0.00006143 + ( 2-9 ) 0.00005505 + ( 2-10 ) 0.00004030 + ( 2-11 ) 0.00003421 + ( 2-12 ) 0.00000277 + ( 3-3 ) 0.00013540 + ( 3-4 ) 0.00010860 + ( 3-5 ) 0.00011292 + ( 3-6 ) 0.00010387 + ( 3-7 ) 0.00008113 + ( 3-8 ) 0.00007977 + ( 3-9 ) 0.00007120 + ( 3-10 ) 0.00005329 + ( 3-11 ) 0.00004479 + ( 3-12 ) 0.00000370 + ( 4-4 ) 0.00011936 + ( 4-5 ) 0.00011004 + ( 4-6 ) 0.00009835 + ( 4-7 ) 0.00007740 + ( 4-8 ) 0.00007642 + ( 4-9 ) 0.00006871 + ( 4-10 ) 0.00005100 + ( 4-11 ) 0.00004172 + ( 4-12 ) 0.00000337 + ( 5-5 ) 0.00014181 + ( 5-6 ) 0.00011018 + ( 5-7 ) 0.00008384 + ( 5-8 ) 0.00008141 + ( 5-9 ) 0.00007384 + ( 5-10 ) 0.00005519 + ( 5-11 ) 0.00004808 + ( 5-12 ) 0.00000369 + ( 6-6 ) 0.00012355 + ( 6-7 ) 0.00008184 + ( 6-8 ) 0.00007513 + ( 6-9 ) 0.00006766 + ( 6-10 ) 0.00005024 + ( 6-11 ) 0.00004274 + ( 6-12 ) 0.00000342 + ( 7-7 ) 0.00008216 + ( 7-8 ) 0.00006588 + ( 7-9 ) 0.00005505 + ( 7-10 ) 0.00004114 + ( 7-11 ) 0.00003449 + ( 7-12 ) 0.00000274 + ( 8-8 ) 0.00008376 + ( 8-9 ) 0.00005940 + ( 8-10 ) 0.00004156 + ( 8-11 ) 0.00003395 + ( 8-12 ) 0.00000261 + ( 9-9 ) 0.00007216 + ( 9-10 ) 0.00004037 + ( 9-11 ) 0.00003235 + ( 9-12 ) 0.00000254 + ( 10-10 ) 0.00004131 + ( 10-11 ) 0.00002616 + ( 10-12 ) 0.00000199 + ( 11-11 ) 0.00003700 + ( 11-12 ) 0.00000201 + ( 12-12 ) 0.00000016 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/phi_adler.root b/data/T2K/CC1pip/CH/phi_adler.root new file mode 100644 index 0000000..7648d41 Binary files /dev/null and b/data/T2K/CC1pip/CH/phi_adler.root differ diff --git a/data/T2K/CC1pip/CH/phi_adler.rootout.root b/data/T2K/CC1pip/CH/phi_adler.rootout.root new file mode 100644 index 0000000..2ca4a16 Binary files /dev/null and b/data/T2K/CC1pip/CH/phi_adler.rootout.root differ diff --git a/data/T2K/CC1pip/CH/phi_adler.txt b/data/T2K/CC1pip/CH/phi_adler.txt new file mode 100644 index 0000000..a550f5c --- /dev/null +++ b/data/T2K/CC1pip/CH/phi_adler.txt @@ -0,0 +1,166 @@ + + Variable : Phi_Adler + +(Low end - Up end ) (rads) cross-section (10^{-38}cm2/rad/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (-3.1415,-2.8000) 0.0070 +Bin 1 (-2.8000,-2.4000) 0.0060 +Bin 2 (-2.4000,-2.0000) 0.0071 +Bin 3 (-2.0000,-1.6000) 0.0063 +Bin 4 (-1.6000,-1.2000) 0.0071 +Bin 5 (-1.2000,-0.8000) 0.0079 +Bin 6 (-0.8000,-0.4000) 0.0077 +Bin 7 (-0.4000,0.0000) 0.0084 +Bin 8 (0.0000,0.4000) 0.0085 +Bin 9 (0.4000,0.8000) 0.0094 +Bin 10 (0.8000,1.2000) 0.0070 +Bin 11 (1.2000,1.6000) 0.0055 +Bin 12 (1.6000,2.0000) 0.0060 +Bin 13 (2.0000,2.4000) 0.0040 +Bin 14 (2.4000,2.8000) 0.0046 +Bin 15 (2.8000,3.1415) 0.0081 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.00000420 + ( 0-1 ) 0.00000276 + ( 0-2 ) 0.00000270 + ( 0-3 ) 0.00000236 + ( 0-4 ) 0.00000255 + ( 0-5 ) 0.00000259 + ( 0-6 ) 0.00000266 + ( 0-7 ) 0.00000277 + ( 0-8 ) 0.00000286 + ( 0-9 ) 0.00000305 + ( 0-10 ) 0.00000232 + ( 0-11 ) 0.00000202 + ( 0-12 ) 0.00000228 + ( 0-13 ) 0.00000169 + ( 0-14 ) 0.00000217 + ( 0-15 ) 0.00000396 + ( 1-1 ) 0.00000295 + ( 1-2 ) 0.00000262 + ( 1-3 ) 0.00000223 + ( 1-4 ) 0.00000225 + ( 1-5 ) 0.00000236 + ( 1-6 ) 0.00000233 + ( 1-7 ) 0.00000241 + ( 1-8 ) 0.00000246 + ( 1-9 ) 0.00000264 + ( 1-10 ) 0.00000202 + ( 1-11 ) 0.00000172 + ( 1-12 ) 0.00000195 + ( 1-13 ) 0.00000142 + ( 1-14 ) 0.00000179 + ( 1-15 ) 0.00000298 + ( 2-2 ) 0.00000357 + ( 2-3 ) 0.00000288 + ( 2-4 ) 0.00000273 + ( 2-5 ) 0.00000283 + ( 2-6 ) 0.00000283 + ( 2-7 ) 0.00000292 + ( 2-8 ) 0.00000295 + ( 2-9 ) 0.00000315 + ( 2-10 ) 0.00000245 + ( 2-11 ) 0.00000216 + ( 2-12 ) 0.00000234 + ( 2-13 ) 0.00000168 + ( 2-14 ) 0.00000201 + ( 2-15 ) 0.00000320 + ( 3-3 ) 0.00000339 + ( 3-4 ) 0.00000295 + ( 3-5 ) 0.00000288 + ( 3-6 ) 0.00000298 + ( 3-7 ) 0.00000310 + ( 3-8 ) 0.00000314 + ( 3-9 ) 0.00000322 + ( 3-10 ) 0.00000251 + ( 3-11 ) 0.00000217 + ( 3-12 ) 0.00000230 + ( 3-13 ) 0.00000171 + ( 3-14 ) 0.00000193 + ( 3-15 ) 0.00000294 + ( 4-4 ) 0.00000384 + ( 4-5 ) 0.00000333 + ( 4-6 ) 0.00000327 + ( 4-7 ) 0.00000344 + ( 4-8 ) 0.00000346 + ( 4-9 ) 0.00000355 + ( 4-10 ) 0.00000276 + ( 4-11 ) 0.00000233 + ( 4-12 ) 0.00000242 + ( 4-13 ) 0.00000180 + ( 4-14 ) 0.00000202 + ( 4-15 ) 0.00000303 + ( 5-5 ) 0.00000419 + ( 5-6 ) 0.00000365 + ( 5-7 ) 0.00000370 + ( 5-8 ) 0.00000374 + ( 5-9 ) 0.00000379 + ( 5-10 ) 0.00000294 + ( 5-11 ) 0.00000244 + ( 5-12 ) 0.00000252 + ( 5-13 ) 0.00000187 + ( 5-14 ) 0.00000208 + ( 5-15 ) 0.00000312 + ( 6-6 ) 0.00000469 + ( 6-7 ) 0.00000419 + ( 6-8 ) 0.00000403 + ( 6-9 ) 0.00000408 + ( 6-10 ) 0.00000315 + ( 6-11 ) 0.00000264 + ( 6-12 ) 0.00000269 + ( 6-13 ) 0.00000202 + ( 6-14 ) 0.00000211 + ( 6-15 ) 0.00000307 + ( 7-7 ) 0.00000559 + ( 7-8 ) 0.00000463 + ( 7-9 ) 0.00000447 + ( 7-10 ) 0.00000337 + ( 7-11 ) 0.00000284 + ( 7-12 ) 0.00000288 + ( 7-13 ) 0.00000215 + ( 7-14 ) 0.00000223 + ( 7-15 ) 0.00000324 + ( 8-8 ) 0.00000569 + ( 8-9 ) 0.00000476 + ( 8-10 ) 0.00000348 + ( 8-11 ) 0.00000283 + ( 8-12 ) 0.00000288 + ( 8-13 ) 0.00000217 + ( 8-14 ) 0.00000228 + ( 8-15 ) 0.00000330 + ( 9-9 ) 0.00000600 + ( 9-10 ) 0.00000376 + ( 9-11 ) 0.00000291 + ( 9-12 ) 0.00000301 + ( 9-13 ) 0.00000224 + ( 9-14 ) 0.00000238 + ( 9-15 ) 0.00000347 + ( 10-10 ) 0.00000371 + ( 10-11 ) 0.00000250 + ( 10-12 ) 0.00000241 + ( 10-13 ) 0.00000178 + ( 10-14 ) 0.00000185 + ( 10-15 ) 0.00000268 + ( 11-11 ) 0.00000266 + ( 11-12 ) 0.00000230 + ( 11-13 ) 0.00000155 + ( 11-14 ) 0.00000164 + ( 11-15 ) 0.00000237 + ( 12-12 ) 0.00000297 + ( 12-13 ) 0.00000183 + ( 12-14 ) 0.00000182 + ( 12-15 ) 0.00000271 + ( 13-13 ) 0.00000166 + ( 13-14 ) 0.00000152 + ( 13-15 ) 0.00000200 + ( 14-14 ) 0.00000221 + ( 14-15 ) 0.00000283 + ( 15-15 ) 0.00000574 +---------------------------------------------------------------------------------------- + diff --git a/data/T2K/CC1pip/CH/readme b/data/T2K/CC1pip/CH/readme new file mode 100644 index 0000000..ae22466 --- /dev/null +++ b/data/T2K/CC1pip/CH/readme @@ -0,0 +1,3 @@ +Initial data release of Raquel Castillo's CC1pi+ CH measurement on CH + +Released for Pittsburgh TENSIONS workshop 2019 diff --git a/data/T2K/CC1pip/CH/theta_adler.root b/data/T2K/CC1pip/CH/theta_adler.root new file mode 100644 index 0000000..cf26e67 Binary files /dev/null and b/data/T2K/CC1pip/CH/theta_adler.root differ diff --git a/data/T2K/CC1pip/CH/theta_adler.rootout.root b/data/T2K/CC1pip/CH/theta_adler.rootout.root new file mode 100644 index 0000000..544ad36 Binary files /dev/null and b/data/T2K/CC1pip/CH/theta_adler.rootout.root differ diff --git a/data/T2K/CC1pip/CH/theta_adler.txt b/data/T2K/CC1pip/CH/theta_adler.txt new file mode 100644 index 0000000..f967d48 --- /dev/null +++ b/data/T2K/CC1pip/CH/theta_adler.txt @@ -0,0 +1,244 @@ + + Variable : Theta_Adler + +(Low end - Up end ) (rads) cross-section (10^{-38}cm2/rad/Nucleon) +--------------------------------------------------------------------------------------- +Bin 0 (-1.0000,-0.9000) 0.0003 +Bin 1 (-0.9000,-0.8000) -0.0008 +Bin 2 (-0.8000,-0.7000) -0.0004 +Bin 3 (-0.7000,-0.6000) 0.0007 +Bin 4 (-0.6000,-0.5000) -0.0003 +Bin 5 (-0.5000,-0.4000) 0.0048 +Bin 6 (-0.4000,-0.3000) 0.0133 +Bin 7 (-0.3000,-0.2000) 0.0152 +Bin 8 (-0.2000,-0.1000) 0.0155 +Bin 9 (-0.1000,0.0000) 0.0177 +Bin 10 (0.0000,0.1000) 0.0248 +Bin 11 (0.1000,0.2000) 0.0260 +Bin 12 (0.2000,0.3000) 0.0290 +Bin 13 (0.3000,0.4000) 0.0284 +Bin 14 (0.4000,0.5000) 0.0355 +Bin 15 (0.5000,0.6000) 0.0379 +Bin 16 (0.6000,0.7000) 0.0374 +Bin 17 (0.7000,0.8000) 0.0358 +Bin 18 (0.8000,0.9000) 0.0331 +Bin 19 (0.9000,1.0000) 0.0607 +---------------------------------------------------------------------------------------- + + Covariance matrix + +(bin1-bin2) Cov(bin1,bin2) [= Cov(bin2,bin1)] +---------------------------------------------------------------------------------------- + + ( 0-0 ) 0.00000003 + ( 0-1 ) 0.00000002 + ( 0-2 ) 0.00000004 + ( 0-3 ) 0.00000018 + ( 0-4 ) 0.00000011 + ( 0-5 ) 0.00000020 + ( 0-6 ) 0.00000036 + ( 0-7 ) 0.00000036 + ( 0-8 ) 0.00000035 + ( 0-9 ) 0.00000041 + ( 0-10 ) 0.00000045 + ( 0-11 ) 0.00000042 + ( 0-12 ) 0.00000047 + ( 0-13 ) 0.00000044 + ( 0-14 ) 0.00000046 + ( 0-15 ) 0.00000047 + ( 0-16 ) 0.00000051 + ( 0-17 ) 0.00000055 + ( 0-18 ) 0.00000063 + ( 0-19 ) 0.00000085 + ( 1-1 ) 0.00000103 + ( 1-2 ) 0.00000084 + ( 1-3 ) 0.00000266 + ( 1-4 ) 0.00000085 + ( 1-5 ) 0.00000110 + ( 1-6 ) 0.00000185 + ( 1-7 ) 0.00000200 + ( 1-8 ) 0.00000204 + ( 1-9 ) 0.00000187 + ( 1-10 ) 0.00000222 + ( 1-11 ) 0.00000160 + ( 1-12 ) 0.00000150 + ( 1-13 ) 0.00000133 + ( 1-14 ) 0.00000131 + ( 1-15 ) 0.00000153 + ( 1-16 ) 0.00000138 + ( 1-17 ) 0.00000162 + ( 1-18 ) 0.00000170 + ( 1-19 ) 0.00000231 + ( 2-2 ) 0.00000141 + ( 2-3 ) 0.00000327 + ( 2-4 ) 0.00000127 + ( 2-5 ) 0.00000176 + ( 2-6 ) 0.00000277 + ( 2-7 ) 0.00000283 + ( 2-8 ) 0.00000277 + ( 2-9 ) 0.00000272 + ( 2-10 ) 0.00000328 + ( 2-11 ) 0.00000271 + ( 2-12 ) 0.00000266 + ( 2-13 ) 0.00000246 + ( 2-14 ) 0.00000251 + ( 2-15 ) 0.00000288 + ( 2-16 ) 0.00000279 + ( 2-17 ) 0.00000294 + ( 2-18 ) 0.00000323 + ( 2-19 ) 0.00000469 + ( 3-3 ) 0.00002023 + ( 3-4 ) 0.00000536 + ( 3-5 ) 0.00000917 + ( 3-6 ) 0.00001527 + ( 3-7 ) 0.00001551 + ( 3-8 ) 0.00001625 + ( 3-9 ) 0.00001479 + ( 3-10 ) 0.00001851 + ( 3-11 ) 0.00001557 + ( 3-12 ) 0.00001545 + ( 3-13 ) 0.00001464 + ( 3-14 ) 0.00001603 + ( 3-15 ) 0.00001726 + ( 3-16 ) 0.00001654 + ( 3-17 ) 0.00001788 + ( 3-18 ) 0.00001912 + ( 3-19 ) 0.00002859 + ( 4-4 ) 0.00000452 + ( 4-5 ) 0.00000479 + ( 4-6 ) 0.00000516 + ( 4-7 ) 0.00000507 + ( 4-8 ) 0.00000516 + ( 4-9 ) 0.00000494 + ( 4-10 ) 0.00000565 + ( 4-11 ) 0.00000500 + ( 4-12 ) 0.00000492 + ( 4-13 ) 0.00000480 + ( 4-14 ) 0.00000499 + ( 4-15 ) 0.00000520 + ( 4-16 ) 0.00000530 + ( 4-17 ) 0.00000583 + ( 4-18 ) 0.00000688 + ( 4-19 ) 0.00000948 + ( 5-5 ) 0.00001092 + ( 5-6 ) 0.00001296 + ( 5-7 ) 0.00001284 + ( 5-8 ) 0.00001175 + ( 5-9 ) 0.00001118 + ( 5-10 ) 0.00001387 + ( 5-11 ) 0.00001281 + ( 5-12 ) 0.00001270 + ( 5-13 ) 0.00001236 + ( 5-14 ) 0.00001314 + ( 5-15 ) 0.00001447 + ( 5-16 ) 0.00001444 + ( 5-17 ) 0.00001496 + ( 5-18 ) 0.00001608 + ( 5-19 ) 0.00002378 + ( 6-6 ) 0.00002681 + ( 6-7 ) 0.00002543 + ( 6-8 ) 0.00002282 + ( 6-9 ) 0.00002108 + ( 6-10 ) 0.00002719 + ( 6-11 ) 0.00002378 + ( 6-12 ) 0.00002483 + ( 6-13 ) 0.00002319 + ( 6-14 ) 0.00002618 + ( 6-15 ) 0.00002843 + ( 6-16 ) 0.00002772 + ( 6-17 ) 0.00002863 + ( 6-18 ) 0.00002847 + ( 6-19 ) 0.00004488 + ( 7-7 ) 0.00002874 + ( 7-8 ) 0.00002610 + ( 7-9 ) 0.00002343 + ( 7-10 ) 0.00002931 + ( 7-11 ) 0.00002611 + ( 7-12 ) 0.00002670 + ( 7-13 ) 0.00002522 + ( 7-14 ) 0.00002782 + ( 7-15 ) 0.00003043 + ( 7-16 ) 0.00002976 + ( 7-17 ) 0.00003096 + ( 7-18 ) 0.00003099 + ( 7-19 ) 0.00004912 + ( 8-8 ) 0.00002830 + ( 8-9 ) 0.00002615 + ( 8-10 ) 0.00003024 + ( 8-11 ) 0.00002671 + ( 8-12 ) 0.00002731 + ( 8-13 ) 0.00002592 + ( 8-14 ) 0.00002917 + ( 8-15 ) 0.00003156 + ( 8-16 ) 0.00003121 + ( 8-17 ) 0.00003237 + ( 8-18 ) 0.00003258 + ( 8-19 ) 0.00005198 + ( 9-9 ) 0.00002967 + ( 9-10 ) 0.00003342 + ( 9-11 ) 0.00002841 + ( 9-12 ) 0.00002879 + ( 9-13 ) 0.00002692 + ( 9-14 ) 0.00003028 + ( 9-15 ) 0.00003259 + ( 9-16 ) 0.00003219 + ( 9-17 ) 0.00003345 + ( 9-18 ) 0.00003438 + ( 9-19 ) 0.00005416 + ( 10-10 ) 0.00004464 + ( 10-11 ) 0.00003844 + ( 10-12 ) 0.00003755 + ( 10-13 ) 0.00003459 + ( 10-14 ) 0.00003924 + ( 10-15 ) 0.00004232 + ( 10-16 ) 0.00004126 + ( 10-17 ) 0.00004235 + ( 10-18 ) 0.00004275 + ( 10-19 ) 0.00006888 + ( 11-11 ) 0.00003967 + ( 11-12 ) 0.00003847 + ( 11-13 ) 0.00003431 + ( 11-14 ) 0.00003817 + ( 11-15 ) 0.00004112 + ( 11-16 ) 0.00004053 + ( 11-17 ) 0.00004111 + ( 11-18 ) 0.00004188 + ( 11-19 ) 0.00006713 + ( 12-12 ) 0.00004369 + ( 12-13 ) 0.00003819 + ( 12-14 ) 0.00004157 + ( 12-15 ) 0.00004389 + ( 12-16 ) 0.00004346 + ( 12-17 ) 0.00004392 + ( 12-18 ) 0.00004400 + ( 12-19 ) 0.00007083 + ( 13-13 ) 0.00003996 + ( 13-14 ) 0.00004284 + ( 13-15 ) 0.00004386 + ( 13-16 ) 0.00004288 + ( 13-17 ) 0.00004314 + ( 13-18 ) 0.00004312 + ( 13-19 ) 0.00006950 + ( 14-14 ) 0.00005629 + ( 14-15 ) 0.00005574 + ( 14-16 ) 0.00005208 + ( 14-17 ) 0.00005119 + ( 14-18 ) 0.00005072 + ( 14-19 ) 0.00008344 + ( 15-15 ) 0.00006577 + ( 15-16 ) 0.00006012 + ( 15-17 ) 0.00005626 + ( 15-18 ) 0.00005516 + ( 15-19 ) 0.00009104 + ( 16-16 ) 0.00006609 + ( 16-17 ) 0.00006111 + ( 16-18 ) 0.00005678 + ( 16-19 ) 0.00009200 + ( 17-17 ) 0.00006786 + ( 17-18 ) 0.00006217 + ( 17-19 ) 0.00009420 + ( 18-18 ) 0.00007251 + ( 18-19 ) 0.00010122 + ( 19-19 ) 0.00019005 +---------------------------------------------------------------------------------------- + diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index d137b57..9b8a0b9 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,1414 +1,1423 @@ #include "SampleList.h" #ifndef __NO_ANL__ #include "ANL_CCQE_Evt_1DQ2_nu.h" #include "ANL_CCQE_XSec_1DEnu_nu.h" // ANL CC1ppip #include "ANL_CC1ppip_Evt_1DQ2_nu.h" #include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" #include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" #include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h" #include "ANL_CC1ppip_Evt_1Dphi_nu.h" #include "ANL_CC1ppip_Evt_1Dppi_nu.h" #include "ANL_CC1ppip_Evt_1Dthpr_nu.h" #include "ANL_CC1ppip_XSec_1DEnu_nu.h" #include "ANL_CC1ppip_XSec_1DQ2_nu.h" // ANL CC1npip #include "ANL_CC1npip_Evt_1DQ2_nu.h" #include "ANL_CC1npip_Evt_1DcosmuStar_nu.h" #include "ANL_CC1npip_Evt_1Dppi_nu.h" #include "ANL_CC1npip_XSec_1DEnu_nu.h" // ANL CC1pi0 #include "ANL_CC1pi0_Evt_1DQ2_nu.h" #include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h" #include "ANL_CC1pi0_XSec_1DEnu_nu.h" // ANL NC1npip (mm, exotic!) #include "ANL_NC1npip_Evt_1Dppi_nu.h" // ANL NC1ppim (mm, exotic!) #include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h" #include "ANL_NC1ppim_XSec_1DEnu_nu.h" // ANL CC2pi 1pim1pip (mm, even more exotic!) #include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pip (mm, even more exotic!) #include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h" #include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pi0 (mm, even more exotic!) #include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h" #endif #ifndef __NO_ArgoNeuT__ // ArgoNeuT CC1Pi #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h" // ArgoNeuT CC-inclusive #include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h" #include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h" #include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h" #endif #ifndef __NO_BNL__ // BNL CCQE #include "BNL_CCQE_Evt_1DQ2_nu.h" #include "BNL_CCQE_XSec_1DEnu_nu.h" // BNL CC1ppip #include "BNL_CC1ppip_Evt_1DQ2_nu.h" #include "BNL_CC1ppip_Evt_1DQ2_nu.h" #include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h" #include "BNL_CC1ppip_Evt_1Dphi_nu.h" #include "BNL_CC1ppip_XSec_1DEnu_nu.h" // BNL CC1npip #include "BNL_CC1npip_Evt_1DQ2_nu.h" #include "BNL_CC1npip_XSec_1DEnu_nu.h" // BNL CC1pi0 #include "BNL_CC1pi0_Evt_1DQ2_nu.h" #include "BNL_CC1pi0_XSec_1DEnu_nu.h" #endif #ifndef __NO_FNAL__ // FNAL CCQE #include "FNAL_CCQE_Evt_1DQ2_nu.h" // FNAL CC1ppip #include "FNAL_CC1ppip_Evt_1DQ2_nu.h" #include "FNAL_CC1ppip_XSec_1DEnu_nu.h" #include "FNAL_CC1ppip_XSec_1DQ2_nu.h" // FNAL CC1ppim #include "FNAL_CC1ppim_XSec_1DEnu_antinu.h" #endif #ifndef __NO_BEBC__ // BEBC CCQE #include "BEBC_CCQE_XSec_1DQ2_nu.h" // BEBC CC1ppip #include "BEBC_CC1ppip_XSec_1DEnu_nu.h" #include "BEBC_CC1ppip_XSec_1DQ2_nu.h" // BEBC CC1npip #include "BEBC_CC1npip_XSec_1DEnu_nu.h" #include "BEBC_CC1npip_XSec_1DQ2_nu.h" // BEBC CC1pi0 #include "BEBC_CC1pi0_XSec_1DEnu_nu.h" #include "BEBC_CC1pi0_XSec_1DQ2_nu.h" // BEBC CC1npim #include "BEBC_CC1npim_XSec_1DEnu_antinu.h" #include "BEBC_CC1npim_XSec_1DQ2_antinu.h" // BEBC CC1ppim #include "BEBC_CC1ppim_XSec_1DEnu_antinu.h" #include "BEBC_CC1ppim_XSec_1DQ2_antinu.h" #endif #ifndef __NO_GGM__ // GGM CC1ppip #include "GGM_CC1ppip_Evt_1DQ2_nu.h" #include "GGM_CC1ppip_XSec_1DEnu_nu.h" #endif #ifndef __NO_MiniBooNE__ // MiniBooNE CCQE #include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_nu.h" // MiniBooNE CC1pi+ 1D #include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h" #include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h" #include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h" #include "MiniBooNE_CC1pip_XSec_1DTu_nu.h" // MiniBooNE CC1pi+ 2D #include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h" // MiniBooNE CC1pi0 #include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h" #include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h" #include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h" // MiniBooNE NC1pi0 //#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h" // MiniBooNE NCEL #include "MiniBooNE_NCEL_XSec_Treco_nu.h" #endif #ifndef __NO_MINERvA__ // MINERvA CCQE #include "MINERvA_CCQE_XSec_1DQ2_antinu.h" #include "MINERvA_CCQE_XSec_1DQ2_joint.h" #include "MINERvA_CCQE_XSec_1DQ2_nu.h" // MINERvA CC0pi #include "MINERvA_CC0pi_XSec_1DEe_nue.h" #include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h" #include "MINERvA_CC0pi_XSec_1DQ2_nue.h" #include "MINERvA_CC0pi_XSec_1DThetae_nue.h" // 2018 MINERvA CC0pi STV #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" // 2018 MINERvA CC0pi 2D #include "MINERvA_CC0pi_XSec_2D_nu.h" #include "MINERvA_CC0pi_XSec_1D_2018_nu.h" // 2018 MINERvA CC0pi 2D antinu #include "MINERvA_CC0pi_XSec_2D_antinu.h" // MINERvA CC1pi+ #include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h" #include "MINERvA_CC1pip_XSec_1DTpi_nu.h" #include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h" #include "MINERvA_CC1pip_XSec_1Dth_nu.h" // 2017 data update #include "MINERvA_CC1pip_XSec_1D_2017Update.h" // MINERvA CCNpi+ #include "MINERvA_CCNpip_XSec_1DEnu_nu.h" #include "MINERvA_CCNpip_XSec_1DQ2_nu.h" #include "MINERvA_CCNpip_XSec_1DTpi_nu.h" #include "MINERvA_CCNpip_XSec_1Dpmu_nu.h" #include "MINERvA_CCNpip_XSec_1Dth_nu.h" #include "MINERvA_CCNpip_XSec_1Dthmu_nu.h" // MINERvA CC1pi0 #include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h" #include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h" #include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dth_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h" // MINERvA CC1pi0 neutrino #include "MINERvA_CC1pi0_XSec_1D_nu.h" // MINERvA CCINC #include "MINERvA_CCinc_XSec_1DEnu_ratio.h" #include "MINERvA_CCinc_XSec_1Dx_ratio.h" #include "MINERvA_CCinc_XSec_2DEavq3_nu.h" // MINERvA CCDIS #include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" #include "MINERvA_CCDIS_XSec_1Dx_ratio.h" // MINERvA CCCOH pion #include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h" #include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" #include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" #include "MINERvA_CCCOHPI_XSec_joint.h" #include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" #include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h" #endif #ifndef __NO_T2K__ // T2K CC0pi 2016 #include "T2K_CC0pi_XSec_2DPcos_nu.h" // T2K CC-inclusive with full acceptance 2018 #include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h" // T2K STV CC0pi 2018 #include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h" #include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h" #include "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h" // T2K CC1pi+ on CH -#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" -#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h" +#include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" +#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" +#include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h" +#include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h" +//#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" +//#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" +//#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" // T2K CC1pi+ on H2O #include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h" #include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" #endif #ifndef __NO_SciBooNE__ // SciBooNE COH studies #include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dpmu_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h" #include "SciBooNE_CCCOH_STOP_NTrks_nu.h" #endif #ifndef __NO_K2K__ // K2K NC1pi0 #include "K2K_NC1pi0_Evt_1Dppi0_nu.h" #endif // MC Studies #include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h" #include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h" #include "MCStudy_CCQEHistograms.h" #include "GenericFlux_Tester.h" #include "GenericFlux_Vectors.h" #include "ElectronFlux_FlatTree.h" #include "ElectronScattering_DurhamData.h" #include "MCStudy_KaonPreSelection.h" #include "MCStudy_MuonValidation.h" #include "OfficialNIWGPlots.h" #include "T2K2017_FakeData.h" #include "Simple_Osc.h" #include "Smear_SVDUnfold_Propagation_Osc.h" #include "FitWeight.h" #include "NuisConfig.h" #include "NuisKey.h" #ifdef __USE_DYNSAMPLES__ #include "TRegexp.h" #include // linux #include DynamicSampleFactory::DynamicSampleFactory() : NSamples(0), NManifests(0) { LoadPlugins(); QLOG(FIT, "Loaded " << NSamples << " from " << NManifests << " shared object libraries."); } DynamicSampleFactory* DynamicSampleFactory::glblDSF = NULL; DynamicSampleFactory::PluginManifest::~PluginManifest() { for (size_t i_it = 0; i_it < Instances.size(); ++i_it) { (*(DSF_DestroySample))(Instances[i_it]); } } std::string EnsureTrailingSlash(std::string const& inp) { if (!inp.length()) { return "/"; } if (inp[inp.length() - 1] == '/') { return inp; } return inp + "/"; } void DynamicSampleFactory::LoadPlugins() { std::vector SearchDirectories; if (Config::HasPar("dynamic_sample.path")) { SearchDirectories = GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":"); } char const* envPath = getenv("NUISANCE_DS_PATH"); if (envPath) { std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":"); for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) { SearchDirectories.push_back(envPaths[ep_it]); } } if (!SearchDirectories.size()) { char const* pwdPath = getenv("PWD"); if (pwdPath) { SearchDirectories.push_back(pwdPath); } } for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) { std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]); QLOG(FIT, "Searching for dynamic sample manifests in: " << dirpath); Ssiz_t len = 0; DIR* dir; struct dirent* ent; dir = opendir(dirpath.c_str()); if (dir != NULL) { TRegexp matchExp("*.so", true); while ((ent = readdir(dir)) != NULL) { if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) { QLOG(FIT, "\tFound shared object: " << ent->d_name << " checking for relevant methods..."); void* dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); char const* dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tDL Load Error: " << dlerr); continue; } PluginManifest plgManif; plgManif.dllib = dlobj; plgManif.soloc = (dirpath + ent->d_name); plgManif.DSF_NSamples = reinterpret_cast(dlsym(dlobj, "DSF_NSamples")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSampleName = reinterpret_cast( dlsym(dlobj, "DSF_GetSampleName")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSample = reinterpret_cast( dlsym(dlobj, "DSF_GetSample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_DestroySample = reinterpret_cast( dlsym(dlobj, "DSF_DestroySample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "Failed to load symbol \"DSF_DestroySample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.NSamples = (*(plgManif.DSF_NSamples))(); QLOG(FIT, "\tSuccessfully loaded dynamic sample manifest: " << plgManif.soloc << ". Contains " << plgManif.NSamples << " samples."); for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) { char const* smp_name = (*(plgManif.DSF_GetSampleName))(smp_it); if (!smp_name) { THROW("Could not load sample " << smp_it << " / " << plgManif.NSamples << " from " << plgManif.soloc); } if (Samples.count(smp_name)) { ERROR(WRN, "Already loaded a sample named: \"" << smp_name << "\". cannot load duplciates. This " "sample will be skipped."); continue; } plgManif.SamplesProvided.push_back(smp_name); Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it); QLOG(FIT, "\t\t" << smp_name); } if (plgManif.SamplesProvided.size()) { Manifests[plgManif.soloc] = plgManif; NSamples += plgManif.SamplesProvided.size(); NManifests++; } else { dlclose(dlobj); } } } closedir(dir); } else { ERROR(WRN, "Tried to open non-existant directory."); } } } DynamicSampleFactory& DynamicSampleFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSampleFactory(); } return *glblDSF; } void DynamicSampleFactory::Print() { std::map > ManifestSamples; for (std::map >::iterator smp_it = Samples.begin(); smp_it != Samples.end(); ++smp_it) { if (!ManifestSamples.count(smp_it->second.first)) { ManifestSamples[smp_it->second.first] = std::vector(); } ManifestSamples[smp_it->second.first].push_back(smp_it->first); } QLOG(FIT, "Dynamic sample manifest: "); for (std::map >::iterator m_it = ManifestSamples.begin(); m_it != ManifestSamples.end(); ++m_it) { QLOG(FIT, "\tLibrary " << m_it->first << " contains: "); for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) { QLOG(FIT, "\t\t" << m_it->second[s_it]); } } } bool DynamicSampleFactory::HasSample(std::string const& name) { return Samples.count(name); } bool DynamicSampleFactory::HasSample(nuiskey& samplekey) { return HasSample(samplekey.GetS("name")); } MeasurementBase* DynamicSampleFactory::CreateSample(nuiskey& samplekey) { if (!HasSample(samplekey)) { ERROR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name") << "\"."); return NULL; } std::pair sample = Samples[samplekey.GetS("name")]; QLOG(SAM, "\tLoading sample " << sample.second << " from " << sample.first); return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey); } DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); } #endif //! Functions to make it easier for samples to be created and handled. namespace SampleUtils { //! Create a given sample given its name, file, type, fakdata(fkdt) file and the //! current rw engine and push it back into the list fChain. MeasurementBase* CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight* rw) { nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", name); samplekey.Set("input", file); samplekey.Set("type", type); return CreateSample(samplekey); } MeasurementBase* CreateSample(nuiskey samplekey) { #ifdef __USE_DYNSAMPLES__ if (DynamicSampleFactory::Get().HasSample(samplekey)) { QLOG(SAM, "Instantiating dynamic sample..."); MeasurementBase* ds = DynamicSampleFactory::Get().CreateSample(samplekey); if (ds) { QLOG(SAM, "Done."); return ds; } THROW("Failed to instantiate dynamic sample."); } #endif FitWeight* rw = FitBase::GetRW(); std::string name = samplekey.GetS("name"); std::string file = samplekey.GetS("input"); std::string type = samplekey.GetS("type"); std::string fkdt = ""; /* ANL CCQE Samples */ #ifndef __NO_ANL__ if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) { return (new ANL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) { return (new ANL_CCQE_Evt_1DQ2_nu(samplekey)); /* ANL CC1ppip samples */ } else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) { return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") || !name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) { return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) { return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) { return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* ANL CC1npip sample */ } else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") || !name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) { return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC1pi0 sample */ } else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") || !name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) { return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey)); /* ANL NC1npip sample */ } else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) { return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey)); /* ANL NC1ppim sample */ } else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) { return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) { return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC2pi sample */ } else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey)); /* ArgoNeut Samples */ } else #endif #ifndef __NO_ArgoNeuT__ if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu(samplekey)); /* BNL Samples */ } else #endif #ifndef __NO_BNL__ if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) { return (new BNL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) { return (new BNL_CCQE_Evt_1DQ2_nu(samplekey)); /* BNL CC1ppip samples */ } else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) { return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") || !name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) { return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* BNL CC1npip samples */ } else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") || !name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) { return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) { return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey)); /* BNL CC1pi0 samples */ } else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) { return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) { return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey)); /* FNAL Samples */ } else #endif #ifndef __NO_FNAL__ if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) { return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppip */ } else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) { return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) { return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) { return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppim */ } else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) { return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey)); /* BEBC Samples */ } else #endif #ifndef __NO_BEBC__ if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) { return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey)); /* BEBC CC1ppip samples */ } else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) { return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) { return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npip samples */ } else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) { return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) { return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1pi0 samples */ } else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) { return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) { return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npim samples */ } else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) { return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) { return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey)); /* BEBC CC1ppim samples */ } else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) { return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) { return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey)); /* GGM CC1ppip samples */ } else #endif #ifndef __NO_GGM__ if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) { return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) { return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey)); /* MiniBooNE Samples */ /* CCQE */ } else #endif #ifndef __NO_MiniBooNE__ if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey)); /* MiniBooNE CC1pi+ */ // 1D } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey)); // 2D } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey)); /* MiniBooNE CC1pi0 */ } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey)); /* MiniBooNE NCEL */ } else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) { return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey)); /* MINERvA Samples */ } else #endif #ifndef __NO_MINERvA__ if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) { return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) { return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) { return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) { return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new MINERvA_CC0pinp_STV_XSec_1D_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) { return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey)); // Dan Ruterbories measurements of late 2018 } else if ( !name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) { return (new MINERvA_CC0pi_XSec_2D_nu(samplekey)); } else if ( !name.compare("MINERvA_CC0pi_XSec_1Dpt_nu") || !name.compare("MINERvA_CC0pi_XSec_1Dpz_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2QE_nu") || !name.compare("MINERvA_CC0pi_XSec_1DEnuQE_nu")) { return (new MINERvA_CC0pi_XSec_1D_2018_nu(samplekey)); // C. Patrick's early 2018 measurements } else if ( !name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) { return (new MINERvA_CC0pi_XSec_2D_antinu(samplekey)); /* CC1pi+ */ // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey)); // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) { return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey)); /* CCNpi+ */ } else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") || !name.compare( "MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) { return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) { return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey)); /* MINERvA CC1pi0 anti-nu */ // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) { return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) { return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) { return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey)); // MINERvA CC1pi0 nu } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) { return (new MINERvA_CC1pi0_XSec_1D_nu(samplekey)); /* CCINC */ } else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) { return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey)); /* CCDIS */ } else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey)); /* CC-COH */ } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); /* T2K Samples */ } else #endif #ifndef __NO_T2K__ if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu") || !name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") || !name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) { return (new T2K_CC0pi_XSec_2DPcos_nu(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CC0pi_XSec_2DPcos_nu_nonuniform(samplekey)); } else if (!name.compare("T2K_CCinc_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CCinc_XSec_2DPcos_nu_nonuniform(samplekey)); /* T2K CC1pi+ CH samples */ // Comment these out for now because we don't have the proper data - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dpmu_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dpmu_nu(samplekey)); + } else if (!name.compare("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu")) { + return (new T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey)); + } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { + return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(samplekey)); + + } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { + return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(samplekey)); + } else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) { - return (new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt)); + return (new T2K_CC1pip_CH_XSec_1DQ2_nu(samplekey)); - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); + } else if (!name.compare("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu")) { + return (new T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(samplekey)); - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt)); + } else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) { + return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey)); - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(file, rw, type, fkdt)); + // Maybe something for the future: were in Raquel's thesis + //} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) { + //return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt)); + //} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) { + //return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt)); - } else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) { - return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt)); + //} else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) { + //return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt)); /* T2K CC1pi+ H2O samples */ } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey)); /* T2K CC0pi + np CH samples */ } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddphit_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddat_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddat_nu(samplekey)); } else if (!name.compare("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform")) { return (new T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfp_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfp_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfa_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfa_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfip_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfip_nu(samplekey)); // SciBooNE COH studies } else #endif #ifndef __NO_SciBooNE__ if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) { return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) { return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dpmu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) { return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey)); /* K2K Samples */ /* NC1pi0 */ } else #endif #ifndef __NO_K2K__ if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) { return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey)); /* Fake Studies */ } else #endif if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("GenericFlux_") != std::string::npos) { return (new GenericFlux_Tester(name, file, rw, type, fkdt)); } else if (name.find("GenericVectors_") != std::string::npos) { return (new GenericFlux_Vectors(name, file, rw, type, fkdt)); } else if (!name.compare("T2K2017_FakeData")) { return (new T2K2017_FakeData(samplekey)); } else if (!name.compare("MCStudy_CCQE")) { return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt)); } else if (!name.compare("ElectronFlux_FlatTree")) { return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt)); } else if (name.find("ElectronData_") != std::string::npos) { return new ElectronScattering_DurhamData(samplekey); } else if (name.find("MuonValidation_") != std::string::npos) { return (new MCStudy_MuonValidation(name, file, rw, type, fkdt)); } else if (!name.compare("NIWGOfficialPlots")) { return (new OfficialNIWGPlots(samplekey)); } else if (!name.compare("Simple_Osc")) { return (new Simple_Osc(samplekey)); } else if (!name.compare("Smear_SVDUnfold_Propagation_Osc")) { return (new Smear_SVDUnfold_Propagation_Osc(samplekey)); } else { THROW("Error: No such sample: " << name << std::endl); } // Return NULL if no sample loaded. return NULL; } } diff --git a/src/T2K/CMakeLists.txt b/src/T2K/CMakeLists.txt index db91ee4..339a1f3 100644 --- a/src/T2K/CMakeLists.txt +++ b/src/T2K/CMakeLists.txt @@ -1,99 +1,102 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ set(IMPLFILES T2K_CC0pi_XSec_2DPcos_nu.cxx T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx -T2K_CC1pip_CH_XSec_1DQ2_nu.cxx -T2K_CC1pip_CH_XSec_1DWrec_nu.cxx -T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx + +T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx T2K_CC1pip_CH_XSec_1Dppi_nu.cxx -T2K_CC1pip_CH_XSec_1Dq3_nu.cxx -T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx -T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx +T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx +T2K_CC1pip_CH_XSec_1DQ2_nu.cxx +T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx +T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx + T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx + T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx T2K_SignalDef.cxx ) set(HEADERFILES T2K_CC0pi_XSec_2DPcos_nu.h T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx -T2K_CC1pip_CH_XSec_1DQ2_nu.h -T2K_CC1pip_CH_XSec_1DWrec_nu.h -T2K_CC1pip_CH_XSec_1Dpmu_nu.h + +T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h T2K_CC1pip_CH_XSec_1Dppi_nu.h -T2K_CC1pip_CH_XSec_1Dq3_nu.h -T2K_CC1pip_CH_XSec_1Dthmupi_nu.h T2K_CC1pip_CH_XSec_1Dthpi_nu.h -T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h +T2K_CC1pip_CH_XSec_1Dthmupi_nu.h +T2K_CC1pip_CH_XSec_1DQ2_nu.h +T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h +T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h + T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h T2K_CC1pip_H2O_XSec_1Dcospi_nu.h T2K_CC1pip_H2O_XSec_1Dpmu_nu.h T2K_CC1pip_H2O_XSec_1Dppi_nu.h T2K_CC0pinp_STV_XSec_1Ddpt_nu.h T2K_CC0pinp_STV_XSec_1Ddphit_nu.h T2K_CC0pinp_STV_XSec_1Ddat_nu.h T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h T2K_SignalDef.h ) set(LIBNAME expT2K) if(CMAKE_BUILD_TYPE MATCHES DEBUG) add_library(${LIBNAME} STATIC ${IMPLFILES}) else(CMAKE_BUILD_TYPE MATCHES RELEASE) add_library(${LIBNAME} SHARED ${IMPLFILES}) endif() include_directories(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? #install(FILES ${HEADERFILES} DESTINATION include) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx new file mode 100644 index 0000000..343fed5 --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx @@ -0,0 +1,98 @@ +#include + +#include "T2K_SignalDef.h" + +#include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h" + +// The constructor +T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("#phi_{Adler} (radians)"); + fSettings.SetYTitle("d#sigma/d#phi_{Adler} (cm^{2}/rad/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.rootout.root", "Phi_Adler"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.rootout.root", "Phi_AdlerCov"); + + SetShapeCovar(); + fDataHist->Scale(1E-38); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + + +void T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::FillEventVariables(FitEvent *event) { + + if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; + + // Reconstruct the neutrino + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; + double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip)*1000.; + // Now make the reconstructed neutrino + // Has the same direction as the predicted neutrino + TLorentzVector PnuReco(Pnu.Vect().Unit().X()*rEnu, Pnu.Vect().Unit().Y()*rEnu, Pnu.Vect().Unit().Z()*rEnu, rEnu); + + // Reconstruct the initial state assuming still nucleon + TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton*1000.); + // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed neutrino instead of true neutrino{ + + // Get the "resonance" lorentz vector (pion proton system) reconstructed from the variables + TLorentzVector Pres = PnuReco+Pinit-Pmu; + // Boost the particles into the resonance rest frame so we can define the x,y,z axis + PnuReco.Boost(Pres.BoostVector()); + Pmu.Boost(-Pres.BoostVector()); + Ppip.Boost(-Pres.BoostVector()); + + // The z-axis is defined as the axis of three-momentum transfer, \vec{k} + // Also unit normalise the axis + TVector3 zAxis = (PnuReco.Vect()-Pmu.Vect())*(1.0/((PnuReco.Vect()-Pmu.Vect()).Mag())); + + // The y-axis is then defined perpendicular to z and muon momentum in the resonance frame + TVector3 yAxis = PnuReco.Vect().Cross(Pmu.Vect()); + yAxis *= 1.0/double(yAxis.Mag()); + + // And the x-axis is then simply perpendicular to z and x + TVector3 xAxis = yAxis.Cross(zAxis); + xAxis *= 1.0/double(xAxis.Mag()); + + double x = Ppip.Vect().Dot(xAxis); + double y = Ppip.Vect().Dot(yAxis); + //double z = Ppi.Vect().Dot(zAxis); + + double newphi = atan2(y, x); + // Convert negative angles to positive + //if (newphi < 0.0) newphi += 360.0; + + fXVar = newphi; + + return; +}; + +//******************************************************************** +bool T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::isSignal(FitEvent *event) { +//******************************************************************** + return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); +} + diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h new file mode 100644 index 0000000..80bd168 --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h @@ -0,0 +1,16 @@ +#ifndef T2K_CC1PIP_CH_XSEC_1DADLERPHI_NU_H_SEEN +#define T2K_CC1PIP_CH_XSEC_1DADLERPHI_NU_H_SEEN + +#include "Measurement1D.h" + +class T2K_CC1pip_CH_XSec_1DAdlerPhi_nu : public Measurement1D { +public: + T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(nuiskey samplekey); + virtual ~T2K_CC1pip_CH_XSec_1DAdlerPhi_nu() {}; + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + +}; + +#endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx new file mode 100644 index 0000000..d21785a --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx @@ -0,0 +1,86 @@ +#include + +#include "T2K_SignalDef.h" + +#include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h" + +// The constructor +T2K_CC1pip_CH_XSec_1DCosThAdler_nu::T2K_CC1pip_CH_XSec_1DCosThAdler_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_1DCosThAdler_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DCosThAdler_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("cos#theta_{Adler} (radians)"); + fSettings.SetYTitle("d#sigma/dcos#theta_{Adler} (cm^{2}/1/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.rootout.root", "Theta_Adler"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.rootout.root", "Theta_AdlerCov"); + + SetShapeCovar(); + fDataHist->Scale(1E-38); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + + +void T2K_CC1pip_CH_XSec_1DCosThAdler_nu::FillEventVariables(FitEvent *event) { + + if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; + + // Essentially the same as FitUtils::CosThAdler but uses the reconsturcted neutrino instead of the true neutrino + // Reconstruct the neutrino + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; + double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip)*1000.; + // Now make the reconstructed neutrino + // Has the same direction as the predicted neutrino + TLorentzVector PnuReco(Pnu.Vect().Unit().X()*rEnu, Pnu.Vect().Unit().Y()*rEnu, Pnu.Vect().Unit().Z()*rEnu, rEnu); + + // Reconstruct the initial state assuming still nucleon + TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton*1000.); + // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed neutrino instead of true neutrino{ + + // Get the "resonance" lorentz vector (pion proton system) reconstructed from the variables + TLorentzVector Pres = PnuReco+Pinit-Pmu; + // Boost the particles into the resonance rest frame so we can define the x,y,z axis + PnuReco.Boost(Pres.BoostVector()); + Pmu.Boost(-Pres.BoostVector()); + Ppip.Boost(-Pres.BoostVector()); + + // The z-axis is defined as the axis of three-momentum transfer, \vec{k} + // Also unit normalise the axis + TVector3 zAxis = (PnuReco.Vect()-Pmu.Vect())*(1.0/((PnuReco.Vect()-Pmu.Vect()).Mag())); + + // Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle + double costhAdler = cos(Ppip.Vect().Angle(zAxis)); + + fXVar = costhAdler; + + return; +}; + +//******************************************************************** +bool T2K_CC1pip_CH_XSec_1DCosThAdler_nu::isSignal(FitEvent *event) { +//******************************************************************** + return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); +} + diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h new file mode 100644 index 0000000..4b29e2a --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h @@ -0,0 +1,16 @@ +#ifndef T2K_CC1PIP_CH_XSEC_1DCOSTHADLER_NU_H_SEEN +#define T2K_CC1PIP_CH_XSEC_1DCOSTHADLER_NU_H_SEEN + +#include "Measurement1D.h" + +class T2K_CC1pip_CH_XSec_1DCosThAdler_nu : public Measurement1D { +public: + T2K_CC1pip_CH_XSec_1DCosThAdler_nu(nuiskey samplekey); + virtual ~T2K_CC1pip_CH_XSec_1DCosThAdler_nu() {}; + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + +}; + +#endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx index ee52182..d3ab741 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx @@ -1,201 +1,118 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" // The constructor -T2K_CC1pip_CH_XSec_1DQ2_nu::T2K_CC1pip_CH_XSec_1DQ2_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - - // Here we can give either MB (kMB), extended MB (keMB) or Delta (kDelta) - if (type.find("eMB") != std::string::npos) { - fT2KSampleType = keMB; - fName = "T2K_CC1pip_CH_XSec_1DQ2eMB_nu"; - fPlotTitles = "; Q^{2}_{eMB} (GeV^{2}); d#sigma/dQ^{2}_{eMB} (cm^{2}/GeV^{2}/nucleon)"; - } else if (type.find("MB") != std::string::npos) { - fT2KSampleType = kMB; - fName = "T2K_CC1pip_CH_XSec_1DQ2MB_nu"; - fPlotTitles = "; Q^{2}_{MB} (GeV^{2}); d#sigma/dQ^{2}_{MB} (cm^{2}/GeV^{2}/nucleon)"; - } else if (type.find("Delta") != std::string::npos) { - fT2KSampleType = kDelta; - fName = "T2K_CC1pip_CH_XSec_1DQ2delta_nu"; - fPlotTitles = "; Q^{2}_{#Delta} (GeV^{2}); d#sigma/dQ^{2}_{#Delta} (cm^{2}/GeV^{2}/nucleon)"; - } else { - LOG(SAM) << "Found no specified type, using MiniBooNE E_nu/Q2 definition" << std::endl; - fT2KSampleType = kMB; - fName = "T2K_CC1pip_CH_XSec_1DQ2MB_nu"; - fPlotTitles = "; Q^{2}_{MB} (GeV^{2}); d#sigma/dQ^{2}_{MB} (cm^{2}/GeV^{2}/nucleon)"; - } - - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - //type = keMB; - //type = kDelta; - - if (fT2KSampleType == kMB) { - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_MB.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_MB.root"); - } else if (fT2KSampleType == keMB) { - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_extendedMB.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_extendedMB.root"); - } else if (fT2KSampleType == kDelta) { - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_Delta.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_Delta.root"); - } else { - ERR(FTL) << "No data type set for " << fName << std::endl; - ERR(FTL) << __FILE__ << ":" << __LINE__ << std::endl; - exit(-1); - } - SetShapeCovar(); - - this->SetupDefaultHist(); +T2K_CC1pip_CH_XSec_1DQ2_nu::T2K_CC1pip_CH_XSec_1DQ2_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_1DQ2_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DQ2_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("Q^{2} (GeV/c)^{2}"); + fSettings.SetYTitle("d#sigma/dQ^{2} (cm^{2}/GeV^{2}/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.rootout.root", "Q2"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.rootout.root", "Q2Cov"); - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1DQ2_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - const int nPoints = dataCopy->GetNbinsX()-6; - LOG(DEB) << nPoints << std::endl; - double *binEdges = new double[nPoints+1]; - for (int i = 0; i < nPoints+1; i++) { - binEdges[i] = dataCopy->GetBinLowEdge(i+1); - } - - for (int i = 0; i < nPoints+1; i++) { - LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), nPoints, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); - - dataFile->Close(); - -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1DQ2_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - LOG(DEB) << nBinsX << std::endl; - LOG(DEB) << fDataHist->GetNbinsX() << std::endl; - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - this->fFullCovar = new TMatrixDSym(nBinsX-7); + SetShapeCovar(); + fDataHist->Scale(1E-38); - // First two entries are BS - // Last entry is BS - for (int i = 0; i < nBinsX-7; i++) { - for (int j = 0; j < nBinsY-7; j++) { - (*this->fFullCovar)(i, j) = covarMatrix->GetBinContent(i+3, j+3); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i << ", " << j << ") = " << (*this->fFullCovar)(i,j) << std::endl; - } - } //should now have set covariance, I hope + // Final setup --------------------------------------------------- + FinaliseMeasurement(); - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - return; }; void T2K_CC1pip_CH_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double q2 = -999; - switch(fT2KSampleType) { + //switch(fT2KSampleType) { // First int refers to how we reconstruct Enu // 0 uses true neutrino energy (not used here but common for other analyses when they unfold to true Enu from reconstructed Enu) // 1 uses "extended MiniBooNE" method // 2 uses "MiniBooNE reconstructed" method // 3 uses Delta resonance mass for reconstruction // // The last bool refers to if pion directional information was used // // Use MiniBooNE reconstructed Enu; uses Michel tag so no pion direction information - case kMB: - q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 2, false); - break; + //case kMB: + //q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 2, false); + //break; // Use Extended MiniBooNE reconstructed Enu // Needs pion information to reconstruct so bool is true (did not include Michel e tag) - case keMB: + //case keMB: q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 1, true); - break; + //break; // Use Delta resonance reconstructed Enu // Uses Michel electron so don't have pion directional information - case kDelta: - q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 3, false); - break; - } + //case kDelta: + //q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 3, false); + //break; + //} fXVar = q2; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DQ2_nu::isSignal(FitEvent *event) { //******************************************************************** // Warning: The CH analysis has different signal definition to the H2O analysis! // Often to do with the Michel tag - switch(fT2KSampleType) { + //switch(fT2KSampleType) { // Using MiniBooNE formula for Enu reconstruction on the Q2 variable // Does have Michel e tag, set bool to true! - case kMB: - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); - break; + //case kMB: + //return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); + //break; // Using extended MiniBooNE formula for Enu reconstruction on the Q2 variable // Does not have Michel e tag because we need directional information to reconstruct Q2 - case keMB: + //case keMB: return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); - break; + //break; // Using Delta resonance for Enu reconstruction on the Q2 variable // Does have Michel e tag, bool to true - case kDelta: - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); - break; - } + //case kDelta: + //return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); + //break; + //} // Default to return false - return false; + //return false; } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.h index 9091cc4..dbc62c6 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.h +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.h @@ -1,25 +1,18 @@ #ifndef T2K_CC1PIP_CH_XSEC_1DQ2_NU_H_SEEN #define T2K_CC1PIP_CH_XSEC_1DQ2_NU_H_SEEN #include "Measurement1D.h" -enum T2K_CC1pip_type {kMB = 0, keMB, kDelta}; +//enum T2K_CC1pip_type {kMB = 0, keMB, kDelta}; class T2K_CC1pip_CH_XSec_1DQ2_nu : public Measurement1D { public: - T2K_CC1pip_CH_XSec_1DQ2_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); + T2K_CC1pip_CH_XSec_1DQ2_nu(nuiskey samplekey); virtual ~T2K_CC1pip_CH_XSec_1DQ2_nu() {}; - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - private: - T2K_CC1pip_type fT2KSampleType; - }; #endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx deleted file mode 100644 index c31b6aa..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx +++ /dev/null @@ -1,116 +0,0 @@ -#include - -#include "T2K_SignalDef.h" - -#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" - -// The constructor -T2K_CC1pip_CH_XSec_1DWrec_nu::T2K_CC1pip_CH_XSec_1DWrec_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - fName = "T2K_CC1pip_CH_XSec_1DWrec_nu"; - fPlotTitles = "; W_{rec} (GeV/c); d#sigma/dW_{rec} (cm^{2}/(GeV/c^{2})/nucleon)"; - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/W.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/W.root"); - SetShapeCovar(); - this->SetupDefaultHist(); - - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1DWrec_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the first and last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - LOG(DEB) << dataCopy->GetNbinsX() << std::endl; - const int dataPoints = dataCopy->GetNbinsX()-2; - double *binEdges = new double[dataPoints+1]; - // Want to skip the first bin here - for (int i = 0; i < dataPoints+1; i++) { - binEdges[i] = dataCopy->GetBinLowEdge(i+2); - } - - for (int i = 0; i < dataPoints+1; i++) { - LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataPoints, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+2)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+2)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - //fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); - - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1DWrec_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - const int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - const int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - LOG(DEB) << nBinsX << std::endl; - LOG(DEB) << fDataHist->GetNbinsX() << std::endl; - this->fFullCovar = new TMatrixDSym(nBinsX-3); - - for (int i = 0; i < nBinsX-3; i++) { - for (int j = 0; j < nBinsY-3; j++) { - (*this->fFullCovar)(i, j) = covarMatrix->GetBinContent(i+4, j+4); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i << ", " << j << ") = " << (*this->fFullCovar)(i,j) << std::endl; - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - - return; -}; - - -void T2K_CC1pip_CH_XSec_1DWrec_nu::FillEventVariables(FitEvent *event) { - - if (event->NumFSParticle(13) == 0 || - event->NumFSParticle(211) == 0) - return; - - TLorentzVector Pnu = event->GetNeutrinoIn()->fP; - TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; - TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; - - double Wrec = FitUtils::WrecCC1pip_T2K_MB(Pnu, Pmu, Ppip); - - fXVar = Wrec; - - return; -}; - -//******************************************************************** -bool T2K_CC1pip_CH_XSec_1DWrec_nu::isSignal(FitEvent *event) { -//******************************************************************** -// This sample includes the Michel e tag so do not have to cut into the pion phase space - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); -} - diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.h deleted file mode 100644 index 0c24978..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef T2K_CC1PIP_CH_XSEC_1DWREC_NU_H_SEEN -#define T2K_CC1PIP_CH_XSEC_1DWREC_NU_H_SEEN - -#include "Measurement1D.h" - -class T2K_CC1pip_CH_XSec_1DWrec_nu : public Measurement1D { -public: - T2K_CC1pip_CH_XSec_1DWrec_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); - virtual ~T2K_CC1pip_CH_XSec_1DWrec_nu() {}; - - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - - void FillEventVariables(FitEvent *event); - bool isSignal(FitEvent *event); - - private: -}; - -#endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx index f787b76..898329c 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx @@ -1,177 +1,90 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" //******************************************************************** T2K_CC1pip_CH_XSec_1Dppi_nu::T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_1Dppi_nu sample. \n" \ "Target: CH \n" \ - "Flux: T2k Forward Horn Current nue + nuebar \n" \ - "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n"; + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dppi_nu"); fSettings.SetDescription(descrip); fSettings.SetXTitle("p_{#pi} (GeV/c)"); - fSettings.SetYTitle("d#sigma/dW_{rec} (cm^{2}/(GeV/c)/nucleon)"); + fSettings.SetYTitle("d#sigma/dp_{#pi} (cm^{2}/(GeV/c)/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- - if (fSettings.GetS("type").find("Michel") != std::string::npos) { - useMichel = true; - fName += "_Michel"; - SetDataValues(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi.root"); - SetCovarMatrix(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi.root"); - } else { - useMichel = false; - // fName += "_kin"; - SetDataValues(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi_noME.root"); - SetCovarMatrix(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi_noME.root"); - } + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.rootout.root", "Momentum_pion"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.rootout.root", "Momentum_pionCov"); SetShapeCovar(); + fDataHist->Scale(1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dppi_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - LOG(DEB) << dataCopy->GetNbinsX() << std::endl; - - double *binEdges = new double[dataCopy->GetNbinsX() - 1]; - for (int i = 0; i < dataCopy->GetNbinsX() - 1; i++) { - binEdges[i] = dataCopy->GetBinLowEdge(i + 1); - } - binEdges[dataCopy->GetNbinsX() - 1] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()); - - fDataHist = new TH1D((fName + "_data").c_str(), (fName + "_data" + fPlotTitles).c_str(), dataCopy->GetNbinsX() - 2, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i + 1, dataCopy->GetBinContent(i + 1) * 1E-38); - fDataHist->SetBinError(i + 1, dataCopy->GetBinError(i + 1) * 1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i + 1) << " " << fDataHist->GetBinContent(i + 1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_MC" + fPlotTitles).c_str()); - - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dppi_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - this->fFullCovar = new TMatrixDSym(nBinsX - 3); - - // First two entries are BS - // Last entry is BS - for (int i = 2; i < nBinsX-1; i++) { - for (int j = 2; j < nBinsY-1; j++) { - LOG(DEB) << "(" << i << ", " << j << ") = " << covarMatrix->GetBinContent(i + 1, j + 1) << std::endl; - (*this->fFullCovar)(i - 2, j - 2) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - return; -}; - - void T2K_CC1pip_CH_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double ppip = FitUtils::p(Ppip); fXVar = ppip; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dppi_nu::isSignal(FitEvent *event) { //******************************************************************** // This distribution uses a somewhat different signal definition so might as well implement it separately here - // If we use Michel tag sample we don't cut into the pion phase space, only the muon phase space - // The last bool refers to if we have Michel e or not - if (useMichel) { - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); - } else { // Custom signal definition if we aren't using Michel tag; cut on muon and cut only on pion angle - - // does the event pass the muon cut? - bool muonPass = SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); - - // If the event doesn't pass the muon cut return false - if (!muonPass) { - return false; - } - - // does the event pass the pion angle cut? - // we already know there's just one muon in the event if it passes muonPass so don't need to make an event loop rejection - // Need the neutrino four-vector to get the angle between pion and neutrino - TLorentzVector Pnu = event->PartInfo(0)->fP; - TLorentzVector Ppip; - for (unsigned int j = 2; j < event->Npart(); j++) { - if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fNEUTStatusCode != 0) continue; //move on if NOT ALIVE and NOT NORMAL - int PID = (event->PartInfo(j))->fPID; - if (PID == 211) { - Ppip = event->PartInfo(j)->fP; // Once the pion is found we can break - break; - } - } - - double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); - // Now check the angle of the pion - if (cos_th_pi <= 0.2) { - return false; - } else { - return true; - } + if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false; + + TLorentzVector Pnu = event->GetHMISParticle(14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; + + // If this event passes the criteria on particle counting, enforce the T2K + // ND280 phase space constraints + // Will be different if Michel tag sample is included or not + // Essentially, if there's a Michel tag we don't cut on the pion variables + + double p_mu = FitUtils::p(Pmu) * 1000; + double cos_th_mu = cos(FitUtils::th(Pnu, Pmu)); + double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); + + if (p_mu <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.2) { + return false; + } else { + return true; } - // Unnecessary default to false return false; } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.h index 8649e4a..07ce298 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.h +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.h @@ -1,22 +1,15 @@ #ifndef T2K_CC1PIP_CH_XSEC_1DPPI_NU_H_SEEN #define T2K_CC1PIP_CH_XSEC_1DPPI_NU_H_SEEN #include "Measurement1D.h" class T2K_CC1pip_CH_XSec_1Dppi_nu : public Measurement1D { public: T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey); virtual ~T2K_CC1pip_CH_XSec_1Dppi_nu() {}; - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - - private: - bool useMichel; }; #endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx deleted file mode 100644 index 6a34455..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx +++ /dev/null @@ -1,109 +0,0 @@ -#include - -#include "T2K_SignalDef.h" - -#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" - -// The constructor -T2K_CC1pip_CH_XSec_1Dq3_nu::T2K_CC1pip_CH_XSec_1Dq3_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - fName = "T2K_CC1pip_CH_XSec_1Dq3_nu"; - fPlotTitles = "; q_{3} (GeV/c); d#sigma/dq_{3} (cm^{2}/(GeV/c)/nucleon)"; - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q3.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q3.root"); - SetShapeCovar(); - this->SetupDefaultHist(); - - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dq3_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - double *binEdges = new double[dataCopy->GetNbinsX()-2]; - LOG(DEB) << dataCopy->GetNbinsX() << std::endl; - for (int i = 1; i < dataCopy->GetNbinsX()-1; i++) { - binEdges[i-1] = dataCopy->GetBinLowEdge(i+1); - LOG(DEB) << i-1 << " " << binEdges[i-1] << " from binLowEdge " << i+1 << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-2, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+2)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+2)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); - - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dq3_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - this->fFullCovar = new TMatrixDSym(nBinsX-3); - - // First two entries are BS - // Last entry is BS - for (int i = 2; i < nBinsX-1; i++) { - for (int j = 2; j < nBinsY-1; j++) { - (*this->fFullCovar)(i-2, j-2) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i-2 << ", " << j-2 << ") = " << (*this->fFullCovar)(i-2,j-2) << std::endl; - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - return; -}; - - -void T2K_CC1pip_CH_XSec_1Dq3_nu::FillEventVariables(FitEvent *event) { - - if (event->NumFSParticle(13) == 0 || - event->NumFSParticle(211) == 0) - return; - - TLorentzVector Pnu = event->GetNeutrinoIn()->fP; - TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; - TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; - - double q3 = FitUtils::q3_CC1pip_T2K(Pnu, Pmu, Ppip); - - fXVar = q3; - - return; -}; - -//******************************************************************** -bool T2K_CC1pip_CH_XSec_1Dq3_nu::isSignal(FitEvent *event) { -//******************************************************************** -// Has a Michel e sample in so no phase space cut on pion required - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); -} - diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.h deleted file mode 100644 index ff2227d..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef T2K_CC1PIP_CH_XSEC_1DQ3_NU_H_SEEN -#define T2K_CC1PIP_CH_XSEC_1DQ3_NU_H_SEEN - -#include "Measurement1D.h" - -class T2K_CC1pip_CH_XSec_1Dq3_nu : public Measurement1D { -public: - T2K_CC1pip_CH_XSec_1Dq3_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); - virtual ~T2K_CC1pip_CH_XSec_1Dq3_nu() {}; - - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - - void FillEventVariables(FitEvent *event); - bool isSignal(FitEvent *event); - - private: -}; - -#endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx index da6cc55..77d2048 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx @@ -1,109 +1,65 @@ #include #include "T2K_SignalDef.h" - #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" // The constructor -T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - fName = "T2K_CC1pip_CH_XSec_1Dthmupi_nu"; - fPlotTitles = "; #theta_{#pi,#mu} (radians); d#sigma/d#theta_{#pi} (cm^{2}/nucleon)"; - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapimu.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapimu.root"); +T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_1Dthmupi_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthmupi_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("#theta_{#mu,#pi} (radians)"); + fSettings.SetYTitle("d#sigma/d#theta_{#mu,#pi} (cm^{2}/radians/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.rootout.root", "Theta(pi,mu)(rads)"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.rootout.root", "Theta(pi,mu)(rads)Cov"); + SetShapeCovar(); - this->SetupDefaultHist(); + fDataHist->Scale(1E-38); - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); + FinaliseMeasurement(); }; -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthmupi_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the first and last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - LOG(DEB) << "dataCopy->GetNbinsX() = " << dataCopy->GetNbinsX() << std::endl; - double *binEdges = new double[dataCopy->GetNbinsX()]; - for (int i = 0; i < dataCopy->GetNbinsX(); i++) { - binEdges[i] = dataCopy->GetBinLowEdge(i+1); - } - binEdges[dataCopy->GetNbinsX()] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()+1); - - for (int i = 0; i < dataCopy->GetNbinsX()+5; i++) { - LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX(), binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX()+1; i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthmupi_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - this->fFullCovar = new TMatrixDSym(nBinsX-1); - - for (int i = 1; i < nBinsX; i++) { - for (int j = 1; j < nBinsY; j++) { - (*this->fFullCovar)(i-1, j-1) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i-1 << ", " << j-1 << ") = " << (*this->fFullCovar)(i-1,j-1) << std::endl; - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - dataFile->Close(); -}; - - void T2K_CC1pip_CH_XSec_1Dthmupi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double thmupi = FitUtils::th(Pmu, Ppip); fXVar = thmupi; + //std::cout << thmupi << std::endl; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthmupi_nu::isSignal(FitEvent *event) { //******************************************************************** -// This sample requires directional information on the pion so can't use Michel tag sample +// This distribution uses a somewhat different signal definition so might as well implement it separately here return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.h index 4627b31..033bec0 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.h +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.h @@ -1,21 +1,17 @@ #ifndef T2K_CC1PIP_CH_XSEC_1DTHMUPI_NU_H_SEEN #define T2K_CC1PIP_CH_XSEC_1DTHMUPI_NU_H_SEEN #include "Measurement1D.h" class T2K_CC1pip_CH_XSec_1Dthmupi_nu : public Measurement1D { public: - T2K_CC1pip_CH_XSec_1Dthmupi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); + T2K_CC1pip_CH_XSec_1Dthmupi_nu(nuiskey samplekey); virtual ~T2K_CC1pip_CH_XSec_1Dthmupi_nu() {}; - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); private: }; #endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx index b8be6c5..7473635 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx @@ -1,110 +1,83 @@ #include #include "T2K_SignalDef.h" - #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" // The constructor -T2K_CC1pip_CH_XSec_1Dthpi_nu::T2K_CC1pip_CH_XSec_1Dthpi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - fName = "T2K_CC1pip_CH_XSec_1Dthpi_nu"; - fPlotTitles = "; #theta_{#pi} (radians); d#sigma/d#theta_{#pi} (cm^{2}/nucleon)"; - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapi.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapi.root"); +T2K_CC1pip_CH_XSec_1Dthpi_nu::T2K_CC1pip_CH_XSec_1Dthpi_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_1Dthpi_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2K Forward Horn Current numu \n" \ + "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthpi_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("#theta_{#pi} (radians)"); + fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/radians/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.rootout.root", "Theta_pion"); + SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.rootout.root", "Theta_pionCov"); + SetShapeCovar(); - this->SetupDefaultHist(); - - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthpi_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the first and last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - LOG(DEB) << "dataCopy->GetNbinsX() = " << dataCopy->GetNbinsX() << std::endl; - double *binEdges = new double[dataCopy->GetNbinsX()-4]; - for (int i = 0; i < dataCopy->GetNbinsX()-4; i++) { - binEdges[i] = dataCopy->GetBinLowEdge(i+1); - } - binEdges[dataCopy->GetNbinsX()-4] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()-3); - - for (int i = 0; i < dataCopy->GetNbinsX(); i++) { - LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-4, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthpi_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); + fDataHist->Scale(1E-38); - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - this->fFullCovar = new TMatrixDSym(nBinsX-5); - - for (int i = 2; i < nBinsX-3; i++) { - for (int j = 2; j < nBinsY-3; j++) { - (*this->fFullCovar)(i-2, j-2) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i-2 << ", " << j-2 << ") = " << (*this->fFullCovar)(i-2,j-2) << std::endl; - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - dataFile->Close(); + FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1Dthpi_nu::FillEventVariables(FitEvent *event) { - if (event->NumFSParticle(13) == 0 || - event->NumFSParticle(211) == 0) - return; + if (event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; - TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double thpi = FitUtils::th(Pnu, Ppip); fXVar = thpi; - - return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthpi_nu::isSignal(FitEvent *event) { //******************************************************************** -// This sample uses directional info on the pion so Michel e tag sample can not be included -// i.e. we need reduce the pion variable phase space - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); +// This distribution uses a somewhat different signal definition so might as well implement it separately here + + if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false; + + TLorentzVector Pnu = event->GetHMISParticle(14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; + + // If this event passes the criteria on particle counting, enforce the T2K + // ND280 phase space constraints + // Will be different if Michel tag sample is included or not + // Essentially, if there's a Michel tag we don't cut on the pion variables + + double p_mu = FitUtils::p(Pmu) * 1000; + double p_pi = FitUtils::p(Ppip) * 1000; + double cos_th_mu = cos(FitUtils::th(Pnu, Pmu)); + double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); + + if (p_mu <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.0 || p_pi <= 200) { + return false; + } else { + return true; + } + + return false; } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.h index 5052bf7..b43e7d1 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.h +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.h @@ -1,21 +1,15 @@ #ifndef T2K_CC1PIP_CH_XSEC_1DTHPI_NU_H_SEEN #define T2K_CC1PIP_CH_XSEC_1DTHPI_NU_H_SEEN #include "Measurement1D.h" class T2K_CC1pip_CH_XSec_1Dthpi_nu : public Measurement1D { public: - T2K_CC1pip_CH_XSec_1Dthpi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); + T2K_CC1pip_CH_XSec_1Dthpi_nu(nuiskey samplekey); virtual ~T2K_CC1pip_CH_XSec_1Dthpi_nu() {}; - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - - private: }; #endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx deleted file mode 100644 index 1426e80..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx +++ /dev/null @@ -1,115 +0,0 @@ -#include - -#include "T2K_SignalDef.h" - -#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" - -// The constructor -T2K_CC1pip_CH_XSec_1Dthq3pi_nu::T2K_CC1pip_CH_XSec_1Dthq3pi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ - - fName = "T2K_CC1pip_CH_XSec_1Dthq3pi_nu"; - fPlotTitles = "; #theta_{q_{3},#pi} (radians); d#sigma/d#theta_{q_{3},#pi} (cm^{2}/(radian)/nucleon)"; - EnuMin = 0.; - EnuMax = 100.; - fIsDiag = false; - Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); - - this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/ThetaQ3Pi.root"); - this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/ThetaQ3Pi.root"); - SetShapeCovar(); - this->SetupDefaultHist(); - - this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::SetDataValues(std::string fileLocation) { - LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - // Don't want the last bin of dataCopy - TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - - double *binEdges = new double[dataCopy->GetNbinsX()-1]; - LOG(DEB) << dataCopy->GetNbinsX() << std::endl; - for (int i = 1; i < dataCopy->GetNbinsX()+1; i++) { - binEdges[i-1] = dataCopy->GetBinLowEdge(i); - LOG(DEB) << i-1 << " " << binEdges[i-1] << std::endl; - } - - for (int i = 0; i < dataCopy->GetNbinsX()+5; i++) { - LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; - } - - fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-1, binEdges); - - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); - fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); - LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; - } - - fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile - fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); - - - dataFile->Close(); -}; - -// Override this for now -// Should really have Measurement1D do this properly though -void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::SetCovarMatrix(std::string fileLocation) { - LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; - TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! - - TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); - - int nBinsX = covarMatrix->GetXaxis()->GetNbins(); - int nBinsY = covarMatrix->GetYaxis()->GetNbins(); - - if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - - // First bin is underflow, last bin is overflow - this->fFullCovar = new TMatrixDSym(nBinsX-2); - - // First two entries are BS - // Last entry is BS - for (int i = 1; i < nBinsX-1; i++) { - for (int j = 1; j < nBinsY-1; j++) { - (*this->fFullCovar)(i-1, j-1) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances - LOG(DEB) << "fFullCovar(" << i-1 << ", " << j-1 << ") = " << (*this->fFullCovar)(i-1,j-1) << std::endl; - } - } //should now have set covariance, I hope - - this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); - this->covar = StatUtils::GetInvert(this->fFullCovar); - return; -}; - - -void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::FillEventVariables(FitEvent *event) { - - if (event->NumFSParticle(13) == 0 || - event->NumFSParticle(211) == 0) - return; - - TLorentzVector Pnu = event->GetNeutrinoIn()->fP; - TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; - TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; - - double th_q3_pi = FitUtils::thq3pi_CC1pip_T2K(Pnu, Pmu, Ppip); - - fXVar = th_q3_pi; - - return; -}; - -//******************************************************************** -bool T2K_CC1pip_CH_XSec_1Dthq3pi_nu::isSignal(FitEvent *event) { -//******************************************************************** -// This sample requires pion directional information so can not include Michel tag sample -// i.e. will need to cut the pion phase space - return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); -} - diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h deleted file mode 100644 index 44e26b0..0000000 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef T2K_CC1PIP_CH_XSEC_1DTHQ3PI_NU_H_SEEN -#define T2K_CC1PIP_CH_XSEC_1DTHQ3PI_NU_H_SEEN - -#include "Measurement1D.h" - -class T2K_CC1pip_CH_XSec_1Dthq3pi_nu : public Measurement1D { -public: - T2K_CC1pip_CH_XSec_1Dthq3pi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); - virtual ~T2K_CC1pip_CH_XSec_1Dthq3pi_nu() {}; - - // Functions to deal with the input data and covariance - void SetDataValues(std::string fileLocation); - void SetCovarMatrix(std::string covarFile); - - void FillEventVariables(FitEvent *event); - bool isSignal(FitEvent *event); - - private: -}; - -#endif diff --git a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx new file mode 100644 index 0000000..9b442f7 --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx @@ -0,0 +1,173 @@ +#include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h" + +//******************************************************************** +T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(nuiskey samplekey) { +//******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu sample. \n" \ + "Target: CH \n" \ + "Flux: T2k Forward Horn Current nue + nuebar \n" \ + "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle(" "); + fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/(GeV/c)/nucleon)"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + + // CCQELike plot information + fSettings.SetTitle("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu"); + fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/PmuThetamu.root" ); + fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/PmuThetamu.root" ); + fSettings.DefineAllowedSpecies("numu"); + + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + //SetDataValues( fSettings.GetDataInput() ); + //SetCovarMatrix( fSettings.GetCovarInput() ); + SetHistograms(); + fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "Covariance_pmu_thetamu"); + covar = StatUtils::GetInvert(fFullCovar); + fDecomp = StatUtils::GetDecomp(fFullCovar); + SetShapeCovar(); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + +void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::SetHistograms() { + + TFile *data = new TFile(fSettings.GetDataInput().c_str(), "open"); + //std::string dataname = fSettings.Get + std::string dataname = "p_mu_theta_mu"; + + // Number of slices we have + const int nslices = 4; + int nbins = 0; + for (int i = 0; i < nslices; ++i) { + TH1D* slice = (TH1D*)data->Get(Form("%s_%i", dataname.c_str(), i)); + slice = (TH1D*)slice->Clone((fName+Form("_data_slice%i", i)).c_str()); + slice->Scale(1E-38); + slice->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str()); + slice->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); + fDataHist_Slices.push_back(slice); + fMCHist_Slices.push_back((TH1D*)slice->Clone((fName+Form("_mc_slice%i", i)).c_str())); + SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); + SetAutoProcessTH1(fMCHist_Slices[i]); + fMCHist_Slices[i]->Reset(); + fMCHist_Slices[i]->SetLineColor(kRed); + nbins += slice->GetXaxis()->GetNbins(); + } + + fDataHist = new TH1D(dataname.c_str(), dataname.c_str(), nbins, 0, nbins); + fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); + int bincount = 1; + for (int i = 0; i < nslices; ++i) { + for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { + fDataHist->SetBinContent(bincount, fDataHist_Slices[i]->GetBinContent(j+1)); + fDataHist->SetBinError(bincount, fDataHist_Slices[i]->GetBinError(j+1)); + TString title; + if (j == 0) { + title = "cos#theta_{#mu}="; + if (i == 0) { + title += "0-0.8, "; + } else if (i == 1) { + title += "0.8-0.85, "; + } else if (i == 2) { + title += "0.85-0.90, "; + } else if (i == 3) { + title += "0.90-1.00, "; + } + } + title += Form("p_{#mu}=%.2f-%.2f", fDataHist_Slices[i]->GetBinLowEdge(j+1), fDataHist_Slices[i]->GetBinLowEdge(j+2)); + fDataHist->GetXaxis()->SetBinLabel(bincount, title); + bincount++; + } + } + fDataHist->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str()); + fDataHist->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); +}; + + +void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillEventVariables(FitEvent *event) { + if (event->NumFSParticle(13) == 0) return; + + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + + double pmu = FitUtils::p(Pmu); + double costhmu = cos(FitUtils::th(Pnu, Pmu)); + + fXVar = pmu; + fYVar = costhmu; +}; + +void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillHistograms() { + Measurement1D::FillHistograms(); + if (Signal) { + FillMCSlice(fXVar, fYVar, Weight); + } +}; + +void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::ConvertEventRates(){ + + const int nslices = 4; + for (int i = 0; i < nslices; i++){ + fMCHist_Slices[i]->GetSumw2(); + } + + // Do standard conversion. + Measurement1D::ConvertEventRates(); + + // First scale MC slices also by their width in Y and Z + fMCHist_Slices[0]->Scale(1.0 / 0.80); + fMCHist_Slices[1]->Scale(1.0 / 0.05); + fMCHist_Slices[2]->Scale(1.0 / 0.05); + fMCHist_Slices[3]->Scale(1.0 / 0.10); + + // Now Convert into 1D list + fMCHist->Reset(); + int bincount = 1; + for (int i = 0; i < nslices; i++){ + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ + fMCHist->SetBinContent(bincount, fMCHist_Slices[i]->GetBinContent(j+1)); + bincount++; + } + } +} + +void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillMCSlice(double pmu, double cosmu, double weight) { + // Hard code the bin edges in here + if (cosmu < 0.8) { + fMCHist_Slices[0]->Fill(pmu, weight); + } else if (cosmu > 0.8 && cosmu < 0.85) { + fMCHist_Slices[1]->Fill(pmu, weight); + } else if (cosmu > 0.85 && cosmu < 0.90) { + fMCHist_Slices[2]->Fill(pmu, weight); + } else if (cosmu > 0.90 && cosmu < 1.00) { + fMCHist_Slices[3]->Fill(pmu, weight); + } +}; + +//******************************************************************** +bool T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::isSignal(FitEvent *event) { +//******************************************************************** + if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false; + + TLorentzVector Pnu = event->GetHMISParticle(14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + + double cos_th_mu = cos(FitUtils::th(Pnu,Pmu)); + if (cos_th_mu >= 0) return true; + return false; +}; + diff --git a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h new file mode 100644 index 0000000..3c8bbc5 --- /dev/null +++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h @@ -0,0 +1,30 @@ +#ifndef T2K_CC1PIP_CH_XSEC_2DPMUCOSMU_NU_H_SEEN +#define T2K_CC1PIP_CH_XSEC_2DPMUCOSMU_NU_H_SEEN + +#include "Measurement1D.h" + +class T2K_CC1pip_CH_XSec_2Dpmucosmu_nu : public Measurement1D { +public: + T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(nuiskey samplekey); + virtual ~T2K_CC1pip_CH_XSec_2Dpmucosmu_nu() {}; + + // Functions to deal with the input data and covariance + //void SetDataValues(std::string fileLocation); + //void SetCovarMatrix(std::string covarFile); + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + + void FillHistograms(); + + // Piggy-back off this one to fill the MC hist in slices + void ConvertEventRates(); + + private: + std::vector fMCHist_Slices; + std::vector fDataHist_Slices; + void FillMCSlice(double x, double y, double w); + void SetHistograms(); +}; + +#endif diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx index 76a92d8..b4be3af 100644 --- a/src/Utils/FitUtils.cxx +++ b/src/Utils/FitUtils.cxx @@ -1,1186 +1,1184 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "FitUtils.h" /* MISC Functions */ //******************************************************************** double *FitUtils::GetArrayFromMap(std::vector invals, std::map inmap) { //******************************************************************** double *outarr = new double[invals.size()]; int count = 0; for (size_t i = 0; i < invals.size(); i++) { outarr[count++] = inmap[invals.at(i)]; } return outarr; } /* MISC Event */ //******************************************************************** // Returns the kinetic energy of a particle in GeV double FitUtils::T(TLorentzVector part) { //******************************************************************** double E_part = part.E() / 1000.; double p_part = part.Vect().Mag() / 1000.; double m_part = sqrt(E_part * E_part - p_part * p_part); double KE_part = E_part - m_part; return KE_part; }; //******************************************************************** // Returns the momentum of a particle in GeV double FitUtils::p(TLorentzVector part) { //******************************************************************** double p_part = part.Vect().Mag() / 1000.; return p_part; }; double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); }; //******************************************************************** // Returns the angle between two particles in radians double FitUtils::th(TLorentzVector part1, TLorentzVector part2) { //******************************************************************** double th = part1.Vect().Angle(part2.Vect()); return th; }; double FitUtils::th(FitParticle *part1, FitParticle *part2) { return FitUtils::th(part1->fP, part2->fP); }; // T2K CC1pi+ helper functions // //******************************************************************** // Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH // paper // Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB //******************************************************************** double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { // Want this in GeV TVector3 p_mu = pmu.Vect() * (1. / 1000.); // Get the reconstructed Enu // We are not using Michel e sample, so we have pion kinematic information double Enu = EnuCC1piprec(pnu, pmu, ppi, true); // Get neutrino unit direction, multiply by reconstructed Enu TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu; TVector3 p_pi = ppi.Vect() * (1. / 1000.); // This is now in GeV TVector3 q3 = (p_nu - p_mu); // Want this in GeV double th_q3_pi = q3.Angle(p_pi); return th_q3_pi; } //******************************************************************** // Returns the q3 defined in Raquel's CC1pi+ on CH paper // Uses "MiniBooNE formula" for Enu //******************************************************************** double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { // Can't use the true Enu here; need to reconstruct it // We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without // pion kinematics // The bool false refers to that we don't have pion kinematics // Last bool refers to if we have pion kinematic information or not double Enu = EnuCC1piprec(pnu, pmu, ppi, false); TVector3 p_mu = pmu.Vect() * (1. / 1000.); TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu; double q3 = (p_nu - p_mu).Mag(); return q3; } //******************************************************************** // Returns the W reconstruction from Raquel CC1pi+ CH thesis // Uses the MiniBooNE formula Enu //******************************************************************** double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { double E_mu = pmu.E() / 1000.; double p_mu = pmu.Vect().Mag() / 1000.; double E_nu = EnuCC1piprec(pnu, pmu, ppi, false); double a1 = (E_nu + PhysConst::mass_neutron) - E_mu; double a2 = E_nu - p_mu; double wrec = sqrt(a1 * a1 - a2 * a2); return wrec; } //******************************************************** double FitUtils::ProtonQ2QErec(double pE, double binding) { //******************************************************** const double V = binding / 1000.; // binding potential const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass const double mn_eff = mn - V; // effective proton mass const double pki = (pE / 1000.0) - mp; double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff); return q2qe; }; //******************************************************************** double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV const double V = binding / 1000.; // binding potential const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double mN_eff = mn - V; double mN_oth = mp; if (!neutrino) { mN_eff = mp - V; mN_oth = mn; } double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double ml = sqrt(el * el - pl * pl); // lepton mass double rEnu = (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / (2 * (mN_eff - el + pl * costh)); return rEnu; }; // Another good old helper function double FitUtils::EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding, bool neutrino) { return EnuQErec(pmu, cos(pnu.Vect().Angle(pmu.Vect())), binding, neutrino); } double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) { double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double ml = sqrt(el * el - pl * pl); // lepton mass double rEnu = EnuQErec(pmu, costh, binding, neutrino); double q2 = -ml * ml + 2. * rEnu * (el - pl * costh); return q2; }; double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) { double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino); return q2qe; } double FitUtils::EnuQErec(double pl, double costh, double binding, bool neutrino) { if (pl < 0) return 0.; // Make sure nobody is silly double mN_eff = PhysConst::mass_neutron - binding / 1000.; double mN_oth = PhysConst::mass_proton; if (!neutrino) { mN_eff = PhysConst::mass_proton - binding / 1000.; mN_oth = PhysConst::mass_neutron; } double ml = PhysConst::mass_muon; double el = sqrt(pl * pl + ml * ml); double rEnu = (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / (2 * (mN_eff - el + pl * costh)); return rEnu; }; double FitUtils::Q2QErec(double pl, double costh, double binding, bool neutrino) { if (pl < 0) return 0.; // Make sure nobody is silly double ml = PhysConst::mass_muon; double el = sqrt(pl * pl + ml * ml); double rEnu = EnuQErec(pl, costh, binding, neutrino); double q2 = -ml * ml + 2. * rEnu * (el - pl * costh); return q2; }; //******************************************************************** // Reconstructs Enu for CC1pi0 // Very similar for CC1pi+ reconstruction double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0) { //******************************************************************** double E_mu = pmu.E() / 1000; double p_mu = pmu.Vect().Mag() / 1000; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double E_pi0 = ppi0.E() / 1000; double p_pi0 = ppi0.Vect().Mag() / 1000; double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0); double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect()); const double m_n = PhysConst::mass_neutron; // neutron mass const double m_p = PhysConst::mass_proton; // proton mass double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect()); double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p - 2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu - 2 * p_pi0 * p_mu * cos(th_pi0_mu)) / (2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) - p_mu * cos(th_nu_mu) - m_n)); return rEnu; }; //******************************************************************** // Reconstruct Q2 for CC1pi0 // Beware: uses true Enu, not reconstructed Enu double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0) { //******************************************************************** double E_mu = pmu.E() / 1000.; // energy of lepton in GeV double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy // Use true neutrino energy double rEnu = pnu.E() / 1000.; double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu)); return q2; }; //******************************************************************** // Reconstruct Enu for CC1pi+ // pionInfo reflects if we're using pion kinematics or not // In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info // is lost and Enu is reconstructed without pion kinematics double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi, bool pionInfo) { //******************************************************************** double E_mu = pmu.E() / 1000.; double p_mu = pmu.Vect().Mag() / 1000.; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double E_pi = ppi.E() / 1000.; double p_pi = ppi.Vect().Mag() / 1000.; double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi); const double m_n = PhysConst::mass_neutron; // neutron/proton mass // should really take proton mass for proton interaction, neutron for neutron // interaction. However, difference is pretty much negligable here! // need this angle too double th_nu_pi = pnu.Vect().Angle(ppi.Vect()); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double th_pi_mu = ppi.Vect().Angle(pmu.Vect()); double rEnu = -999; // If experiment doesn't have pion kinematic information (T2K CC1pi+ CH // example of this; Michel e sample has no directional information on pion) // We'll need to modify the reconstruction variables if (!pionInfo) { // Assumes that pion angle contribution contributes net zero // Also assumes the momentum of reconstructed pions when using Michel // electrons is 130 MeV // So need to redefine E_pi and p_pi // It's a little unclear what happens to the pmu . ppi term (4-vector); do // we include the angular contribution? // This below doesn't th_nu_pi = M_PI / 2.; p_pi = 0.130; E_pi = sqrt(p_pi * p_pi + m_pi * m_pi); // rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi - // 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n)); } rEnu = (m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu - 2 * p_pi * p_mu * cos(th_pi_mu)) / (2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n)); return rEnu; }; //******************************************************************** // Reconstruct neutrino energy from outgoing particles; will differ from the // actual neutrino energy. Here we use assumption of a Delta resonance double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) { //******************************************************************** const double m_Delta = PhysConst::mass_delta; // PDG value for Delta mass in GeV const double m_n = PhysConst::mass_neutron; // neutron/proton mass // should really take proton mass for proton interaction, neutron for neutron // interaction. However, difference is pretty much negligable here! double E_mu = pmu.E() / 1000; double p_mu = pmu.Vect().Mag() / 1000; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) / (2 * (m_n - E_mu + p_mu * cos(th_nu_mu))); return rEnu; }; // MOVE TO T2K UTILS! //******************************************************************** // Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN // // Supposedly includes pion direction and binding energy of target nucleon // I'm not convinced (yet), maybe double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { //******************************************************************** // Unit vector for neutrino momentum TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E()); double E_mu = pmu.E() / 1000.; TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.); double E_pi = ppi.E() / 1000.; TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.); - double E_bind = - 27. / 1000.; // This should be roughly correct for CH; but not clear! + double E_bind = 25. / 1000.; double m_p = PhysConst::mass_proton; // Makes life a little easier, gonna square this one double a1 = m_p - E_bind - E_mu - E_pi; // Just gets the magnitude square of the muon and pion momentum vectors double a2 = (p_mu_vect + p_pi_vect).Mag2(); // Gets the somewhat complicated scalar product between neutrino and // (p_mu+p_pi), scaled to Enu double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect); double rEnu = (m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3)); return rEnu; } //******************************************************************** // Reconstructed Q2 for CC1pi+ // // enuType describes how the neutrino energy is reconstructed // 0 uses true neutrino energy; case for MINERvA and MiniBooNE // 1 uses "extended MiniBooNE" reconstructed (T2K CH) // 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH) // "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the // case for T2K when using Michel tag) (T2K CH) // 3 uses Delta for reconstruction (T2K CH) double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi, int enuType, bool pionInfo) { //******************************************************************** double E_mu = pmu.E() / 1000.; // energy of lepton in GeV double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // Different ways of reconstructing the neutrino energy; default is just to // use the truth (case 0) double rEnu = -999; switch (enuType) { case 0: // True neutrino energy, default; this is the case for when Q2 // defined as Q2 true (MiniBooNE, MINERvA) rEnu = pnu.E() / 1000.; break; case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+ // CH T2K analysis - // Definitely uses pion info :) rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi); break; case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's // CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O // Can have this with and without pion info, depending on if Michel tag // used (Raquel) rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo); break; case 3: rEnu = EnuCC1piprecDelta(pnu, pmu); break; } // No need for default here since default value of enuType = 0, defined in // header double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu)); return q2; }; //******************************************************************** // Returns the reconstructed W from a nucleon and an outgoing pion // // Could do this a lot more clever (pp + ppi).Mag() would do the job, but this // would be less instructive //******************************************************************** double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) { double E_p = pp.E(); double p_p = pp.Vect().Mag(); double m_p = sqrt(E_p * E_p - p_p * p_p); double E_pi = ppi.E(); double p_pi = ppi.Vect().Mag(); double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi); double th_p_pi = pp.Vect().Angle(ppi.Vect()); // fairly easy thing to derive since bubble chambers measure the proton! double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi - 2 * p_pi * p_p * cos(th_p_pi)); return invMass; }; //******************************************************** // Reconstruct the hadronic mass using neutrino and muon // Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be // reconstructed Enu; so gives reconstructed Wrec which most of the time isn't // used! // // Only MINERvA uses this so far; and the Enu is Enu_true // If we want W_true need to take initial state nucleon motion into account // Return value is in MeV!!! double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) { //******************************************************** double E_mu = pmu.E(); double p_mu = pmu.Vect().Mag(); double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // The factor of 1000 is necessary for downstream functions const double m_p = PhysConst::mass_proton * 1000; // MINERvA cut on W_exp which is tuned to W_true; so use true Enu from // generators double E_nu = pnu.E(); double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + 2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu))); return w_rec; }; //******************************************************** // Reconstruct the true hadronic mass using the initial state and muon // Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be // reconstructed Enu; so gives reconstructed Wrec which most of the time isn't // used! // // No one seems to use this because it's fairly MC dependent! // Return value is in MeV!!! double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc) { //******************************************************** // Could simply do the TLorentzVector operators here but this is more // instructive? // ... and prone to errors ... double E_mu = pmu.E(); double p_mu = pmu.Vect().Mag(); double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double E_nuc = pnuc.E(); double p_nuc = pnuc.Vect().Mag(); double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc); double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect()); double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect()); double E_nu = pnu.E(); double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu + 2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu + 2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc - 2 * E_nu * p_nuc * cos(th_nu_nuc)); return w_rec; }; double FitUtils::SumKE_PartVect(std::vector const fps) { double sum = 0.0; for (size_t p_it = 0; p_it < fps.size(); ++p_it) { sum += fps[p_it]->KE(); } return sum; } double FitUtils::SumTE_PartVect(std::vector const fps) { double sum = 0.0; for (size_t p_it = 0; p_it < fps.size(); ++p_it) { sum += fps[p_it]->E(); } return sum; } /* E Recoil */ double FitUtils::GetErecoil_TRUE(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1) continue; // Add Up KE of protons and TE of everything else if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) { Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); } else { Erecoil += event->PartInfo(i)->fP.E(); } } return Erecoil; } double FitUtils::GetErecoil_CHARGED(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1) continue; // Skip Neutral particles if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 || event->PartInfo(i)->fPID == 22) continue; // Add Up KE of protons and TE of everything else if (event->PartInfo(i)->fPID == 2212) { Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); } else { Erecoil += event->PartInfo(i)->fP.E(); } } return Erecoil; } // MOVE TO MINERVA Utils! double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == 13) continue; // Skip Neutrons particles if (event->PartInfo(i)->fPID == 2112) continue; int PID = event->PartInfo(i)->fPID; // KE of Protons and charged pions if (PID == 2212 or PID == 211 or PID == -211) { // Erecoil += FitUtils::T(event->PartInfo(i)->fP); Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); // Total Energy of non-neutrons // } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) != // 14) { } else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) { Erecoil += (event->PartInfo(i)->fP.E()); } } return Erecoil; } // MOVE TO MINERVA Utils! // The alternative Eavailble definition takes true q0 and subtracts the kinetic energy of neutrons and pion masses // returns in MeV double FitUtils::Eavailable(FitEvent *event) { double Eav = 0.0; // Now take q0 and subtract Eav double q0 = event->GetNeutrinoIn()->fP.E(); // Get the pdg of incoming neutrino int ISPDG = event->GetBeamNeutrinoPDG(); // For CC if (event->IsCC()) q0 -= event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.E(); else q0 -= event->GetHMFSParticle(ISPDG)->fP.E(); for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; int PID = event->PartInfo(i)->fPID; // Neutrons if (PID == 2112) { // Adding kinetic energy of neutron Eav += FitUtils::T(event->PartInfo(i)->fP)*1000.; // All pion masses } else if (abs(PID) == 211 || PID == 111) { Eav += event->PartInfo(i)->fP.M(); } } return q0-Eav; } TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { TVector3 pnUnit = planarNormal.Unit(); double inpProjectPN = inp.Dot(pnUnit); TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit); return InPlanarInput; } TVector3 GetUnitVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { return GetVectorInTPlane(inp, planarNormal).Unit(); } Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false) { TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal); return PiMinus ? acos(V_lepton_T.Dot(V_other_T)) : (M_PI - acos(V_lepton_T.Dot(V_other_T))); } TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal) { TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetVectorInTPlane(V_other, Normal); return V_lepton_T + V_other_T; } Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false) { TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal); return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus); } double FitUtils::Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 TLorentzVector Pnu = event->GetNeutrinoIn()->fP; int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HighestMomentum = Protons[i]->p(); HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; // Get highest momentum proton in allowed proton range TVector3 HadronP = Pprot.Vect(); // If we don't have a CC0pi signal definition we also add in pion momentum if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } // Count up pion momentum TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaPT(LeptonP, HadronP, NuP).Mag(); } double FitUtils::Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 TLorentzVector Pnu = event->GetNeutrinoIn()->fP; int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HighestMomentum = Protons[i]->p(); HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; // Get highest momentum proton in allowed proton range TVector3 HadronP = Pprot.Vect(); if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaPhiT(LeptonP, HadronP, NuP); } double FitUtils::Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 TLorentzVector Pnu = event->GetNeutrinoIn()->fP; int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HighestMomentum = Protons[i]->p(); HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; // Get highest momentum proton in allowed proton range TVector3 HadronP = Pprot.Vect(); if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaAlphaT(LeptonP, HadronP, NuP); } // As defined in PhysRevC.95.065501 // Using prescription from arXiv 1805.05486 // Returns in GeV double FitUtils::Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) { const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 TLorentzVector Pnu = event->GetNeutrinoIn()->fP; int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { // Update the highest momentum particle if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HighestMomentum = Protons[i]->p(); HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; // Get highest momentum proton in allowed proton range TVector3 HadronP = Pprot.Vect(); //TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.; double const eh = Pprot.E()/1000.; if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP); double dptMag = dpt.Mag()/1000.; double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501) double map = ma - mn + 0.02713; // remnant mass double pmul = LeptonP.Dot(NuP.Unit())/1000.; double phl = HadronP.Dot(NuP.Unit())/1000.; //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.; //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.; double R = ma + pmul + phl - el - eh; double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R); //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl)); //std::cout << "Diagnostics: " << std::endl; //std::cout << "mn: " << mn << std::endl; //std::cout << "ma: " << ma << std::endl; //std::cout << "map: " << map << std::endl; //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl; //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl; //std::cout << "pmul: " << pmul << std::endl; //std::cout << "phl: " << phl << std::endl; //std::cout << "el: " << el << std::endl; //std::cout << "eh: " << eh << std::endl; //std::cout << "R: " << R << std::endl; //std::cout << "dptMag: " << dptMag << std::endl; //std::cout << "dpl: " << dpl << std::endl; //std::cout << "pn_reco: " << pn_reco << std::endl; return pn_reco; } double FitUtils::Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi) { const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 TLorentzVector Pnu = event->GetNeutrinoIn()->fP; int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { // Update the highest momentum particle if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HighestMomentum = Protons[i]->p(); HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; // Get highest momentum proton in allowed proton range TVector3 HadronP = Pprot.Vect(); //TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.; double const eh = Pprot.E()/1000.; if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP); double dptMag = dpt.Mag()/1000.; //double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501) //double map = ma - mn + 0.02713; // remnant mass double ma = 6*mn + 6*mp - 0.34381; // target mass (E is from PhysRevC.95.065501) double map = ma - mn + 0.02713; // remnant mass double pmul = LeptonP.Dot(NuP.Unit())/1000.; double phl = HadronP.Dot(NuP.Unit())/1000.; //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.; //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.; double R = ma + pmul + phl - el - eh; double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R); //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl)); //std::cout << "Diagnostics: " << std::endl; //std::cout << "mn: " << mn << std::endl; //std::cout << "ma: " << ma << std::endl; //std::cout << "map: " << map << std::endl; //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl; //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl; //std::cout << "pmul: " << pmul << std::endl; //std::cout << "phl: " << phl << std::endl; //std::cout << "el: " << el << std::endl; //std::cout << "eh: " << eh << std::endl; //std::cout << "R: " << R << std::endl; //std::cout << "dptMag: " << dptMag << std::endl; //std::cout << "dpl: " << dpl << std::endl; //std::cout << "pn_reco: " << pn_reco << std::endl; return pn_reco; } // Get Cos theta with Adler angles double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) { // Get the "resonance" lorentz vector (pion proton system) TLorentzVector Pres = Pprot + Ppi; // Boost the particles into the resonance rest frame so we can define the x,y,z axis Pnu.Boost(Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppi.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag())); // Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle double costhAdler = cos(Ppi.Vect().Angle(zAxis)); return costhAdler; } // Get phi with Adler angles, a bit more complicated... double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) { // Get the "resonance" lorentz vector (pion proton system) TLorentzVector Pres = Pprot + Ppi; // Boost the particles into the resonance rest frame so we can define the x,y,z axis Pnu.Boost(Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppi.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag())); // The y-axis is then defined perpendicular to z and muon momentum in the resonance frame TVector3 yAxis = Pnu.Vect().Cross(Pmu.Vect()); yAxis *= 1.0/double(yAxis.Mag()); // And the x-axis is then simply perpendicular to z and x TVector3 xAxis = yAxis.Cross(zAxis); xAxis *= 1.0/double(xAxis.Mag()); double x = Ppi.Vect().Dot(xAxis); double y = Ppi.Vect().Dot(yAxis); //double z = Ppi.Vect().Dot(zAxis); double newphi = atan2(y, x)*(180./M_PI); // Convert negative angles to positive if (newphi < 0.0) newphi += 360.0; // Old silly method before atan2 /* // Then finally construct phi as the angle between pion projection and x axis // Get the project of the pion momentum on to the zaxis TVector3 PiVectZ = zAxis*Ppi.Vect().Dot(zAxis); // The subtract the projection off the pion vector to get to get the plane TVector3 PiPlane = Ppi.Vect() - PiVectZ; double phi = -999.99; if (PiPlane.Y() > 0) { phi = (180./M_PI)*PiPlane.Angle(xAxis); } else if (PiPlane.Y() < 0) { phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis)); } else if (PiPlane.Y() == 0) { TRandom3 rand; double randNo = rand.Rndm(); if (randNo > 0.5) { phi = (180./M_PI)*PiPlane.Angle(xAxis); } else { phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis)); } } */ return newphi; } //******************************************************************** double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double el = pmu.E() / 1000.; //double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double enu = EnuQErec(pmu, costh, binding, neutrino); double ep_inf = enu - el + mp; double pp_inf = sqrt(ep_inf * ep_inf - mp * mp); return pp_inf; }; //******************************************************************** TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass //const double mp = PhysConst::mass_proton; // proton mass double pl_x = pmu.X() / 1000.; double pl_y = pmu.Y() / 1000.; double pl_z= pmu.Z() / 1000.; double enu = EnuQErec(pmu, costh, binding, neutrino); TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu); return tpp_inf; }; //******************************************************************** double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double enu = EnuQErec(pmu, costh, binding, neutrino); double ep_inf = enu - el + mp; double pp_inf = sqrt(ep_inf * ep_inf - mp * mp); double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf); return cth_inf; };