Hello everybody !
I’m trying to fit a multiple gaussian (7 peaks) with a 3rd grade polynomial background,
I’m sure for the polynomial’s parameters.
For every gaussian I used TSpectrum to set Mean and Constant, and I used FWHM in another sheet to get sigma.
Anybody knows why it doesn’t work ?
Thank you in advance!!
//
// MGF_Background_3_0.cpp
//
//
// Created by Simone Azeglio on 01/03/17.
//
//
// Third Grade Polynomial background function
Double_t background(Double_t *x, Double_t *par) {
return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]+ par[3]*x[0]*x[0]*x[0];
}
// Multi Gaussian function
Double_t MultiGaussianPeaks(Double_t *x, Double_t *par) {
return (par[4]*TMath::Exp(-0.5*((x[0]-par[5])*(x[0]-par[5]))/(par[6]*par[6]))+
par[7]*TMath::Exp(-0.5*((x[0]-par[8])*(x[0]-par[8]))/(par[9]*par[9]))+
par[10]*TMath::Exp(-0.5*((x[0]-par[11])*(x[0]-par[11]))/(par[12]*par[12]))+
par[13]*TMath::Exp(-0.5*((x[0]-par[14])*(x[0]-par[14]))/(par[15]*par[15]))+
par[16]*TMath::Exp(-0.5*((x[0]-par[17])*(x[0]-par[17]))/(par[18]*par[18]))+
par[19]*TMath::Exp(-0.5*((x[0]-par[20])*(x[0]-par[20]))/(par[21]*par[21]))+
par[22]*TMath::Exp(-0.5*((x[0]-par[23])*(x[0]-par[23]))/(par[24]*par[24]))) ;
}
// Sum of background and peak function
Double_t fitFunction(Double_t *x, Double_t *par) {
return background(x,&par[0]) + MultiGaussianPeaks(x,&par[4]);
}
void MGF_Background_3_0() {
// Channels
const Int_t np = 2000 ;
// Counts
Double_t x[np] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 2, 0, 2, 3, 1, 0, 0, 1, 2, 1, 4, 1, 4, 10, 2, 6, 10, 9, 13, 25, 29, 35, 68, 61, 116, 178, 276, 441, 629, 870, 1071, 1319, 1425, 1672, 1563, 1513, 1357, 1178, 902, 680, 552, 376, 228, 165, 97, 81, 56, 32, 47, 45, 38, 35, 30, 34, 35, 35, 32, 30, 23, 24, 25, 26, 31, 33, 24, 28, 22, 38, 33, 38, 31, 25, 36, 38, 50, 58, 62, 54, 72, 92, 114, 162, 197, 273, 383, 506, 703, 952, 1263, 1678, 2155, 2491, 2979, 3286, 3502, 3777, 3709, 3480, 3122, 2730, 2308, 1897, 1435, 1117, 835, 611, 476, 308, 310, 233, 229, 224, 205, 212, 252, 224, 218, 238, 239, 222, 245, 232, 256, 221, 263, 237, 263, 260, 246, 232, 252, 248, 262, 296, 279, 290, 326, 326, 359, 375, 435, 470, 507, 654, 782, 950, 1197, 1562, 1947, 2373, 2891, 3458, 3885, 4317, 4790, 4980, 5309, 5214, 5111, 4665, 4220, 3745, 3166, 2671, 2129, 1683, 1332, 1137, 878, 756, 645, 599, 609, 558, 547, 547, 563, 576, 662, 643, 596, 642, 650, 668, 649, 652, 694, 665, 715, 758, 739, 703, 783, 769, 759, 791, 886, 928, 874, 975, 987, 1055, 1141, 1272, 1512, 1587, 1817, 2062, 2595, 2957, 3529, 3823, 4430, 4908, 5363, 5643, 5819, 5987, 5908, 5611, 5452, 5100, 4561, 3843, 3315, 2892, 2458, 2090, 1752, 1513, 1361, 1266, 1117, 1086, 1124, 1053, 1107, 1148, 1073, 1150, 1184, 1110, 1192, 1195, 1254, 1218, 1283, 1229, 1252, 1316, 1359, 1333, 1384, 1397, 1420, 1512, 1481, 1574, 1633, 1537, 1673, 1763, 1901, 2015, 2110, 2403, 2576, 2920, 3099, 3380, 3853, 4037, 4509, 4985, 5252, 5571, 5764, 5790, 5942, 5860, 5621, 5235, 5020, 4704, 4351, 3691, 3423, 2831, 2621, 2257, 2090, 1900, 1695, 1645, 1682, 1611, 1682, 1592, 1670, 1664, 1729, 1683, 1796, 1694, 1877, 1704, 1827, 1971, 1958, 1936, 1939, 1922, 1979, 1993, 2064, 2085, 2168, 2114, 2205, 2307, 2266, 2438, 2489, 2577, 2722, 2898, 3103, 3164, 3438, 3607, 3969, 4128, 4412, 4770, 5168, 4987, 5244, 5470, 5389, 5455, 5362, 5245, 4868, 4666, 4393, 3933, 3769, 3406, 3229, 2845, 2621, 2520, 2415, 2271, 2300, 2179, 2061, 2155, 2092, 2156, 2166, 2184, 2201, 2243, 2272, 2278, 2433, 2305, 2449, 2289, 2399, 2441, 2520, 2545, 2627, 2564, 2523, 2618, 2623, 2754, 2756, 2781, 2977, 2905, 3047, 3115, 3402, 3487, 3769, 3740, 3878, 4043, 4169, 4387, 4406, 4755, 4800, 4843, 4729, 4924, 4772, 4770, 4531, 4275, 4124, 3788, 3759, 3530, 3337, 3114, 3048, 2774, 2779, 2650, 2615, 2508, 2392, 2449, 2407, 2581, 2420, 2556, 2572, 2561, 2577, 2640, 2561, 2617, 2678, 2703, 2666, 2673, 2785, 2754, 2859, 2878, 2781, 2908, 2962, 3005, 3059, 3060, 3059, 3100, 3273, 3342, 3409, 3469, 3530, 3560, 3831, 3655, 3882, 3977, 4112, 4104, 4259, 4287, 4256, 4121, 4199, 4024, 3994, 3883, 3830, 3510, 3556, 3350, 3171, 3229, 2947, 2864, 2869, 2763, 2709, 2620, 2613, 2661, 2771, 2575, 2666, 2743, 2687, 2676, 2694, 2808, 2762, 2731, 2742, 2765, 2775, 2955, 2816, 2877, 2915, 2982, 2894, 2814, 3067, 3014, 3042, 2999, 3117, 3041, 3120, 3128, 3247, 3301, 3263, 3299, 3412, 3427, 3579, 3549, 3515, 3622, 3585, 3675, 3651, 3622, 3533, 3546, 3514, 3340, 3325, 3208, 3056, 3103, 3061, 2899, 2877, 2781, 2803, 2702, 2724, 2597, 2638, 2640, 2588, 2509, 2541, 2581, 2652, 2586, 2538, 2603, 2594, 2574, 2661, 2652, 2728, 2716, 2729, 2688, 2724, 2711, 2723, 2806, 2771, 2793, 2863, 2820, 2814, 2894, 2890, 2917, 3038, 2955, 3044, 3036, 3046, 3115, 3097, 3041, 3055, 3067, 3144, 3103, 3127, 3043, 3036, 3023, 2988, 2915, 2916, 2852, 2706, 2810, 2794, 2650, 2602, 2602, 2557, 2485, 2556, 2502, 2485, 2511, 2430, 2492, 2453, 2327, 2418, 2365, 2482, 2463, 2569, 2522, 2548, 2341, 2489, 2510, 2446, 2556, 2437, 2497, 2426, 2519, 2426, 2530, 2517, 2617, 2479, 2568, 2542, 2535, 2604, 2617, 2640, 2539, 2710, 2612, 2669, 2615, 2638, 2651, 2512, 2690, 2666, 2540, 2587, 2658, 2531, 2494, 2502, 2493, 2489, 2432, 2342, 2375, 2333, 2324, 2198, 2189, 2153, 2261, 2164, 2103, 2225, 2168, 2171, 2157, 2197, 2167, 2141, 2256, 2190, 2153, 2220, 2275, 2214, 2149, 2105, 2230, 2244, 2223, 2213, 2157, 2122, 2234, 2126, 2147, 2212, 2254, 2105, 2269, 2206, 2255, 2183, 2195, 2238, 2311, 2250, 2155, 2217, 2181, 2062, 2104, 2128, 2083, 2091, 2223, 2099, 2080, 2060, 2088, 2028, 2031, 1981, 1932, 2091, 1938, 1895, 1834, 1936, 1885, 1862, 1947, 1949, 1839, 1909, 1897, 1831, 1845, 1909, 1815, 1884, 1891, 1816, 1907, 1833, 1880, 1865, 1784, 1854, 1896, 1873, 1907, 1840, 1910, 1863, 1855, 1827, 1889, 1803, 1798, 1822, 1878, 1758, 1859, 1815, 1833, 1825, 1734, 1843, 1728, 1740, 1799, 1814, 1766, 1734, 1723, 1700, 1816, 1751, 1744, 1626, 1643, 1604, 1641, 1685, 1649, 1534, 1625, 1630, 1559, 1624, 1631, 1634, 1611, 1570, 1577, 1538, 1488, 1647, 1557, 1629, 1558, 1601, 1551, 1603, 1504, 1587, 1516, 1522, 1526, 1552, 1464, 1483, 1593, 1562, 1477, 1519, 1520, 1523, 1447, 1462, 1518, 1475, 1514, 1443, 1432, 1561, 1398, 1461, 1474, 1543, 1493, 1377, 1492, 1389, 1445, 1459, 1355, 1420, 1446, 1385, 1351, 1328, 1264, 1275, 1384, 1371, 1316, 1383, 1308, 1281, 1284, 1281, 1262, 1340, 1220, 1266, 1290, 1254, 1275, 1272, 1254, 1281, 1301, 1256, 1197, 1299, 1261, 1238, 1204, 1317, 1204, 1223, 1268, 1214, 1211, 1200, 1278, 1217, 1270, 1194, 1142, 1214, 1153, 1201, 1137, 1212, 1202, 1148, 1164, 1138, 1152, 1147, 1122, 1143, 1110, 1122, 1163, 1069, 1124, 1122, 1114, 1062, 1091, 1052, 1094, 1042, 1113, 1032, 1058, 1006, 1066, 1041, 1016, 1019, 1029, 997, 980, 1005, 1039, 1017, 1015, 964, 1025, 997, 981, 992, 1029, 980, 1002, 1003, 981, 997, 973, 927, 970, 926, 996, 974, 957, 885, 917, 999, 965, 967, 956, 964, 863, 944, 895, 901, 948, 886, 877, 883, 891, 907, 840, 912, 812, 892, 822, 877, 854, 790, 802, 805, 844, 833, 792, 820, 787, 830, 827, 835, 805, 810, 801, 764, 770, 788, 828, 732, 832, 808, 829, 779, 780, 717, 774, 753, 770, 764, 765, 774, 759, 716, 720, 781, 710, 715, 707, 752, 744, 763, 773, 734, 702, 721, 699, 735, 738, 696, 702, 725, 657, 717, 708, 690, 701, 683, 681, 659, 635, 641, 690, 676, 618, 636, 619, 614, 597, 639, 617, 624, 585, 605, 623, 614, 579, 605, 601, 599, 579, 566, 639, 640, 587, 561, 556, 552, 598, 561, 568, 579, 564, 576, 532, 487, 574, 560, 541, 567, 586, 530, 576, 589, 532, 557, 496, 519, 517, 506, 524, 487, 466, 476, 540, 499, 524, 551, 507, 492, 498, 484, 508, 522, 489, 453, 493, 475, 508, 484, 496, 480, 476, 486, 456, 492, 433, 501, 452, 450, 446, 449, 470, 439, 481, 424, 428, 455, 450, 456, 431, 406, 428, 379, 428, 435, 413, 426, 442, 399, 416, 420, 420, 424, 393, 390, 378, 421, 395, 380, 405, 401, 413, 424, 374, 387, 357, 412, 382, 365, 381, 372, 362, 356, 395, 362, 363, 383, 322, 349, 360, 371, 357, 342, 354, 351, 333, 346, 348, 344, 316, 343, 322, 327, 316, 335, 336, 367, 325, 329, 352, 314, 305, 311, 298, 297, 314, 349, 330, 305, 261, 315, 292, 300, 319, 282, 285, 296, 265, 281, 292, 260, 272, 276, 273, 263, 298, 268, 294, 264, 294, 289, 253, 260, 265, 266, 255, 260, 272, 224, 289, 252, 263, 240, 276, 245, 269, 223, 230, 236, 243, 242, 259, 217, 273, 218, 224, 254, 262, 254, 242, 194, 236, 241, 238, 218, 216, 257, 219, 207, 207, 225, 217, 229, 197, 224, 223, 248, 232, 199, 230, 177, 241, 208, 192, 234, 213, 203, 219, 200, 220, 222, 189, 175, 183, 188, 217, 199, 201, 194, 177, 204, 202, 176, 215, 186, 209, 185, 196, 192, 189, 185, 181, 162, 194, 191, 175, 163, 164, 162, 185, 144, 155, 138, 163, 166, 153, 147, 163, 154, 159, 155, 182, 153, 169, 154, 172, 170, 141, 149, 164, 144, 145, 165, 141, 139, 159, 127, 154, 140, 151, 146, 150, 133, 142, 162, 130, 144, 128, 139, 142, 129, 126, 129, 128, 116, 110, 128, 138, 132, 135, 130, 126, 151, 129, 107, 130, 130, 126, 125, 104, 122, 131, 111, 122, 118, 133, 127, 121, 116, 106, 127, 114, 119, 106, 132, 114, 103, 117, 109, 107, 105, 95, 113, 90, 110, 105, 105, 100, 106, 105, 117, 101, 118, 92, 101, 109, 106, 114, 106, 105, 94, 88, 83, 87, 89, 89, 78, 94, 96, 84, 79, 88, 87, 88, 100, 90, 98, 105, 90, 88, 102, 91, 79, 72, 77, 83, 77, 80, 89, 81, 73, 92, 80, 81, 81, 78, 81, 90, 65, 67, 82, 84, 72, 72, 67, 64, 78, 80, 93, 82, 56, 82, 81, 74, 72, 81, 62, 73, 65, 68, 78, 66, 45, 73, 74, 63, 53, 71, 76, 64, 65, 58, 72, 58, 55, 65, 63, 60, 86, 67, 63, 62, 72, 58, 51, 71, 61, 72, 51, 62, 75, 57, 42, 52, 52, 62, 57, 62, 66, 46, 52, 62, 63, 55, 51, 62, 62, 39, 46, 45, 65, 49, 37, 64, 65, 43, 45, 47, 35, 48, 56, 60, 45, 55, 51, 48, 43, 43, 39, 62, 48, 48, 37, 36, 41, 50, 45, 44, 50, 39, 44, 36, 44, 28, 44, 51, 38, 36, 34, 47, 40, 44, 37, 40, 30, 42, 45, 40, 41, 34, 38, 46, 38, 39, 35, 32, 36, 37, 32, 44, 47, 31, 40, 35, 30, 35, 34, 33, 23, 38, 22, 31, 29, 30, 39, 34, 32, 33, 36, 27, 39, 35, 24, 32, 30, 41, 34, 27, 31, 23, 30, 25, 31, 26, 22, 32, 29, 23, 21, 31, 31, 27, 31, 38, 26, 23, 27, 29, 26, 22, 24, 24, 26, 28, 24, 28, 24, 32, 27, 23, 21, 27, 23, 29, 28, 31, 23, 21, 16, 18, 21, 26, 23, 14, 20, 26, 16, 24, 19, 18, 15, 26, 17, 22, 21, 26, 27, 27, 20, 20, 14, 25, 16, 17, 13, 31, 18, 21, 16, 21, 12, 20, 13, 14, 24, 22, 22, 17, 20, 27, 23, 23, 11, 8, 20, 28, 15, 20, 13, 22, 15, 14, 10, 14, 19, 17, 17, 13, 22, 13, 15, 15, 19, 20, 16, 7, 21, 19, 20, 8, 20, 8, 12, 11, 11, 22, 18, 16, 17, 12, 14, 13, 11, 15, 11, 17, 13, 11, 19, 9, 14, 12, 8, 17, 14, 14, 17, 14, 12, 10, 14, 16, 15, 9, 10, 10, 12, 12, 8, 18, 16, 14, 8, 12, 5, 13, 13, 13, 15, 5, 10, 9, 7, 20, 6, 5, 11, 6, 4, 15, 5, 6, 9, 8, 10, 6, 7, 9, 8, 10, 13, 5, 10, 11, 7, 6, 8, 8, 7, 7, 12, 10, 5, 14, 5, 6, 7, 10, 7, 6, 9, 6, 10, 8, 4, 14, 6, 6, 14, 7, 12, 6, 11, 11, 7, 8, 6, 12, 7, 7, 5, 5, 8, 6, 3, 4, 4, 8, 9, 8, 8, 6, 7, 8, 5, 8, 5, 5, 8, 4, 6, 9, 6, 5, 4, 3, 5, 6, 6, 7, 4, 5, 11, 5, 6, 5, 6, 4, 3, 3, 9, 3, 3, 9, 4, 5, 8, 5, 6, 1, 4, 4, 8, 4, 5, 2, 4, 4, 3, 11, 6, 8, 4, 4, 7, 7, 7, 3, 0, 4, 7, 4, 8, 1, 5, 3, 5, 1, 1, 1, 4, 4, 0, 6, 9, 5, 2, 1, 1, 2, 5, 5, 2, 2, 3, 3, 5, 4, 1, 1, 2, 4, 5, 4};
Double_t xmin = -996.5;
Double_t xmax = 14957.0;
TH1F *histo = new TH1F("MGF","Multi Gaussian Fit on Background", np, xmin, xmax);
histo->SetMaximum(7000);
histo->GetXaxis()->SetTitle("Carica(Canali-ADC)");
histo->GetYaxis()->SetTitle("Conteggi");
for(int i=0; i < np; i++) {
// we use these methods to explicitly set the content
// and error instead of using the fill method.
histo->SetBinContent(i+1,x[i]);
histo->SetBinError(i+1,TMath::Sqrt(x[i]));
}
//----- USING TSPECTRUM
// Searching Peaks
TSpectrum *s = new TSpectrum(7);
s->Search(histo, 2,"new", 1);
Double_t *xpeaks = s->GetPositionX();
Double_t *ypeaks = s->GetPositionY();
cout << "Peaks X's and Y's :"<< endl ;
for (int i = 0 ; i < 7; i++){
cout << " X"<< i <<": " << xpeaks[i] << endl;
cout << " Y"<< i <<": " << ypeaks[i] << endl;
}
histo->Draw();
// create a TF1 with the range from 0 to 3 and 6 parameters
TF1 *fitFcn = new TF1("fitFcn",fitFunction, -180, 3200 , 25);
// Polynomial's parameters from background fitting
fitFcn->SetParLimits( 0, 343 , 419);
fitFcn->SetParLimits( 1, -1.18, -0.92);
fitFcn->SetParLimits( 2, 0.00122, 0.001463);
fitFcn->SetParLimits( 3, -2.83e-07, -2.23e-07);
// Setting MG parameters from par[4] to par[24]
fitFcn->FixParameter(0*3+4, s->GetPositionY()[6]); // Constant
fitFcn->SetParameter(0*3 + 1 +4, s->GetPositionX()[6]); // Mean
fitFcn->FixParameter(1*3+4, s->GetPositionY()[3]); // Constant
fitFcn->SetParameter(1*3 + 1+4, s->GetPositionX()[3]); // Mean
fitFcn->FixParameter(2*3+4, s->GetPositionY()[1]); // Constant
fitFcn->SetParameter(2*3 + 1 +4, s->GetPositionX()[1]); // Mean
fitFcn->FixParameter(3*3+4, s->GetPositionY()[0]); // Constant
fitFcn->SetParameter(3*3 + 1+4, s->GetPositionX()[0]); // Mean
fitFcn->FixParameter(4*3+4, s->GetPositionY()[2]); // Constant
fitFcn->SetParameter(4*3 + 1+4, s->GetPositionX()[2]); // Mean
fitFcn->FixParameter(5*3+4, s->GetPositionY()[4]); // Constant
fitFcn->SetParameter(5*3 + 1+4, s->GetPositionX()[4]); // Mean
fitFcn->FixParameter(6*3+4, s->GetPositionY()[5]); // Constant
fitFcn->SetParameter(6*3 + 1+4, s->GetPositionX()[5]); // Mean
fitFcn->SetParLimits(0*3 + 2+4, 27.10 , 32.90); // Width - rough order of mag estimate
fitFcn->SetParLimits(1*3 + 2+4, 33.64 , 41.44); // Width - rough order of mag estimate
fitFcn->SetParLimits(2*3 + 2+4, 40.47, 47.27); // Width - rough order of mag estimate
fitFcn->SetParLimits(3*3 + 2+4, 44.78 , 50.58); // Width - rough order of mag estimate
fitFcn->SetParLimits(4*3 + 2+4, 58.35, 65.15); // Width - rough order of mag estimate
fitFcn->SetParLimits(5*3 + 2+4, 77.2 ,85.3); // Width - rough order of mag estimate
fitFcn->SetParLimits(6*3 + 2+4, 188.2, 196.2);
// first try without starting values for the parameters
// this defaults to 1 for each param.
histo->Fit(fitFcn,"R+");
// histo->Fit("background");
//histo->Fit("MultiGaussianPeaks");
cout << "Chi^2:" << fitFcn->GetChisquare() << ", number of DoF: " << fitFcn->GetNDF() << " (Probability: " << fitFcn->GetProb() << ")." << endl;
}