]>
Commit | Line | Data |
---|---|---|
9e040a65 | 1 | #include <math.h> |
2 | #include <time.h> | |
3 | #include <stdio.h> | |
4 | #include <stdlib.h> | |
5 | #include <Riostream.h> | |
6 | #include <assert.h> | |
7 | ||
8 | #include "TVector2.h" | |
9 | #include "TFile.h" | |
10 | #include "TString.h" | |
11 | #include "TF1.h" | |
12 | #include "TF3.h" | |
13 | #include "TH1.h" | |
14 | #include "TH2.h" | |
15 | #include "TH3.h" | |
16 | #include "TProfile.h" | |
17 | #include "TProfile2D.h" | |
18 | #include "TMath.h" | |
19 | #include "TText.h" | |
20 | #include "TRandom3.h" | |
21 | #include "TArray.h" | |
22 | #include "TLegend.h" | |
23 | #include "TStyle.h" | |
24 | #include "TMinuit.h" | |
25 | #include "TCanvas.h" | |
26 | #include "TLatex.h" | |
27 | #include "TPaveStats.h" | |
5130881c | 28 | #include "TASImage.h" |
0e610a08 | 29 | #include "TSpline.h" |
30 | ||
9e040a65 | 31 | |
32 | #define BohrR 1963.6885 // Bohr Radius for pions | |
33 | #define FmToGeV 0.19733 // conversion to fm | |
34 | #define PI 3.1415926 | |
35 | #define masspiC 0.1395702 // pi+ mass (GeV/c^2) | |
0e610a08 | 36 | double kappa3 = 0.1; // 0.16 (default), 0.05 or 0.3 as variations |
37 | double kappa4 = 0.5; // 0.4 (default), | |
9e040a65 | 38 | |
39 | using namespace std; | |
40 | ||
0e610a08 | 41 | int CollisionType_def=1;// 0=PbPb, 1=pPb, 2=pp |
42 | bool MCcase_def=0;// MC data? | |
9e040a65 | 43 | int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case |
b2c0d4ba | 44 | bool SameCharge_def=kTRUE;// 3-pion same-charge? |
9e040a65 | 45 | bool AddCC=kTRUE; |
0e610a08 | 46 | bool Gaussian=1;// Gaussian or Exponential? |
47 | bool UseC2Bkg=1;// use Pythia/DPMJET bkg? | |
9e040a65 | 48 | // |
5be394f6 | 49 | bool MuonCorrection=1;// apply Muon correction? |
0e610a08 | 50 | bool IncludeExpansion=kTRUE;// Include EdgeWorth coefficients? |
51 | bool FixExpansionAvg=1; | |
52 | int Mbin_def=18;// 0-19 (0=1050-2000 pions, 19=0-5 pions) | |
9e040a65 | 53 | int EDbin_def=0;// 0-2: Kt3 bin |
0e610a08 | 54 | int Ktbin_def=1;// 1-6, 10=full range |
9e040a65 | 55 | double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase |
0e610a08 | 56 | bool FitRangeShift=0;// 30% reduction in Fit Range |
9e040a65 | 57 | // |
58 | // | |
59 | bool SaveToFile_def=kFALSE;// Save outputs to file? | |
60 | int SourceType=0;// 0=Therminator, 1=Therminator with 50fm cut (keep at 0 for default), 2=Gaussian | |
61 | bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin? | |
62 | bool GofP=kFALSE;// Include momentum dependence of coherent fraction? | |
63 | bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states? | |
64 | bool LinkN=kFALSE;// Make N(++)=N(+-)? | |
65 | bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC? | |
66 | bool Tricubic=kFALSE;// Tricubic or Trilinear interpolation? Tricubic was found to be unstable. | |
67 | bool IncludeSS=kTRUE;// Include same-charge two-pion correlations in the fit? | |
68 | bool IncludeOS=kFALSE;// Include mixed-charge two-pion correlations in the fit? | |
69 | // | |
70 | // | |
0e610a08 | 71 | double ExpPower=2;// default Gaussian(2), Exponential(1) |
9e040a65 | 72 | // |
73 | // | |
74 | const int Sbins_2=1;// 2-particle Species bins. 1=Pion-Pion only. max=6 (pi-pi, pi-k, pi-p, k-p, k-k, p-p) | |
75 | const int Sbins_3=1;// 3-particle Species bins. 1=Pion-Pion-Pion only. max=10 | |
76 | const int fitlimitLowBin = 3; | |
77 | const int fitlimitHighBin = 14;// max = 20 (100MeV) | |
78 | bool LHC11h=kTRUE;// always kTRUE | |
79 | bool ExplicitLoop=kFALSE;// always kFALSE | |
80 | bool PairCut=kTRUE;// always kTRUE | |
81 | int FSIindex; | |
82 | int Ktlowbin; | |
83 | int Kthighbin; | |
84 | int MomResCentBin;// momentum resolution centrality bin | |
85 | float TwoFrac;// Lambda parameter | |
86 | float OneFrac;// Lambda^{1/2} | |
87 | float ThreeFrac;// Lambda^{3/2} | |
88 | float Q2Limit; | |
89 | float Q3Limit; | |
90 | ||
0e610a08 | 91 | double lambda_PM = 0.7; |
92 | ||
7c6e4906 | 93 | |
9e040a65 | 94 | TH1D *CoulCorr2SS; |
95 | TH1D *CoulCorr2OS; | |
96 | TH1D *CoulCorr2SS_2;// for interpolation in transition from Therminator to Gauss K factors | |
97 | TH1D *CoulCorr2OS_2;// for interpolation in transition from Therminator to Gauss K factors | |
98 | // | |
99 | ||
100 | void ReadCoulCorrections(int); | |
101 | void ReadMomResFile(int); | |
5be394f6 | 102 | void ReadMuonCorrections(int,int); |
9e040a65 | 103 | double CoulCorr2(int, double); |
104 | double CoulCorrGRS(bool, double, double, double); | |
105 | double Dfitfunction_c3(Double_t *x, Double_t *par); | |
106 | double Gamov(int, double); | |
107 | double C2ssFitFunction(Double_t *x, Double_t *par); | |
108 | double C2osFitFunction(Double_t *x, Double_t *par); | |
109 | double cubicInterpolate(double[4], double); | |
110 | double bicubicInterpolate(double[4][4], double, double); | |
111 | double tricubicInterpolate(double[4][4][4], double, double, double); | |
5130881c | 112 | void DrawALICELogo(Bool_t, Float_t, Float_t, Float_t, Float_t); |
113 | ||
9e040a65 | 114 | // |
115 | void fcnC2_Global(int&, double*, double&, double[], int); | |
116 | ||
7c6e4906 | 117 | const int BINRANGE_C2global=400;// 40 |
9e040a65 | 118 | double C2ssFitting[BINRANGE_C2global]; |
119 | double C2osFitting[BINRANGE_C2global]; | |
120 | double C2ssFitting_e[BINRANGE_C2global]; | |
121 | double C2osFitting_e[BINRANGE_C2global]; | |
122 | double K2SS[BINRANGE_C2global]; | |
123 | double K2OS[BINRANGE_C2global]; | |
124 | double K2SS_e[BINRANGE_C2global]; | |
125 | double K2OS_e[BINRANGE_C2global]; | |
126 | double Chi2_C2global; | |
127 | double NFitPoints_C2global; | |
128 | double Chi2_C3global; | |
129 | double NFitPoints_C3global; | |
130 | ||
131 | const int BINS_OSL=40;// out,side,long bins | |
132 | double A[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
133 | double B[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
134 | double A_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
135 | double B_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
136 | double avg_q[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
137 | double K_OSL[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
138 | double K_OSL_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long | |
139 | double Chi2_OSL; | |
140 | double NFitPoints_OSL; | |
141 | ||
142 | const int BINRANGE_3=40;// q12,q13,q23 | |
143 | int BINLIMIT_3=20; | |
144 | double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 | |
145 | double A_3_e[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 | |
146 | double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 | |
147 | double BinCenters[400]; | |
148 | double BinWidthQ2; | |
149 | double Chi2_c3; | |
150 | double NFitPoints_c3; | |
151 | void fcn_c3(int&, double*, double&, double[], int); | |
152 | ||
153 | double C3_base[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
154 | double N_term1[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
155 | double N_term2[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
156 | double N_term3[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
157 | double N_term4[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
158 | double N_term5[BINRANGE_3][BINRANGE_3][BINRANGE_3]; | |
159 | ||
160 | ||
161 | TH3D *MomRes3d[2][5];// SC/MC, term# | |
162 | TH1D *MomRes1d[2][5];// SC/MC, term# | |
163 | TH1D *MomResC2[2];// SC/MC | |
164 | int MRC2index;// index for C2 MRC | |
165 | ||
166 | double AvgP[6][20];// kt bin, qinv bin | |
167 | ||
5be394f6 | 168 | TH1D *C2muonCorrection_SC; |
169 | TH1D *C2muonCorrection_MC; | |
170 | TH1D *C3muonCorrection; | |
9e040a65 | 171 | |
172 | TF1 *StrongSC;// same-charge pion strong FSI | |
173 | TF1 *MixedChargeSysFit;// mixed-charge 3-pion cumulant residue obtained from Plot_plotsTPR.C | |
0e610a08 | 174 | TF1 *BkgMCFit;// Pythia or DPMJET fit |
9e040a65 | 175 | // |
176 | ||
177 | ||
b2c0d4ba | 178 | void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, bool MCcase=MCcase_def, bool SameCharge=SameCharge_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){ |
9e040a65 | 179 | |
180 | EDbin_def=EDbin; | |
181 | Ktbin_def=Ktbin; | |
182 | CollisionType_def=CollisionType; | |
183 | SaveToFile_def=SaveToFile; | |
184 | MCcase_def=MCcase; | |
185 | CHARGE_def=CHARGE; | |
9e040a65 | 186 | SameCharge_def=SameCharge;// 3-pion same-charge |
187 | Mbin_def=Mbin; | |
188 | // | |
5be394f6 | 189 | |
0e610a08 | 190 | // PbPb then pPb then pp |
191 | //double Nch_means[3][20]={{1828,1407,1071,896,757,619,496,402,324,252,171,117,83,63,50,36, 36,36,36,36}, {71,71,71,71,71,71,71,71,71,71,71,71, 71,57,45,32,23,16,9.7,5}, {47,47,47,47,47,47,47,47,47,47,47,47,47, 47,37,27,20,15,8.6,4.6}}; | |
192 | ||
193 | kappa3 = 0.1; | |
194 | kappa4 = 0.5; | |
195 | //kappa3 = 0.05 + 0.25*(Mbin/19.);// linear dependence of kappa3 | |
196 | ||
197 | //double MainFitParams[8]={0}; | |
5be394f6 | 198 | // |
0e610a08 | 199 | //double RefMainFitParams[8]={3.58014, 2.84635, 0.652247, 0.831958, 3.13505, 2.33717, 0.761147, 0.909211}; |
9e040a65 | 200 | |
0e610a08 | 201 | if(Gaussian) ExpPower=2.; |
202 | else ExpPower=1.; | |
9e040a65 | 203 | |
0e610a08 | 204 | // old way v5 and before |
205 | /*if(Mbin==0) FSIindex = 0;// 0-5% Therm | |
5130881c | 206 | else if(Mbin==1) FSIindex = 1;// 0-10% Therm |
207 | else if(Mbin<=3) FSIindex = 2;// 10-20% Therm | |
208 | else if(Mbin<=5) FSIindex = 3;// 20-30% Therm | |
209 | else if(Mbin<=7) FSIindex = 4;// 30-40% Therm | |
210 | else if(Mbin<=9) FSIindex = 5;// 40-50% Therm | |
211 | else if(Mbin<=11) FSIindex = 6;// Gauss R=3.5 fm | |
212 | else if(Mbin<=14) FSIindex = 7;// Gauss R=3.0 fm | |
213 | else if(Mbin<=16) FSIindex = 8;// Gauss R=2.5 fm | |
0e610a08 | 214 | else FSIindex = 9;// Gauss R=2.0 fm*/ |
215 | ||
216 | if(CollisionType==0){ | |
217 | if(Mbin==0) FSIindex = 0;// 0-5% Therm | |
218 | else if(Mbin==1) FSIindex = 1;// 0-10% Therm | |
219 | else if(Mbin<=3) FSIindex = 2;// 10-20% Therm | |
220 | else if(Mbin<=5) FSIindex = 3;// 20-30% Therm | |
221 | else if(Mbin<=7) FSIindex = 4;// 30-40% Therm | |
222 | else if(Mbin<=9) FSIindex = 5;// 40-50% Therm | |
223 | else FSIindex = 6;// EW R=4.0 fm | |
224 | // next bin is 3.0 fm which is used for interpolation | |
225 | }else{ | |
226 | if(Mbin<=14) FSIindex = 8;// EW R=2.5 fm | |
227 | else if(Mbin<=16) FSIindex = 9;// EW R=2.0 fm | |
228 | else FSIindex = 10;// EW R=1.7 fm | |
229 | // last bin is 1.25 fm which is used for interpolation | |
230 | } | |
5130881c | 231 | |
232 | // chooses lambda=0.7 for R=10,8,6,4,2 | |
9e040a65 | 233 | if(Mbin<=2) MRC2index=138; |
234 | else if(Mbin<=5) MRC2index=106; | |
235 | else if(Mbin<=9) MRC2index=74; | |
5130881c | 236 | else if(Mbin<=13) MRC2index=42;// was 9 by mistake, now 13 (corrected on 22/10/2013) |
9e040a65 | 237 | else MRC2index=10; |
238 | ||
239 | ||
240 | ||
0e610a08 | 241 | TwoFrac = lambda_PM; |
9e040a65 | 242 | OneFrac = sqrt(TwoFrac); |
243 | ThreeFrac = pow(TwoFrac,3/2.); | |
244 | ||
245 | if(CollisionType==0 && Mbin<10) BINLIMIT_3=20; | |
246 | else BINLIMIT_3=40; | |
247 | ||
248 | Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...) | |
249 | Kthighbin=(Ktbin)*2+4; | |
5be394f6 | 250 | if(Ktbin==2) {cout<<"Kthighbin increased"<<endl; Kthighbin = 20;}// to make <KT3> comparible to <kT> |
9e040a65 | 251 | |
0e610a08 | 252 | |
253 | BkgMCFit=new TF1("BkgMCFit","[0]*([1] + [2]*exp(-pow([3]*x,2)))",0,2.0); | |
254 | BkgMCFit->SetParameter(0,1); | |
255 | BkgMCFit->SetParameter(1,1.001); | |
256 | BkgMCFit->SetParameter(2,0.04102); | |
257 | BkgMCFit->SetParameter(3,1.874); | |
258 | BkgMCFit->SetLineColor(1); | |
7c6e4906 | 259 | if(!MCcase){ |
0e610a08 | 260 | // kt bin(0.2<kt<0.3, 0.3<kt<1.0), Mult bin(low to high) |
261 | float paramsPP_0[2][10]={{1.00014,1.00019,9.99733e-01,9.99379e-01,9.99488e-01, 9.99488e-01,9.99488e-01,9.99488e-01,9.99488e-01,9.99488e-01},{9.97267e-01,9.98918e-01,9.99117e-01,9.98878e-01,9.99457e-01, 9.99457e-01,9.99457e-01,9.99457e-01,9.99457e-01,9.99457e-01}}; | |
262 | float paramsPP_1[2][10]={{1.00198,1.00156,1.00109,1.00068,1.00079, 1.00079,1.00079,1.00079,1.00079,1.00079}, | |
263 | {9.98091e-01,9.98893e-01,9.99782e-01,9.99595e-01,1.00037, 1.00037,1.00037,1.00037,1.00037,1.00037}}; | |
264 | float paramsPP_2[2][10]={{4.51532e-02,4.30415e-02,4.10306e-02,3.42739e-02,2.95607e-02, 2.95607e-02,2.95607e-02,2.95607e-02,2.95607e-02,2.95607e-02}, | |
265 | {9.89882e-02,1.00657e-01,8.93560e-02,7.55225e-02,6.52689e-02, 6.52689e-02,6.52689e-02,6.52689e-02,6.52689e-02,6.52689e-02}}; | |
266 | float paramsPP_3[2][10]={{1.57172,1.91807,1.87353,1.88549,2.22286, 2.22286,2.22286,2.22286,2.22286,2.22286}, | |
267 | {1.49834,1.65562,1.66501,1.75097,1.78597, 1.78597,1.78597,1.78597,1.78597,1.78597}}; | |
268 | ////////////// | |
269 | float paramsPPb_0[2][10]={{9.99479e-01,9.98772e-01,9.99295e-01,9.99532e-01,9.99430e-01,9.99642e-01,9.99424e-01, 9.99424e-01,9.99424e-01,9.99424e-01}, | |
270 | {9.91941e-01,9.96569e-01,9.97183e-01,9.98106e-01,9.98591e-01,9.98717e-01,9.98517e-01, 9.98517e-01,9.98517e-01,9.98517e-01}}; | |
271 | float paramsPPb_1[2][10]={{1.00093,1.00023,1.00068,1.00082,1.00077,1.00087,1.00069, 1.00069,1.00069,1.00069}, | |
272 | {9.93727e-01,9.98118e-01,9.98630e-01,9.99584e-01,1.00005,1.00017,9.99972e-01, 9.99972e-01,9.99972e-01,9.99972e-01}}; | |
273 | float paramsPPb_2[2][10]={{3.05502e-02,2.21034e-02,1.39907e-02,9.82480e-03,6.48672e-03,5.37259e-03,4.54498e-03, 4.54498e-03,4.54498e-03,4.54498e-03}, | |
274 | {8.86336e-02,5.53244e-02,3.59877e-02,2.57145e-02,1.82072e-02,1.44253e-02,1.28923e-02, 1.28923e-02,1.28923e-02, 1.28923e-02}}; | |
275 | float paramsPPb_3[2][10]={{1.69217,1.51289,1.63134,1.71567,1.71918,1.78037,2.10736, 2.10736,2.10736,2.10736}, | |
276 | {1.21996,1.34570,1.35660,1.37801,1.39951,1.39093,1.33747, 1.33747,1.33747,1.33747}}; | |
277 | ||
278 | int ktindexBkg = Ktbin - 1; | |
279 | if(ktindexBkg > 1) ktindexBkg = 1; | |
280 | if(CollisionType==2){ | |
281 | BkgMCFit->FixParameter(0, paramsPP_0[ktindexBkg][19-Mbin]); | |
282 | BkgMCFit->FixParameter(1, paramsPP_1[ktindexBkg][19-Mbin]); | |
283 | BkgMCFit->FixParameter(2, paramsPP_2[ktindexBkg][19-Mbin]); | |
284 | BkgMCFit->FixParameter(3, paramsPP_3[ktindexBkg][19-Mbin]); | |
285 | }else if(CollisionType==1){ | |
286 | BkgMCFit->FixParameter(0, paramsPPb_0[ktindexBkg][19-Mbin]); | |
287 | BkgMCFit->FixParameter(1, paramsPPb_1[ktindexBkg][19-Mbin]); | |
288 | BkgMCFit->FixParameter(2, paramsPPb_2[ktindexBkg][19-Mbin]); | |
289 | BkgMCFit->FixParameter(3, paramsPPb_3[ktindexBkg][19-Mbin]); | |
290 | }else{// no Bkg for PbPb | |
291 | BkgMCFit->FixParameter(0, 1); | |
292 | BkgMCFit->FixParameter(1, 1); | |
293 | BkgMCFit->FixParameter(2, 0); | |
294 | BkgMCFit->FixParameter(3, 1); | |
295 | } | |
296 | ||
7c6e4906 | 297 | } |
9e040a65 | 298 | // Core-Halo modeling of single-particle and triple-particle core fraction |
299 | //float AvgN[10]={1.99266, 3.97789, 6.4624, 8.94042, 11.4194, 13.8987, 16.385, 18.8756, 21.3691, 26.0742};// 13c (avg total pion mult)/2. 2.0sigma | |
9e040a65 | 300 | // |
301 | ||
302 | // bin centers from QS+FSI | |
303 | double BinCenterPbPbCentral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509}; | |
304 | double BinCenterPbPbPeripheral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509}; | |
305 | double BinCenterpPbAndpp[40]={0.00359275, 0.016357, 0.0257109, 0.035445, 0.045297, 0.0552251, 0.0651888, 0.0751397, 0.0851088, 0.0951108, 0.105084, 0.115079, 0.12507, 0.135059, 0.145053, 0.155049, 0.16505, 0.175038, 0.185039, 0.195034, 0.205023, 0.215027, 0.225024, 0.235023, 0.245011, 0.255017, 0.265017, 0.275021, 0.285021, 0.295017, 0.305018, 0.315018, 0.325013, 0.335011, 0.345016, 0.355019, 0.365012, 0.375016, 0.385017, 0.395016}; | |
306 | if(CollisionType==0){ | |
307 | Q3Limit = 0.1 + 0.1*Mbin/16.; | |
308 | for(int i=0; i<40; i++){ | |
309 | if(Mbin<10) BinCenters[i] = BinCenterPbPbCentral[i]; | |
310 | else BinCenters[i] = BinCenterPbPbPeripheral[i]; | |
311 | } | |
312 | BinWidthQ2 = 0.005; | |
313 | }else{ | |
314 | Q3Limit = 0.3 + 0.2*fabs(Mbin-10)/9.; | |
315 | for(int i=0; i<40; i++){ | |
316 | BinCenters[i] = BinCenterpPbAndpp[i]; | |
317 | } | |
318 | BinWidthQ2 = 0.01; | |
319 | } | |
5be394f6 | 320 | if(FitRangeShift) Q3Limit -= 0.3*Q3Limit; |
9e040a65 | 321 | Q2Limit = Q3Limit/sqrt(2.); |
0e610a08 | 322 | //Q2Limit = 1.2; |
323 | //cout<<Q2Limit<<" "<<Q3Limit<<endl; | |
5130881c | 324 | |
9e040a65 | 325 | // extend BinCenters for high q |
326 | for(int index=40; index<400; index++){ | |
327 | if(CollisionType==0) BinCenters[index] = (index+0.5)*(0.005); | |
328 | else BinCenters[index] = (index+0.5)*(0.010); | |
329 | } | |
b2c0d4ba | 330 | // Set 0's to 3-particle fit arrays |
331 | for(int i=1; i<=BINLIMIT_3; i++){// bin number | |
332 | for(int j=1; j<=BINLIMIT_3; j++){// bin number | |
333 | for(int k=1; k<=BINLIMIT_3; k++){// bin number | |
334 | A_3[i-1][j-1][k-1]=0; | |
335 | A_3_e[i-1][j-1][k-1]=0; | |
336 | B_3[i-1][j-1][k-1]=0; | |
337 | } | |
338 | } | |
339 | } | |
9e040a65 | 340 | |
341 | // same-charge pion strong FSI (obtained from ratio of CoulStrong to Coul at R=7 Gaussian from Lednicky's code) | |
342 | StrongSC=new TF1("StrongSC","[0]+[1]*exp(-pow([2]*x,2))",0,1000); | |
343 | StrongSC->FixParameter(0,9.99192e-01); | |
344 | StrongSC->FixParameter(1,-8.01113e-03); | |
345 | StrongSC->FixParameter(2,5.35912e-02); | |
346 | // mixed-charge 3-pion cumulant residue obtained from Plot_plotsTPR.C | |
347 | MixedChargeSysFit=new TF1("MixedChargeSysFit","[0]+[1]*exp(-pow([2]*x/0.19733,2))",0,.5); | |
348 | MixedChargeSysFit->FixParameter(0,1); | |
349 | MixedChargeSysFit->FixParameter(1,.06); | |
350 | MixedChargeSysFit->FixParameter(2,5.5 - 6*(Mbin/19.)); | |
351 | ||
352 | ||
353 | gStyle->SetOptStat(0); | |
354 | gStyle->SetOptDate(0); | |
355 | gStyle->SetOptFit(0111); | |
356 | ||
5130881c | 357 | TFile *_file0; |
9e040a65 | 358 | if(CollisionType==0){// PbPb |
359 | if(MCcase) { | |
360 | if(Mbin<=1){ | |
b2c0d4ba | 361 | _file0 = new TFile("Results/PDC_12a17a_R10_2.root","READ"); |
9e040a65 | 362 | }else{ |
b2c0d4ba | 363 | _file0 = new TFile("Results/PDC_12a17e_R10.root","READ"); |
9e040a65 | 364 | } |
365 | }else{ | |
5130881c | 366 | //_file0 = new TFile("Results/PDC_10h_11h_0to50_50to100.root","READ"); |
7c6e4906 | 367 | //_file0 = new TFile("Results/PDC_10h_11h_0to50_50to100_3Ktbins.root","READ"); |
b2c0d4ba | 368 | //_file0 = new TFile("Results/PDC_10h_2Kt3bins.root","READ"); |
0e610a08 | 369 | //_file0 = new TFile("Results/PDC_10h_NewNorm.root","READ"); |
370 | //_file0 = new TFile("Results/PDC_10h_dEta0p03dPhi0p03.root","READ"); | |
b2c0d4ba | 371 | //_file0 = new TFile("Results/PDC_10h_noPID.root","READ"); |
372 | //_file0 = new TFile("Results/PDC_10h_1percentCentral.root","READ"); | |
0e610a08 | 373 | //_file0 = new TFile("Results/PDC_10h_11h_2Kt3bins.root","READ");// v5 and before |
5be394f6 | 374 | //_file0 = new TFile("Results/PDC_10h_dEta0p025dPhi0p055.root","READ"); |
5130881c | 375 | //_file0 = new TFile("Results/PDC_11h_3Kt3bins_FB5and7overlap.root","READ"); |
376 | //_file0 = new TFile("Results/PDC_10h_11h_V0binning.root","READ"); | |
0e610a08 | 377 | _file0 = new TFile("Results/PDC_10h_11h_NclsFix.root","READ");// standard |
9e040a65 | 378 | } |
379 | }else if(CollisionType==1){// pPb | |
380 | if(!MCcase){ | |
5130881c | 381 | //_file0 = new TFile("Results/PDC_13c_kMB.root","READ"); |
382 | //_file0 = new TFile("Results/PDC_13bc.root","READ"); | |
383 | //_file0 = new TFile("Results/PDC_13bc_3Ktbins.root","READ"); | |
7c6e4906 | 384 | //_file0 = new TFile("Results/PDC_13bc_kAnyINT_kINT7.root","READ"); |
385 | //_file0 = new TFile("Results/PDC_13bc_2Kt3bins.root","READ"); | |
5be394f6 | 386 | //_file0 = new TFile("Results/PDC_13bc_kINT7.root","READ");// paper v1 |
387 | //_file0 = new TFile("Results/PDC_13c_posEta.root","READ");// eta sub-range | |
388 | //_file0 = new TFile("Results/PDC_13c_negEta.root","READ");// eta sub-range | |
389 | //_file0 = new TFile("Results/PDC_13c_0p4to0p8eta.root","READ");// eta sub-range | |
390 | //_file0 = new TFile("Results/PDC_13c_neg0p4toneg0p8eta.root","READ");// eta sub-range | |
391 | //_file0 = new TFile("Results/PDC_13bc_dEta0p025dPhi0p055_kINT7_kHighMult.root","READ"); | |
5130881c | 392 | //_file0 = new TFile("Results/PDC_13bc_FB5and7overlap_V0binning.root","READ"); |
7c6e4906 | 393 | //_file0 = new TFile("Results/PDC_13bc_kAnyINT_V0binning.root","READ"); |
5be394f6 | 394 | //_file0 = new TFile("Results/PDC_13bcd_kHighMult.root","READ"); |
395 | //_file0 = new TFile("Results/PDC_13bcd_kINT7.root","READ"); | |
0e610a08 | 396 | //_file0 = new TFile("Results/PDC_13f_kHighMult.root","READ"); |
397 | //_file0 = new TFile("Results/PDC_13e_kHighMult.root","READ"); | |
398 | //_file0 = new TFile("Results/PDC_13bcd_kINT7_plus_kHighMult.root","READ");// v5 and before | |
399 | //_file0 = new TFile("Results/PDC_13bcde_kINT7_NclsFix.root","READ"); | |
400 | _file0 = new TFile("Results/PDC_13bcde_kINT7_kHighMult_NclsFix.root","READ");// standard | |
9e040a65 | 401 | }else{ |
0e610a08 | 402 | _file0 = new TFile("Results/PDC_13b2_efix_p1234_R2_2Ktbins.root","READ");// standard (gen level) |
403 | //_file0 = new TFile("Results/PDC_13b2_efix_p1234_R2_RecLevel.root","READ");// Reconstructed tracks | |
404 | //_file0 = new TFile("Results/PDC_13b2_p1_noTTC.root","READ"); | |
9e040a65 | 405 | } |
406 | }else{// pp | |
407 | if(!MCcase){ | |
5130881c | 408 | //_file0 = new TFile("Results/PDC_10cde.root","READ"); |
409 | //_file0 = new TFile("Results/PDC_10cde_3Ktbins.root","READ"); | |
410 | //_file0 = new TFile("Results/PDC_10e_kHighMult.root","READ"); | |
411 | //_file0 = new TFile("Results/PDC_10e_kAnyINT.root","READ"); | |
412 | //_file0 = new TFile("Results/PDC_10cde_kAnyINT.root","READ"); | |
7c6e4906 | 413 | //_file0 = new TFile("Results/PDC_10cde_kAnyINT_Plus_kHighMult.root","READ"); |
5be394f6 | 414 | //_file0 = new TFile("Results/K0sRange_10cde.root","READ");// for K0s peak removal |
7c6e4906 | 415 | //_file0 = new TFile("Results/PDC_10cde_2Kt3bins.root","READ"); |
5be394f6 | 416 | //_file0 = new TFile("Results/PDC_10e_kHighMult.root","READ"); |
b2c0d4ba | 417 | //_file0 = new TFile("Results/PDC_10cde_kMB.root","READ"); |
0e610a08 | 418 | //_file0 = new TFile("Results/PDC_10cde_kMB_plus_kHighMult_2Kt3bins.root","READ");// v2 paper |
419 | //_file0 = new TFile("Results/PDC_10cde_kMB_plus_kHighMult_AOD147.root","READ");// v5 and before | |
5130881c | 420 | //_file0 = new TFile("Results/PDC_10cde_FB5and7overlap.root","READ"); |
0e610a08 | 421 | //_file0 = new TFile("Results/PDC_10d_noTTC.root","READ"); |
422 | //_file0 = new TFile("Results/PDC_10cde_kMB.root","READ"); | |
423 | _file0 = new TFile("Results/PDC_10cde_kMB_kHighMult_NclsFix.root","READ");// standard | |
9e040a65 | 424 | }else{ |
0e610a08 | 425 | _file0 = new TFile("Results/PDC_10f6a_R2_2Ktbins.root","READ");// standard |
426 | //_file0 = new TFile("Results/PDC_10f6a_R10_genLevel_2Ktbins.root","READ"); | |
9e040a65 | 427 | } |
428 | } | |
429 | ||
430 | ReadCoulCorrections(FSIindex); | |
431 | ReadMomResFile(Mbin); | |
5be394f6 | 432 | ReadMuonCorrections(CollisionType,Mbin); |
9e040a65 | 433 | // |
434 | ///////////////////////////////////////////////////////// | |
435 | ||
436 | ||
437 | double NormQcutLow; | |
438 | double NormQcutHigh; | |
439 | if(CollisionType==0) {// PbPb | |
440 | if(Mbin<10){ | |
0e610a08 | 441 | NormQcutLow = 0.15;// was 0.15 |
442 | NormQcutHigh = 0.175;// was 0.175 | |
9e040a65 | 443 | }else{ |
0e610a08 | 444 | NormQcutLow = 0.3;// was 0.3 |
445 | NormQcutHigh = 0.35;// was 0.35 | |
9e040a65 | 446 | } |
447 | }else{// pPb and pp | |
0e610a08 | 448 | NormQcutLow = 1.0;// was 1.0 |
449 | NormQcutHigh = 1.2;// was 1.2 | |
9e040a65 | 450 | } |
451 | ||
452 | ||
453 | TList *MyList; | |
454 | TDirectoryFile *tdir; | |
5130881c | 455 | if(!MCcase) tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputThreePionRadiiAnalysis.root"); |
9e040a65 | 456 | |
0e610a08 | 457 | |
9e040a65 | 458 | if(CollisionType==0){ |
459 | if(!MCcase){ | |
460 | if(Mbin<6) MyList=(TList*)tdir->Get("ThreePionRadiiOutput_1"); | |
461 | else MyList=(TList*)tdir->Get("ThreePionRadiiOutput_2"); | |
462 | }else{ | |
5130881c | 463 | MyList=(TList*)_file0->Get("MyList"); |
9e040a65 | 464 | } |
465 | }else { | |
5be394f6 | 466 | if(!MCcase) { |
5be394f6 | 467 | if(CollisionType==2 && Mbin<15) MyList=(TList*)tdir->Get("ThreePionRadiiOutput_2"); |
468 | else MyList=(TList*)tdir->Get("ThreePionRadiiOutput_1"); | |
469 | }else MyList=(TList*)_file0->Get("MyList"); | |
9e040a65 | 470 | } |
5130881c | 471 | _file0->Close(); |
9e040a65 | 472 | |
473 | TH1D *Events = (TH1D*)MyList->FindObject("fEvents2"); | |
474 | // | |
9e040a65 | 475 | cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl; |
476 | ||
0e610a08 | 477 | // below lines determine fractional cross-sections without correction for efficiency |
478 | //TH1F *fMultDist3 = (TH1F*)MyList->FindObject("fMultDist3");// total event count before track cuts | |
479 | //for(int mm=19; mm>=10; mm--){ | |
480 | //cout<<"fracton of cross-section = "<<Events->GetBinContent(mm+1)/fMultDist3->GetEntries()<<endl; | |
481 | //} | |
482 | ||
483 | TF1 *Unity = new TF1("Unity","1",0,10); | |
484 | Unity->SetLineStyle(2); Unity->SetLineColor(1); | |
9e040a65 | 485 | |
9e040a65 | 486 | // Explicit Loop Histos |
487 | TH2D *Two_ex_2d[2][2][Sbins_2][2]; | |
488 | TH1D *Two_ex[2][2][Sbins_2][2]; | |
489 | ||
490 | // Normalizations | |
491 | double NormH_pc[5]={0}; | |
492 | ||
493 | double norm_ex_2[6][2]={{0}}; | |
494 | ||
495 | // 3d method histograms | |
496 | TH3D *Three_3d[2][2][2][Sbins_3][5]; | |
497 | TH1D *Three_1d[2][2][2][Sbins_3][5]; | |
498 | /////////////////////////////////// | |
499 | // Create Histograms | |
500 | ||
501 | // 2-particle | |
502 | for(int c1=0; c1<2; c1++){ | |
503 | for(int c2=0; c2<2; c2++){ | |
504 | for(int sc=0; sc<Sbins_2; sc++){ | |
505 | for(int term=0; term<2; term++){ | |
506 | ||
507 | TString *name2_ex = new TString("Explicit2_Charge1_"); | |
508 | *name2_ex += c1; | |
509 | name2_ex->Append("_Charge2_"); | |
510 | *name2_ex += c2; | |
511 | name2_ex->Append("_SC_"); | |
512 | *name2_ex += sc; | |
513 | name2_ex->Append("_M_"); | |
514 | *name2_ex += Mbin; | |
515 | name2_ex->Append("_ED_"); | |
516 | *name2_ex += 0;// EDbin | |
517 | name2_ex->Append("_Term_"); | |
518 | *name2_ex += term+1; | |
519 | ||
520 | if(sc==0 || sc==3 || sc==5){ | |
521 | if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram | |
522 | } | |
523 | ||
524 | Two_ex_2d[c1][c2][sc][term] = (TH2D*)MyList->FindObject(name2_ex->Data()); | |
525 | Two_ex_2d[c1][c2][sc][term]->Sumw2(); | |
526 | TString *proname = new TString(name2_ex->Data()); | |
527 | proname->Append("_pro"); | |
528 | ||
529 | if(Ktbin==10) {Ktlowbin=1; Kthighbin=Two_ex_2d[c1][c2][sc][term]->GetNbinsX();}//full kt interval | |
530 | Two_ex[c1][c2][sc][term] = (TH1D*)Two_ex_2d[c1][c2][sc][term]->ProjectionY(proname->Data(),Ktlowbin,Kthighbin);// bins:5-6 (kt:0.2-0.3) | |
531 | ||
532 | norm_ex_2[sc][term] = Two_ex[c1][c2][sc][term]->Integral(Two_ex[c1][c2][sc][term]->GetXaxis()->FindBin(NormQcutLow),Two_ex[c1][c2][sc][term]->GetXaxis()->FindBin(NormQcutHigh)); | |
533 | Two_ex[c1][c2][sc][term]->Scale(norm_ex_2[sc][0]/norm_ex_2[sc][term]); | |
534 | ||
535 | Two_ex[c1][c2][sc][term]->SetMarkerStyle(20); | |
5130881c | 536 | Two_ex[c1][c2][sc][term]->GetXaxis()->SetTitle("#font[12]{q}_{inv}"); |
537 | Two_ex[c1][c2][sc][term]->GetYaxis()->SetTitle("#font[12]{C}_{2}"); | |
9e040a65 | 538 | Two_ex[c1][c2][sc][term]->SetTitle(""); |
539 | ||
540 | }// term | |
541 | }// sc | |
542 | ||
543 | // 3-particle | |
544 | for(int c3=0; c3<2; c3++){ | |
545 | for(int sc=0; sc<Sbins_3; sc++){ | |
546 | for(int term=0; term<5; term++){ | |
547 | ||
548 | TString *name3_ex = new TString("Explicit3_Charge1_"); | |
549 | *name3_ex += c1; | |
550 | name3_ex->Append("_Charge2_"); | |
551 | *name3_ex += c2; | |
552 | name3_ex->Append("_Charge3_"); | |
553 | *name3_ex += c3; | |
554 | name3_ex->Append("_SC_"); | |
555 | *name3_ex += sc; | |
556 | name3_ex->Append("_M_"); | |
557 | *name3_ex += Mbin; | |
558 | name3_ex->Append("_ED_"); | |
559 | *name3_ex += EDbin; | |
560 | name3_ex->Append("_Term_"); | |
561 | *name3_ex += term+1; | |
562 | ||
563 | ||
564 | TString *name3_pc = new TString("PairCut3_Charge1_"); | |
565 | *name3_pc += c1; | |
566 | name3_pc->Append("_Charge2_"); | |
567 | *name3_pc += c2; | |
568 | name3_pc->Append("_Charge3_"); | |
569 | *name3_pc += c3; | |
570 | name3_pc->Append("_SC_"); | |
571 | *name3_pc += sc; | |
572 | name3_pc->Append("_M_"); | |
573 | *name3_pc += Mbin; | |
574 | name3_pc->Append("_ED_"); | |
575 | *name3_pc += EDbin; | |
576 | name3_pc->Append("_Term_"); | |
577 | *name3_pc += term+1; | |
578 | ||
579 | /////////////////////////////////////// | |
580 | // skip degenerate histograms | |
581 | if(sc==0 || sc==6 || sc==9){// Identical species | |
582 | if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;} | |
583 | if( (c1+c2+c3)==2) {if(c1!=0) continue;} | |
584 | }else if(sc!=5){ | |
585 | if( (c1+c2)==1) {if(c1!=0) continue;} | |
586 | }else {}// do nothing for pi-k-p case | |
587 | ||
588 | ///////////////////////////////////////// | |
589 | ||
590 | ||
591 | if(PairCut){ | |
592 | ||
593 | TString *nameNorm=new TString("PairCut3_Charge1_"); | |
594 | *nameNorm += c1; | |
595 | nameNorm->Append("_Charge2_"); | |
596 | *nameNorm += c2; | |
597 | nameNorm->Append("_Charge3_"); | |
598 | *nameNorm += c3; | |
599 | nameNorm->Append("_SC_"); | |
600 | *nameNorm += sc; | |
601 | nameNorm->Append("_M_"); | |
602 | *nameNorm += Mbin; | |
603 | nameNorm->Append("_ED_"); | |
604 | *nameNorm += 0;// EDbin | |
605 | nameNorm->Append("_Term_"); | |
606 | *nameNorm += term+1; | |
607 | nameNorm->Append("_Norm"); | |
608 | NormH_pc[term] = ((TH1D*)MyList->FindObject(nameNorm->Data()))->Integral(); | |
609 | ||
610 | if(NormH_pc[term] > 0){ | |
611 | ||
612 | if(sc<=2) { | |
613 | ||
614 | TString *name_3d = new TString(name3_pc->Data()); | |
615 | name_3d->Append("_3D"); | |
616 | Three_3d[c1][c2][c3][sc][term] = (TH3D*)MyList->FindObject(name_3d->Data()); | |
617 | Three_3d[c1][c2][c3][sc][term]->Sumw2(); | |
618 | Three_3d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]); | |
619 | TString *name_Q3 = new TString(name3_pc->Data()); | |
620 | name_Q3->Append("_Q3"); | |
621 | Three_1d[c1][c2][c3][sc][term] = (TH1D*)MyList->FindObject(name_Q3->Data()); | |
622 | Three_1d[c1][c2][c3][sc][term]->Sumw2(); | |
623 | Three_1d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]); | |
624 | //cout<<"Scale factor = "<<NormH_pc[0]/NormH_pc[term]<<endl; | |
625 | ||
626 | ||
627 | } | |
628 | ||
629 | }else{ | |
630 | cout<<"normalization = 0. Skipping this SC. SC="<<sc<<" c1="<<c1<<" c2="<<c2<<" c3="<<c3<<endl; | |
631 | } | |
632 | ||
633 | ||
634 | ||
635 | }// pair cut | |
636 | }// term | |
637 | ||
638 | }// SC | |
639 | ||
640 | ||
641 | }// c3 | |
642 | ||
643 | }// c2 | |
644 | }// c1 | |
5130881c | 645 | |
0e610a08 | 646 | |
9e040a65 | 647 | cout<<"Done getting Histograms"<<endl; |
648 | ||
649 | ||
650 | ||
0e610a08 | 651 | TCanvas *can1 = new TCanvas("can1", "can1",11,53,800,800); |
9e040a65 | 652 | can1->SetHighLightColor(2); |
0e610a08 | 653 | //can1->Range(-0.7838115,-1.033258,0.7838115,1.033258); |
9e040a65 | 654 | gStyle->SetOptFit(0111); |
655 | can1->SetFillColor(10); | |
656 | can1->SetBorderMode(0); | |
657 | can1->SetBorderSize(2); | |
9e040a65 | 658 | can1->SetFrameFillColor(0); |
659 | can1->SetFrameBorderMode(0); | |
660 | can1->SetFrameBorderMode(0); | |
0e610a08 | 661 | gPad->SetRightMargin(0.025); gPad->SetLeftMargin(0.1); gPad->SetTopMargin(0.02); |
662 | gPad->SetGridx(0); | |
663 | gPad->SetGridy(0); | |
9e040a65 | 664 | TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC"); |
665 | legend1->SetBorderSize(1); | |
666 | legend1->SetTextSize(.04);// small .03; large .036 | |
667 | legend1->SetFillColor(0); | |
668 | ||
0e610a08 | 669 | TLatex *Specif_System = new TLatex(0.23, 0.9,"ALICE p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV, #LT#font[12]{N}_{ch}#GT = 9.8 #pm 0.5"); |
670 | TLatex *Specif_KT3 = new TLatex(0.23, 0.8,"0.16<#font[12]{K}_{T,3}<0.3 GeV/#font[12]{c}"); | |
671 | TLatex *Specif_kT = new TLatex(0.23, 0.8,"0.2<#font[12]{k}_{T}<0.3 GeV/#font[12]{c}"); | |
672 | TLatex *Specif_FSI = new TLatex(0.23, 0.7,"FSI uncorrected"); | |
673 | TLatex *Specif_Disclaimer = new TLatex(0.23, 0.6,"Statistical errors only"); | |
674 | Specif_FSI->SetTextFont(42); Specif_System->SetTextFont(42); Specif_KT3->SetTextFont(42); Specif_Disclaimer->SetTextFont(42); | |
675 | Specif_kT->SetTextFont(42); | |
676 | Specif_FSI->SetTextSize(0.04); Specif_System->SetTextSize(0.04); Specif_KT3->SetTextSize(0.04); Specif_Disclaimer->SetTextSize(0.04); | |
677 | Specif_kT->SetTextSize(0.04); | |
678 | Specif_FSI->SetNDC(); Specif_System->SetNDC(); Specif_KT3->SetNDC(); Specif_Disclaimer->SetNDC(); | |
679 | Specif_kT->SetNDC(); | |
680 | ||
9e040a65 | 681 | ///////////////////////// |
682 | // 2-particle legend | |
683 | // 0 = pi-pi | |
684 | // 1 = pi-k | |
685 | // 2 = pi-p | |
686 | // 3 = k-k | |
687 | // 4 = k-p | |
688 | // 5 = p-p | |
689 | ///////////////////////// | |
690 | TH1D *Two_pipi_mp = (TH1D*)Two_ex[0][1][0][0]->Clone(); | |
691 | Two_pipi_mp->Divide(Two_ex[0][1][0][1]); | |
692 | ||
693 | // Normalization range counting. Just to evaluate required normalization width in qinv. | |
694 | //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.1))<<endl; | |
695 | //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.06), Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.1))<<endl; | |
696 | //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.15), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.175))<<endl; | |
697 | ||
698 | ||
699 | ||
700 | const int SCOI_2=0; | |
701 | ||
702 | TH1D *Two_ex_clone_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone(); | |
703 | Two_ex_clone_mm->Divide(Two_ex[0][0][SCOI_2][1]); | |
704 | TH1D *Two_ex_clone_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone(); | |
705 | Two_ex_clone_mp->Divide(Two_ex[0][1][SCOI_2][1]); | |
706 | TH1D *Two_ex_clone_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone(); | |
707 | Two_ex_clone_pp->Divide(Two_ex[1][1][SCOI_2][1]); | |
708 | ||
0e610a08 | 709 | |
5130881c | 710 | Two_ex_clone_mm->GetYaxis()->SetTitle("#font[12]{C}_{2}"); |
9e040a65 | 711 | Two_ex_clone_mm->SetTitle(""); |
5130881c | 712 | Two_ex_clone_mp->GetYaxis()->SetTitle("#font[12]{C}_{2}"); |
9e040a65 | 713 | Two_ex_clone_mm->SetMarkerColor(2); |
714 | Two_ex_clone_mm->SetLineColor(2); | |
715 | Two_ex_clone_mp->SetMarkerColor(1); | |
716 | Two_ex_clone_pp->SetMarkerColor(4); | |
717 | Two_ex_clone_pp->SetLineColor(4); | |
718 | Two_ex_clone_mm->GetXaxis()->SetRangeUser(0,0.6); | |
719 | Two_ex_clone_mm->SetMinimum(0.95); | |
0e610a08 | 720 | Two_ex_clone_mm->SetMaximum(2.4); |
9e040a65 | 721 | |
0e610a08 | 722 | |
9e040a65 | 723 | if(MCcase){ |
724 | Two_ex_clone_mp->SetMarkerColor(4); | |
725 | Two_ex_clone_mp->SetLineColor(4); | |
726 | Two_ex_clone_mm->Add(Two_ex_clone_pp); | |
727 | Two_ex_clone_mm->Scale(0.5); | |
5130881c | 728 | Two_ex_clone_mm->GetYaxis()->SetTitle("#font[12]{C}_{2}^{MC}"); |
9e040a65 | 729 | Two_ex_clone_mm->GetYaxis()->SetTitleOffset(1.3); |
730 | Two_ex_clone_mm->SetMinimum(0.95); | |
731 | Two_ex_clone_mm->SetMaximum(1.3); | |
732 | Two_ex_clone_mm->DrawCopy(); | |
7c6e4906 | 733 | gStyle->SetOptFit(1111); |
0e610a08 | 734 | Two_ex_clone_mm->Fit(BkgMCFit,"","",0,1.8); |
735 | //BkgMCFit->Draw("same"); | |
9e040a65 | 736 | Two_ex_clone_mp->DrawCopy("same"); |
737 | legend1->AddEntry(Two_ex_clone_mm,"same-charge","p"); | |
738 | //legend1->AddEntry(Two_ex_clone_pp,"++","p"); | |
739 | legend1->AddEntry(Two_ex_clone_mp,"mixed-charge","p"); | |
740 | legend1->Draw("same"); | |
741 | } | |
7c6e4906 | 742 | |
0e610a08 | 743 | //for(int i=0; i<Two_ex_clone_mm->GetNbinsX(); i++){ |
744 | //cout<<Two_ex_clone_mm->GetBinError(i+1)<<", "; | |
745 | //} | |
9e040a65 | 746 | |
9e040a65 | 747 | ///////////////////////////////////////////////////// |
748 | // Global fitting C2os and C2ss | |
749 | double C2ssSys_e[BINRANGE_C2global]={0}; | |
750 | double C2osSys_e[BINRANGE_C2global]={0}; | |
751 | // | |
752 | const int npar=10;// 10 | |
753 | TMinuit MyMinuit(npar); | |
754 | MyMinuit.SetFCN(fcnC2_Global); | |
5130881c | 755 | double arglist_C2 = 0; int dummy=0; |
756 | MyMinuit.mnexcm("SET NOWarnings",&arglist_C2,1, dummy); | |
757 | arglist_C2 = -1; | |
758 | MyMinuit.mnexcm("SET PRint",&arglist_C2,1, dummy); | |
9e040a65 | 759 | double OutputPar[npar]={0}; |
760 | double OutputPar_e[npar]={0}; | |
761 | ||
762 | double par[npar]; // the start values | |
763 | double stepSize[npar]; // step sizes | |
764 | double minVal[npar]; // minimum bound on parameter | |
765 | double maxVal[npar]; // maximum bound on parameter | |
766 | string parName[npar]; | |
767 | ||
768 | TH1D *temp_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone(); | |
769 | temp_mm->Divide(Two_ex[0][0][SCOI_2][1]); | |
770 | TH1D *temp_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone(); | |
771 | temp_mp->Divide(Two_ex[0][1][SCOI_2][1]); | |
772 | TH1D *temp_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone(); | |
773 | temp_pp->Divide(Two_ex[1][1][SCOI_2][1]); | |
774 | TH1D *C2ssRaw=(TH1D*)temp_mm->Clone();// MRC uncorrected | |
775 | TH1D *C2osRaw=(TH1D*)temp_mp->Clone();// MRC uncorrected | |
776 | C2ssRaw->SetMarkerStyle(24); | |
777 | C2osRaw->SetMarkerStyle(21);//21 | |
778 | C2ssRaw->SetMarkerColor(2); | |
779 | C2osRaw->SetMarkerColor(4); | |
780 | ||
781 | for(int i=0; i<BINRANGE_C2global; i++){ | |
782 | if(i >= temp_mm->GetNbinsX()) continue;// bin limit | |
783 | ||
784 | C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.; | |
785 | double MRC_SC = 1.0; | |
786 | double MRC_MC = 1.0; | |
5130881c | 787 | if(MomResC2[0]->GetBinContent(i+1) > 0) MRC_SC = (MomResC2[0]->GetBinContent(i+1)-1)*MRCShift + 1; |
788 | if(MomResC2[1]->GetBinContent(i+1) > 0) MRC_MC = (MomResC2[1]->GetBinContent(i+1)-1)*MRCShift + 1; | |
9e040a65 | 789 | if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= MRC_SC; |
790 | C2ssFitting_e[i] = pow(MRC_SC*sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2); | |
791 | C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.); | |
792 | C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2)); | |
9e040a65 | 793 | C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]); |
794 | C2osFitting[i] = temp_mp->GetBinContent(i+1); | |
795 | if(!GeneratedSignal && !MCcase) C2osFitting[i] *= MRC_MC; | |
796 | C2osFitting_e[i] = pow(MRC_MC*temp_mp->GetBinError(i+1),2); | |
797 | C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1)); | |
798 | C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1)); | |
799 | C2osFitting_e[i] += pow((MRC_MC-1)*0.1 * temp_mp->GetBinContent(i+1),2); | |
800 | C2osFitting_e[i] = sqrt(C2osFitting_e[i]); | |
801 | // | |
802 | C2ssSys_e[i] = pow((MRC_SC-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2); | |
803 | C2ssSys_e[i] = sqrt(C2ssSys_e[i]); | |
804 | C2osSys_e[i] = pow((MRC_MC-1)*0.1 * temp_mp->GetBinContent(i+1),2); | |
805 | C2osSys_e[i] = sqrt(C2osSys_e[i]); | |
806 | // | |
0e610a08 | 807 | if(UseC2Bkg) C2ssFitting[i] /= BkgMCFit->Eval(BinCenters[i]); |
808 | // | |
809 | // method with undilution for plotting purposes | |
810 | //C2ssFitting[i] = (C2ssFitting[i] - (1-lambda_PM)) / (CoulCorr2(+1, BinCenters[i]) * lambda_PM); | |
7c6e4906 | 811 | // |
9e040a65 | 812 | K2SS[i] = CoulCorr2(+1, BinCenters[i]); |
813 | K2OS[i] = CoulCorr2(-1, BinCenters[i]); | |
814 | //K2SS[i] = 1; | |
815 | //K2OS[i] = 1; | |
816 | // | |
9e040a65 | 817 | K2SS_e[i] = 0.0; |
818 | K2OS_e[i] = 0.0; | |
819 | } | |
0e610a08 | 820 | |
9e040a65 | 821 | |
822 | ||
823 | ||
7c6e4906 | 824 | |
9e040a65 | 825 | C2ssFitting[0]=-1000; C2osFitting[0]=-1000; |
826 | C2ssFitting_e[0]=1000; C2osFitting_e[0]=1000; | |
827 | K2SS_e[0]=1000; K2OS_e[0]=1000; | |
9e040a65 | 828 | int ierflg=0; |
829 | double arglist[10]; | |
830 | //arglist[0]=2;// improve Minimization Strategy (1 is default) | |
831 | //MyMinuit.mnexcm("SET STR",arglist,1,ierflg); | |
832 | //arglist[0] = 0; | |
833 | //MyMinuit.mnexcm("SCAN", arglist,1,ierflg); | |
0e610a08 | 834 | arglist[0] = 5000;// 5000 |
9e040a65 | 835 | MyMinuit.mnexcm("MIGRAD", arglist ,1,ierflg); |
9e040a65 | 836 | |
0e610a08 | 837 | TF1 *fitC2ss_Base = new TF1("fitC2ss_Base",C2ssFitFunction, 0.005,2, npar);//0.2 |
838 | TF1 *fitC2ss_Expan = new TF1("fitC2ss_Expan",C2ssFitFunction, 0.005,2, npar);//0.2 | |
839 | fitC2ss_Base->SetLineStyle(2); | |
b2c0d4ba | 840 | TH1D *fitC2ss_h = new TH1D("fitC2ss_h","",Two_ex_clone_mm->GetNbinsX(),Two_ex_clone_mm->GetXaxis()->GetBinLowEdge(1), Two_ex_clone_mm->GetXaxis()->GetBinUpEdge(Two_ex_clone_mm->GetNbinsX())); |
841 | ||
842 | for(int ft=0; ft<2; ft++){// Gaussian or EW | |
0e610a08 | 843 | if(ft==1 && !IncludeExpansion) continue; |
844 | par[0] = 1.0; par[1]=0.7; par[2]=0.5; par[3]=7.2; par[4] = .1; par[5] = .0; par[6] = .0; par[7] = 0.; par[8] = 0.; par[9] = 0.;// was par[1]=0.5 | |
b2c0d4ba | 845 | stepSize[0] = 0.01; stepSize[1] = 0.01; stepSize[2] = 0.02; stepSize[3] = 0.2; stepSize[4] = 0.01; stepSize[5] = 0.001; stepSize[6] = 0.001; stepSize[7] = 0.001; stepSize[8]=0.001; stepSize[9]=0.01; |
0e610a08 | 846 | minVal[0] = 0.955; minVal[1] = 0.5; minVal[2] = 0.; minVal[3] = 0.1; minVal[4] = 0.001; minVal[5] = -10.; minVal[6] = -10.; minVal[7] = -10.; minVal[8]=-10; minVal[9] = 0.995;// minVal[1] was 0.2 |
847 | maxVal[0] = 1.1; maxVal[1] = 1.3; maxVal[2] = 0.99; maxVal[3] = 15.; maxVal[4] = 2.; maxVal[5] = 10.; maxVal[6] = 10.; maxVal[7] = 10.; maxVal[8]=10.; maxVal[9]=1.1;// maxVal[1] was 1.0 | |
848 | if(Gaussian==kFALSE) maxVal[1]=2.0; | |
849 | ||
b2c0d4ba | 850 | parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh"; |
0e610a08 | 851 | parName[5] = "coeff_3"; parName[6] = "coeff_4"; parName[7] = "coeff_5"; parName[8]="coeff_6"; parName[9]="Norm_2"; |
b2c0d4ba | 852 | |
853 | for (int i=0; i<npar; i++){ | |
854 | MyMinuit.DefineParameter(i, parName[i].c_str(), par[i], stepSize[i], minVal[i], maxVal[i]); | |
855 | } | |
856 | ||
857 | //MyMinuit.DefineParameter(1, parName[1].c_str(), 0.7, stepSize[1], minVal[1], maxVal[1]); MyMinuit.FixParameter(1);// lambda | |
858 | //MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N | |
859 | MyMinuit.DefineParameter(2, parName[2].c_str(), 0., stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G | |
860 | MyMinuit.DefineParameter(4, parName[4].c_str(), 0., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh | |
861 | MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5 | |
862 | MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6 | |
863 | ||
864 | // | |
865 | if(ft==0){ | |
866 | MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 | |
867 | MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 | |
0e610a08 | 868 | }else{// IncludeExpansion |
869 | if(FixExpansionAvg){ | |
b2c0d4ba | 870 | MyMinuit.DefineParameter(5, parName[5].c_str(), kappa3, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 |
871 | MyMinuit.DefineParameter(6, parName[6].c_str(), kappa4, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 | |
0e610a08 | 872 | |
b2c0d4ba | 873 | }else{ |
874 | Int_t err=0; | |
875 | Double_t tmp[1]; | |
876 | tmp[0] = 5+1; | |
877 | MyMinuit.mnexcm( "RELease", tmp, 1, err );// kappa3 | |
878 | tmp[0] = 6+1; | |
879 | MyMinuit.mnexcm( "RELease", tmp, 1, err );// kappa4 | |
0e610a08 | 880 | // |
b2c0d4ba | 881 | } |
882 | } | |
883 | ||
884 | if(IncludeSS==kFALSE){ | |
885 | MyMinuit.DefineParameter(3, parName[3].c_str(), 9.1, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch | |
886 | MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N | |
887 | MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 | |
888 | MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 | |
889 | } | |
890 | if(IncludeOS==kFALSE){ | |
891 | MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2 | |
892 | } | |
893 | ||
894 | ||
895 | // Do the minimization! | |
896 | if(!MCcase){ | |
897 | cout<<"Start C2 Global fit"<<endl; | |
898 | MyMinuit.Migrad();// generally the best minimization algorithm | |
899 | cout<<"End C2 Global fit"<<endl; | |
900 | } | |
901 | for (int i=0; i<npar; i++){ | |
902 | MyMinuit.GetParameter(i,OutputPar[i],OutputPar_e[i]); | |
903 | } | |
904 | MyMinuit.mnexcm("SHOw PARameters", &arglist_C2, 1, ierflg); | |
905 | cout<<"C2 fit: Chi2/NDF = "<<Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars())<<" Chi^2="<<Chi2_C2global<<endl; | |
906 | ||
0e610a08 | 907 | if(ft==0){// Gaussian or Exponential |
908 | for(int i=0; i<npar; i++) { | |
909 | fitC2ss_Base->FixParameter(i,OutputPar[i]); | |
910 | fitC2ss_Base->SetParError(i,OutputPar_e[i]); | |
911 | } | |
912 | ||
913 | }else{// Edgeworth or Laguerre | |
b2c0d4ba | 914 | for(int i=0; i<npar; i++) { |
0e610a08 | 915 | fitC2ss_Expan->FixParameter(i,OutputPar[i]); |
916 | fitC2ss_Expan->SetParError(i,OutputPar_e[i]); | |
b2c0d4ba | 917 | } |
918 | for(int bin=1; bin<=Two_ex_clone_mm->GetNbinsX(); bin++){ | |
919 | double qinv = Two_ex_clone_mm->GetXaxis()->GetBinCenter(bin); | |
0e610a08 | 920 | if(!MCcase) fitC2ss_h->SetBinContent(bin, fitC2ss_Expan->Eval(qinv));// uncomment to show expansion |
b2c0d4ba | 921 | } |
0e610a08 | 922 | |
923 | double lambdaStar=0, lambdaStar_e=0; | |
924 | if(Gaussian){ | |
925 | lambdaStar = OutputPar[1]*pow(1 + OutputPar[6]/8.,2); | |
926 | lambdaStar_e = pow(OutputPar_e[1]*(1 + OutputPar[6]/8.),2); | |
927 | lambdaStar_e = sqrt(lambdaStar_e); | |
928 | }else{ | |
929 | lambdaStar = OutputPar[1]*pow(1 - OutputPar[5] + OutputPar[6],2); | |
930 | lambdaStar_e = pow(OutputPar_e[1]*(1 - OutputPar[5] + OutputPar[6]),2); | |
931 | lambdaStar_e = sqrt(lambdaStar_e); | |
b2c0d4ba | 932 | } |
0e610a08 | 933 | cout<<"lambda* = "<<lambdaStar<<" +- "<<lambdaStar_e<<endl; |
b2c0d4ba | 934 | } |
0e610a08 | 935 | //if(ft==0) {MainFitParams[0]=OutputPar[3]; MainFitParams[2]=OutputPar[1];} |
936 | //else {MainFitParams[4]=OutputPar[3]; MainFitParams[6]=OutputPar[1];} | |
b2c0d4ba | 937 | }// ft |
9e040a65 | 938 | |
7c6e4906 | 939 | |
9e040a65 | 940 | TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone(); |
941 | TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone(); | |
942 | TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone(); | |
943 | TH1D *C2_os_Momsys=(TH1D*)Two_ex_clone_mp->Clone(); | |
944 | TH1D *C2_ss_Ksys=(TH1D*)Two_ex_clone_mm->Clone(); | |
945 | TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone(); | |
946 | TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone(); | |
947 | TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone(); | |
948 | ||
949 | for(int i=0; i<BINRANGE_C2global; i++){ | |
950 | C2_ss->SetBinContent(i+1, C2ssFitting[i]); | |
951 | C2_os->SetBinContent(i+1, C2osFitting[i]); | |
952 | C2_ss_Momsys->SetBinContent(i+1, C2ssFitting[i]); | |
953 | C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]); | |
954 | C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]); | |
955 | C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]); | |
956 | double Toterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2))); | |
957 | double Toterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2))); | |
958 | C2_ss->SetBinError(i+1, Toterror_ss); | |
959 | C2_os->SetBinError(i+1, Toterror_os); | |
960 | C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]); | |
961 | C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]); | |
962 | C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]); | |
963 | C2_os_Ksys->SetBinError(i+1, OutputPar[1]*K2OS_e[i]); | |
964 | // | |
965 | K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0); | |
966 | K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0); | |
967 | } | |
968 | ||
969 | C2_ss_Momsys->SetMarkerSize(0); | |
970 | C2_ss_Momsys->SetFillStyle(1000); | |
971 | C2_ss_Momsys->SetFillColor(kRed-10); | |
972 | C2_os_Momsys->SetMarkerSize(0); | |
973 | C2_os_Momsys->SetFillStyle(1000); | |
974 | C2_os_Momsys->SetFillColor(kBlue-10); | |
975 | C2_ss_Ksys->SetMarkerSize(0); | |
976 | C2_ss_Ksys->SetFillStyle(3004); | |
977 | C2_ss_Ksys->SetFillColor(kRed); | |
978 | C2_os_Ksys->SetMarkerSize(0); | |
979 | C2_os_Ksys->SetFillStyle(3004); | |
980 | C2_os_Ksys->SetFillColor(kBlue); | |
981 | ||
982 | ||
983 | C2_ss->GetXaxis()->SetRangeUser(0,Q2Limit);// 0,0.098 | |
0e610a08 | 984 | C2_ss->GetYaxis()->SetRangeUser(0.95,2.4);//0.98,1.3 |
9e040a65 | 985 | C2_ss->GetXaxis()->SetTitleOffset(.8); |
0e610a08 | 986 | C2_ss->GetYaxis()->SetTitleOffset(.85); |
987 | C2_ss->GetXaxis()->SetTitle("#font[12]{q} (GeV/c)"); | |
988 | C2_ss->GetYaxis()->SetTitle("#font[12]{C}_{2}^{#pm#pm}"); | |
9e040a65 | 989 | C2_ss->GetXaxis()->SetTitleSize(0.05); |
990 | C2_ss->GetYaxis()->SetTitleSize(0.05); | |
991 | C2_os->GetXaxis()->SetRangeUser(0,0.6); | |
992 | C2_os->GetYaxis()->SetRangeUser(0.98,1.3); | |
993 | C2_os->GetXaxis()->SetTitleOffset(.8); | |
994 | C2_os->GetYaxis()->SetTitleOffset(.8); | |
0e610a08 | 995 | C2_os->GetXaxis()->SetTitle("#font[12]{q} (GeV/c)"); |
9e040a65 | 996 | C2_os->GetXaxis()->SetTitleSize(0.05); |
997 | C2_os->GetYaxis()->SetTitleSize(0.05); | |
998 | ||
0e610a08 | 999 | //C2_ss->SetMarkerSize(1.5); |
1000 | //C2_os->SetMarkerSize(1.5); | |
9e040a65 | 1001 | C2_os->SetMarkerStyle(25); |
1002 | C2_os->SetMarkerColor(4); | |
1003 | ||
1004 | ||
9e040a65 | 1005 | fitC2ss_h->SetLineWidth(2); |
1006 | fitC2ss_h->SetLineColor(2); | |
1007 | ||
5130881c | 1008 | |
9e040a65 | 1009 | if(!MCcase){ |
1010 | ||
1011 | C2_ss->DrawCopy(); | |
1012 | legend1->AddEntry(C2_ss,"same-charge","p"); | |
0e610a08 | 1013 | //C2_os->DrawCopy("same"); |
9e040a65 | 1014 | legend1->AddEntry(C2_os,"mixed-charge","p"); |
5130881c | 1015 | |
0e610a08 | 1016 | Specif_System->Draw("same"); |
1017 | Specif_kT->Draw("same"); | |
1018 | Specif_FSI->Draw("same"); | |
1019 | Specif_Disclaimer->Draw("same"); | |
1020 | ||
1021 | TLatex *Stats1 = new TLatex(.5,.55, "#lambda = 0.41 #pm 0.004"); | |
5130881c | 1022 | Stats1->SetNDC(); |
1023 | Stats1->SetTextSize(0.06); | |
7c6e4906 | 1024 | //Stats1->Draw("same"); |
0e610a08 | 1025 | TLatex *Stats2 = new TLatex(.5,.45, "R_{inv} = 1.59 #pm 0.01 fm"); |
5130881c | 1026 | Stats2->SetNDC(); |
1027 | Stats2->SetTextSize(0.06); | |
7c6e4906 | 1028 | //Stats2->Draw("same"); |
9e040a65 | 1029 | /*TF1 *BkgFitC2 = new TF1("BkgFitC2","1+[0]*exp(-pow([1]*x/0.19733,2))",0,1); |
1030 | BkgFitC2->SetParameter(0,0.08); | |
1031 | BkgFitC2->SetParameter(1,0.5); | |
1032 | BkgFitC2->SetLineColor(1); | |
1033 | C2_ss->Fit(BkgFitC2,"IME","",0.2,0.8); | |
1034 | BkgFitC2->Draw("same");*/ | |
0e610a08 | 1035 | |
1036 | Unity->Draw("same"); | |
1037 | //fitC2ss_h->Draw("C same"); | |
1038 | //fitC2ss_Base->Draw("C same"); | |
1039 | //fitC2ss_Expan->Draw("C same"); | |
1040 | //legend1->AddEntry(fitC2ss_Base,"Gaussian","l"); | |
1041 | //legend1->AddEntry(fitC2ss_Expan,"Edgeworth","l"); | |
1042 | //legend1->Draw("same"); | |
9e040a65 | 1043 | } |
1044 | ||
1045 | ||
9e040a65 | 1046 | |
1047 | cout<<"============================================"<<endl; | |
1048 | cout<<"Start 3-pion section"<<endl; | |
1049 | ||
0e610a08 | 1050 | TCanvas *can2 = new TCanvas("can2", "can2",800,0,800,800);//800,0,600,900 |
9e040a65 | 1051 | can2->SetHighLightColor(2); |
1052 | gStyle->SetOptFit(0111); | |
1053 | can2->SetFillColor(10); | |
1054 | can2->SetBorderMode(0); | |
1055 | can2->SetBorderSize(2); | |
1056 | can2->SetGridx(); | |
1057 | can2->SetGridy(); | |
1058 | can2->SetFrameFillColor(0); | |
1059 | can2->SetFrameBorderMode(0); | |
1060 | can2->SetFrameBorderMode(0); | |
1061 | can2->cd(); | |
9e040a65 | 1062 | gPad->SetRightMargin(0.02); |
1063 | ||
1064 | TLegend *legend2 = new TLegend(.58,.55, .97,.65,NULL,"brNDC"); | |
1065 | legend2->SetBorderSize(1); | |
1066 | legend2->SetTextSize(.03);// small .03; large .06 | |
1067 | legend2->SetFillColor(0); | |
1068 | ||
1069 | ///////////////////////////////////////////// | |
1070 | TH3D *C3QS_3d = new TH3D("C3QS_3d","",BINRANGE_3,0,.4, BINRANGE_3,0,.4, BINRANGE_3,0,.4); | |
1071 | TH3D *Combinatorics_3d = new TH3D("Combinatorics_3d","",BINRANGE_3,0,.4, BINRANGE_3,0,.4, BINRANGE_3,0,.4); | |
9e040a65 | 1072 | |
1073 | // | |
1074 | const float Q3HistoLimit=0.5; | |
1075 | const int Q3BINS = int((Q3HistoLimit+0.0001)/0.01); | |
1076 | TH1D *num_withRS = new TH1D("num_withRS","",Q3BINS,0,Q3HistoLimit); | |
1077 | TH1D *num_withGRS = new TH1D("num_withGRS","",Q3BINS,0,Q3HistoLimit); | |
1078 | TH1D *num_withTM = new TH1D("num_withTM","",Q3BINS,0,Q3HistoLimit); | |
1079 | TH1D *Cterm1 = new TH1D("Cterm1","",Q3BINS,0,Q3HistoLimit); | |
1080 | TH1D *Cterm1_noMRC = new TH1D("Cterm1_noMRC","",Q3BINS,0,Q3HistoLimit); | |
1081 | TH1D *Cterm2 = new TH1D("Cterm2","",Q3BINS,0,Q3HistoLimit); | |
1082 | TH1D *Cterm3 = new TH1D("Cterm3","",Q3BINS,0,Q3HistoLimit); | |
1083 | TH1D *Cterm4 = new TH1D("Cterm4","",Q3BINS,0,Q3HistoLimit); | |
1084 | TH1D *num_QS = new TH1D("num_QS","",Q3BINS,0,Q3HistoLimit); | |
1085 | TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,Q3HistoLimit); | |
1086 | TH1D *Combinatorics_1d_noMRC = new TH1D("Combinatorics_1d_noMRC","",Q3BINS,0,Q3HistoLimit); | |
1087 | TH1D *Coul_Riverside = new TH1D("Coul_Riverside","",Q3BINS,0,Q3HistoLimit); | |
1088 | TH1D *Coul_GRiverside = new TH1D("Coul_GRiverside","",Q3BINS,0,Q3HistoLimit); | |
1089 | TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,Q3HistoLimit); | |
0e610a08 | 1090 | |
9e040a65 | 1091 | if(SameCharge) {Cterm1_noMRC->Sumw2(); Combinatorics_1d_noMRC->Sumw2();} |
1092 | // | |
1093 | double num_QS_e[Q3BINS]; | |
1094 | double c3_e[Q3BINS]; | |
b2c0d4ba | 1095 | |
9e040a65 | 1096 | for(int ii=0; ii<Q3BINS; ii++){ |
1097 | num_QS_e[ii]=0.; | |
1098 | c3_e[ii]=0.; | |
9e040a65 | 1099 | } |
1100 | ||
1101 | // CB=Charge Bin. 0 means pi-, 1 means pi+ | |
1102 | int CB1=0, CB2=0, CB3=0; | |
1103 | int CP12=1, CP13=1, CP23=1; | |
1104 | int CB1_2=0, CB2_2=0, CB3_2=0; | |
1105 | if(CHARGE==-1) { | |
1106 | if(SameCharge) {CB1=0; CB2=0; CB3=0;} | |
1107 | else {CB1=0; CB2=0; CB3=1; CP12=+1, CP13=-1, CP23=-1;} | |
1108 | }else { | |
1109 | if(SameCharge) {CB1=1; CB2=1; CB3=1;} | |
1110 | else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;} | |
1111 | } | |
1112 | if(AddCC){ | |
1113 | if(CHARGE==-1) { | |
1114 | if(SameCharge) {CB1_2=1; CB2_2=1; CB3_2=1;} | |
1115 | else {CB1_2=0; CB2_2=1; CB3_2=1;}; | |
1116 | }else { | |
1117 | if(SameCharge) {CB1_2=0; CB2_2=0; CB3_2=0;} | |
1118 | else {CB1_2=0; CB2_2=0; CB3_2=1;} | |
1119 | } | |
1120 | } | |
1121 | ||
1122 | ||
1123 | // SC = species combination. Always 0 meaning pi-pi-pi. | |
1124 | int SCBin=0; | |
1125 | // | |
7c6e4906 | 1126 | |
9e040a65 | 1127 | TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",Q3BINS,0,Q3HistoLimit); |
1128 | TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",Q3BINS,0,Q3HistoLimit); | |
1129 | // | |
1130 | double value_num; | |
1131 | double value_num_e; | |
1132 | double N3_QS; | |
1133 | double N3_QS_e; | |
9e040a65 | 1134 | |
1135 | for(int i=2; i<=BINLIMIT_3; i++){// bin number | |
9e040a65 | 1136 | double Q12 = BinCenters[i-1];// true center |
1137 | //int Q12bin = int(Q12/0.01)+1; | |
1138 | // | |
1139 | for(int j=2; j<=BINLIMIT_3; j++){// bin number | |
9e040a65 | 1140 | double Q13 = BinCenters[j-1];// true center |
1141 | //int Q13bin = int(Q13/0.01)+1; | |
1142 | // | |
1143 | for(int k=2; k<=BINLIMIT_3; k++){// bin number | |
9e040a65 | 1144 | double Q23 = BinCenters[k-1];// true center |
1145 | //int Q23bin = int(Q23/0.01)+1; | |
1146 | // | |
b2c0d4ba | 1147 | |
9e040a65 | 1148 | double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2)); |
1149 | int Q3bin = Cterm1->GetXaxis()->FindBin(Q3); | |
0e610a08 | 1150 | |
9e040a65 | 1151 | if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible |
1152 | if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible | |
1153 | ||
1154 | ||
1155 | // | |
1156 | double G_12 = Gamov(CP12, Q12);// point-source Coulomb correlation | |
1157 | double G_13 = Gamov(CP13, Q13); | |
1158 | double G_23 = Gamov(CP23, Q23); | |
1159 | double K2_12 = CoulCorr2(CP12, Q12);// finite-source Coulomb+Strong correlation from Therminator. | |
1160 | double K2_13 = CoulCorr2(CP13, Q13); | |
1161 | double K2_23 = CoulCorr2(CP23, Q23); | |
1162 | double K3 = 1.;// 3-body Coulomb+Strong correlation | |
1163 | if(SameCharge || CHARGE==-1) K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23); | |
1164 | else K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13); | |
1165 | ||
1166 | if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;} | |
1167 | if(K3==0) continue; | |
1168 | ||
1169 | double TERM1=Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k);// N3 (3-pion yield per q12,q13,q23 cell). 3-pions from same-event | |
1170 | double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event | |
1171 | double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event | |
1172 | double TERM4=Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event | |
1173 | double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics) | |
1174 | if(AddCC){ | |
b2c0d4ba | 1175 | if(CHARGE==-1){// base is (SC,MC,MC) |
1176 | TERM1 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][0]->GetBinContent(j,k,i); | |
1177 | TERM2 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][1]->GetBinContent(j,k,i); | |
1178 | TERM3 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][2]->GetBinContent(j,k,i); | |
1179 | TERM4 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][3]->GetBinContent(j,k,i); | |
1180 | TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(j,k,i); | |
1181 | }else {// base is (MC,MC,SC) | |
1182 | TERM1 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][0]->GetBinContent(k,i,j); | |
1183 | TERM2 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][1]->GetBinContent(k,i,j); | |
1184 | TERM3 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][2]->GetBinContent(k,i,j); | |
1185 | TERM4 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][3]->GetBinContent(k,i,j); | |
1186 | TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(k,i,j); | |
1187 | } | |
9e040a65 | 1188 | } |
1189 | ||
1190 | if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue; | |
1191 | if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue; | |
1192 | ||
5be394f6 | 1193 | double muonCorr_C3=1.0, muonCorr_C2_12=1.0, muonCorr_C2_13=1.0, muonCorr_C2_23=1.0; |
0e610a08 | 1194 | if(MuonCorrection && !MCcase){ |
5be394f6 | 1195 | if(SameCharge){ |
1196 | muonCorr_C3 = C3muonCorrection->GetBinContent(C3muonCorrection->GetXaxis()->FindBin(Q3)); | |
1197 | muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q12)); | |
1198 | muonCorr_C2_13 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q13)); | |
1199 | muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q23)); | |
1200 | }else{ | |
1201 | muonCorr_C3 = C3muonCorrection->GetBinContent(C3muonCorrection->GetXaxis()->FindBin(Q3)); | |
1202 | if(CHARGE==-1){// pi-pi-pi+ | |
1203 | muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q12)); | |
1204 | muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q13)); | |
1205 | muonCorr_C2_23 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q23)); | |
1206 | }else{// pi-pi+pi+ | |
1207 | muonCorr_C2_12 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q12)); | |
1208 | muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q13)); | |
1209 | muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q23)); | |
1210 | } | |
1211 | } | |
1212 | } | |
9e040a65 | 1213 | |
9e040a65 | 1214 | // apply momentum resolution correction |
1215 | if(!MCcase && !GeneratedSignal){ | |
b2c0d4ba | 1216 | int MRC_Q3bin = int(Q3/0.01) + 1; |
9e040a65 | 1217 | if(SameCharge){ |
b2c0d4ba | 1218 | TERM1 *= (MomRes1d[0][0]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; |
1219 | TERM2 *= (MomRes1d[0][1]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1220 | TERM3 *= (MomRes1d[0][2]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1221 | TERM4 *= (MomRes1d[0][3]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1222 | TERM5 *= (MomRes1d[0][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
9e040a65 | 1223 | }else{ |
5130881c | 1224 | // removed due to over-correction of 1-D MRC approximation |
b2c0d4ba | 1225 | TERM1 *= (MomRes1d[1][0]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; |
1226 | TERM2 *= (MomRes1d[1][1]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1227 | TERM3 *= (MomRes1d[1][2]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1228 | TERM4 *= (MomRes1d[1][3]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1229 | TERM5 *= (MomRes1d[1][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
9e040a65 | 1230 | } |
1231 | } | |
1232 | // | |
1233 | ||
0e610a08 | 1234 | TwoFrac=lambda_PM;// 0.7 |
9e040a65 | 1235 | OneFrac=pow(TwoFrac,.5); ThreeFrac=pow(TwoFrac,1.5); |
9e040a65 | 1236 | |
1237 | ||
1238 | // Purify. Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations) | |
1239 | N3_QS = TERM1; | |
1240 | N3_QS -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5; | |
1241 | N3_QS -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5); | |
1242 | N3_QS /= ThreeFrac; | |
1243 | N3_QS /= K3; | |
5be394f6 | 1244 | N3_QS *= muonCorr_C3; |
1245 | ||
9e040a65 | 1246 | num_QS->Fill(Q3, N3_QS); |
1247 | ||
1248 | // Isolate 3-pion cumulant | |
1249 | value_num = N3_QS; | |
5be394f6 | 1250 | value_num -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac * muonCorr_C2_12; |
1251 | value_num -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac * muonCorr_C2_13; | |
1252 | value_num -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac * muonCorr_C2_23; | |
9e040a65 | 1253 | value_num += 2*TERM5; |
1254 | ||
1255 | ||
1256 | ||
1257 | ||
1258 | // errors | |
1259 | N3_QS_e = TERM1; | |
1260 | N3_QS_e += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2); | |
1261 | N3_QS_e += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2)); | |
1262 | N3_QS_e /= pow(ThreeFrac,2); | |
1263 | N3_QS_e /= pow(K3,2); | |
1264 | num_QS_e[Q3bin-1] += N3_QS_e; | |
1265 | // errors | |
1266 | value_num_e = N3_QS_e; | |
1267 | value_num_e += (pow(1/K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(TERM5),2)); | |
1268 | value_num_e += (pow(1/K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(TERM5),2)); | |
1269 | value_num_e += (pow(1/K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(TERM5),2)); | |
1270 | value_num_e += pow(2*sqrt(TERM5),2); | |
b2c0d4ba | 1271 | |
1272 | c3_e[Q3bin-1] += value_num_e + TERM5;// add baseline stat error | |
1273 | ||
9e040a65 | 1274 | |
1275 | // Fill histograms | |
1276 | c3_hist->Fill(Q3, value_num + TERM5);// for cumulant correlation function | |
9e040a65 | 1277 | C3QS_3d->SetBinContent(i,j,k, N3_QS); |
1278 | C3QS_3d->SetBinError(i,j,k, N3_QS_e); | |
1279 | // | |
1280 | Coul_Riverside->Fill(Q3, G_12*G_13*G_23*TERM5); | |
1281 | Coul_GRiverside->Fill(Q3, TERM5*CoulCorrGRS(SameCharge, Q12, Q13, Q23)); | |
1282 | num_withGRS->Fill(Q3, N3_QS); | |
1283 | Cterm1->Fill(Q3, TERM1); | |
1284 | Cterm2->Fill(Q3, TERM2); | |
1285 | Cterm3->Fill(Q3, TERM3); | |
1286 | Cterm4->Fill(Q3, TERM4); | |
1287 | Combinatorics_1d->Fill(Q3, TERM5); | |
1288 | Combinatorics_3d->Fill(Q12,Q13,Q23, TERM5); | |
1289 | ||
1290 | double Q3_m = sqrt(pow((i-0.5)*(0.005),2) + pow((j-0.5)*(0.005),2) + pow((k-0.5)*(0.005),2)); | |
1291 | Cterm1_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k)); | |
1292 | Combinatorics_1d_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k)); | |
1293 | ||
1294 | // | |
1295 | A_3[i-1][j-1][k-1] = value_num + TERM5; | |
1296 | B_3[i-1][j-1][k-1] = TERM5; | |
b2c0d4ba | 1297 | A_3_e[i-1][j-1][k-1] = sqrt(value_num_e + TERM5); |
9e040a65 | 1298 | |
1299 | /////////////////////////////////////////////////////////// | |
9e040a65 | 1300 | |
1301 | } | |
1302 | } | |
1303 | } | |
1304 | ||
1305 | ||
5130881c | 1306 | |
9e040a65 | 1307 | //////////////////////////// |
1308 | ||
1309 | // Intermediate steps | |
1310 | num_withRS->Sumw2(); | |
1311 | num_withGRS->Sumw2(); | |
1312 | num_withTM->Sumw2(); | |
1313 | Cterm1->Sumw2(); | |
1314 | Cterm2->Sumw2(); | |
1315 | Cterm3->Sumw2(); | |
1316 | Cterm4->Sumw2(); | |
1317 | Combinatorics_1d->Sumw2(); | |
1318 | Combinatorics_3d->Sumw2(); | |
1319 | for(int i=0; i<Q3BINS; i++) {c3_hist->SetBinError(i+1, sqrt(c3_e[i])); num_QS->SetBinError(i+1, sqrt(num_QS_e[i]));} | |
b2c0d4ba | 1320 | |
9e040a65 | 1321 | |
1322 | num_withRS->Divide(Combinatorics_1d); | |
1323 | num_withGRS->Divide(Combinatorics_1d); | |
1324 | num_withTM->Divide(Combinatorics_1d); | |
1325 | Cterm1->Divide(Combinatorics_1d); | |
1326 | Cterm1_noMRC->Divide(Combinatorics_1d_noMRC); | |
1327 | Cterm2->Divide(Combinatorics_1d); | |
1328 | Cterm3->Divide(Combinatorics_1d); | |
1329 | Cterm4->Divide(Combinatorics_1d); | |
1330 | c3_hist->Divide(Combinatorics_1d); | |
1331 | num_QS->Divide(Combinatorics_1d); | |
0e610a08 | 1332 | |
1333 | ||
1334 | ||
9e040a65 | 1335 | /////////////////////////////////////////////////////////////////////////////////////////////////// |
1336 | ||
1337 | ||
1338 | Cterm1->SetMarkerStyle(20); | |
1339 | Cterm1->SetMarkerColor(4); | |
1340 | Cterm1->SetLineColor(4); | |
5130881c | 1341 | Cterm1->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})"); |
1342 | Cterm1->GetYaxis()->SetTitle("#font[12]{C}_{3}"); | |
9e040a65 | 1343 | //Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}"); |
1344 | Cterm1->SetMinimum(0.97); | |
1345 | Cterm1->SetMaximum(5.7);// 6.1 for same-charge | |
1346 | Cterm1->GetXaxis()->SetRangeUser(0,Q3Limit); | |
1347 | Cterm1->GetYaxis()->SetTitleOffset(1.4); | |
1348 | Cterm1->Draw(); | |
5130881c | 1349 | legend2->AddEntry(Cterm1,"#font[12]{C}_{3} raw","p"); |
9e040a65 | 1350 | // |
1351 | Cterm2->SetMarkerStyle(20); | |
1352 | Cterm2->SetMarkerColor(7); | |
1353 | ||
1354 | Cterm3->SetMarkerStyle(25); | |
1355 | Cterm3->SetMarkerColor(8); | |
1356 | Cterm4->SetMarkerStyle(26); | |
1357 | Cterm4->SetMarkerColor(9); | |
1358 | //Cterm2->Draw("same"); | |
1359 | //Cterm3->Draw("same"); | |
1360 | //Cterm4->Draw("same"); | |
1361 | //legend2->AddEntry(Cterm1,"++-","p"); | |
1362 | ||
1363 | ||
1364 | if(MCcase){ | |
1365 | double C3points_HIJING_mmm[10]={0, 0.961834, 1.01827, 0.999387, 1.00202, 1.00081, 1.00082, 1.00037, 0.999074, 0.999099}; | |
1366 | double C3points_HIJING_mmm_e[10]={0, 0.0833895, 0.015158, 0.0047978, 0.00235067, 0.00138155, 0.00087485, 0.000612203, 0.000450162, 0.000346943}; | |
1367 | double C3points_HIJING_ppp[10]={0, 1.13015, 1.00623, 1.00536, 1.00293, 0.999964, 1.00116, 1.0007, 1.00037, 1.00105}; | |
1368 | double C3points_HIJING_ppp_e[10]={0, 0.0977058, 0.0150694, 0.0048196, 0.00235518, 0.00138376, 0.000877562, 0.000614069, 0.000452051, 0.00034856}; | |
1369 | TH1D *C3_HIJING_mmm=(TH1D*)Cterm1->Clone(); | |
1370 | TH1D *C3_HIJING_ppp=(TH1D*)Cterm1->Clone(); | |
1371 | for(int i=0; i<10; i++){ | |
1372 | C3_HIJING_mmm->SetBinContent(i+1, C3points_HIJING_mmm[i]); | |
1373 | C3_HIJING_mmm->SetBinError(i+1, C3points_HIJING_mmm_e[i]); | |
1374 | C3_HIJING_ppp->SetBinContent(i+1, C3points_HIJING_ppp[i]); | |
1375 | C3_HIJING_ppp->SetBinError(i+1, C3points_HIJING_ppp_e[i]); | |
1376 | } | |
1377 | C3_HIJING_mmm->SetMarkerColor(2); | |
1378 | C3_HIJING_mmm->SetLineColor(2); | |
1379 | C3_HIJING_ppp->SetMarkerColor(4); | |
1380 | C3_HIJING_ppp->SetLineColor(4); | |
1381 | //C3_HIJING_mmm->Draw("same"); | |
1382 | //C3_HIJING_ppp->Draw("same"); | |
1383 | //legend2->AddEntry(C3_HIJING_mmm,"---","p"); | |
1384 | //legend2->AddEntry(C3_HIJING_ppp,"+++","p"); | |
1385 | } | |
1386 | ||
1387 | num_QS->SetMarkerStyle(24); | |
1388 | num_QS->SetMarkerColor(4); | |
1389 | num_QS->SetLineColor(4); | |
5130881c | 1390 | num_QS->GetXaxis()->SetTitle("#font[12]{Q}_{3}"); |
1391 | num_QS->GetYaxis()->SetTitle("#font[12]{C}_{3}^{QS}"); | |
9e040a65 | 1392 | num_QS->GetXaxis()->SetRangeUser(0,.12); |
1393 | num_QS->SetMaximum(6); | |
1394 | //num_QS->Draw("same"); | |
1395 | //legend2->AddEntry(num_QS,"C_{3}^{QS}","p"); | |
1396 | ||
5130881c | 1397 | c3_hist->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})"); |
1398 | c3_hist->GetYaxis()->SetTitle("#font[12]{C}_{3} or #font[12]{#bf{c}}_{3}"); | |
9e040a65 | 1399 | c3_hist->GetYaxis()->SetTitleOffset(1.4); |
1400 | c3_hist->GetXaxis()->SetRangeUser(0,Q3Limit); | |
1401 | c3_hist->GetYaxis()->SetRangeUser(0.9,4); | |
1402 | c3_hist->SetMarkerStyle(25); | |
1403 | c3_hist->SetMarkerColor(2); | |
1404 | c3_hist->SetLineColor(2); | |
1405 | c3_hist->SetMaximum(3); | |
5130881c | 1406 | c3_hist->SetMinimum(.7); |
9e040a65 | 1407 | if(!MCcase) c3_hist->Draw("same"); |
5130881c | 1408 | legend2->AddEntry(c3_hist,"#font[12]{#bf{c}}_{3} (cumulant correlation)","p"); |
9e040a65 | 1409 | c3_hist->Draw(); |
1410 | ||
0e610a08 | 1411 | //for(int i=1; i<=15; i++) cout<<c3_hist->GetBinContent(i)<<", "; |
1412 | //cout<<endl; | |
b2c0d4ba | 1413 | //for(int i=1; i<=15; i++) cout<<c3_hist->GetBinError(i)<<", "; |
1414 | //cout<<endl; | |
1415 | ||
9e040a65 | 1416 | const int npar_c3=5;// 10 |
1417 | TMinuit MyMinuit_c3(npar_c3); | |
1418 | MyMinuit_c3.SetFCN(fcn_c3); | |
b2c0d4ba | 1419 | int ierflg_c3=0; |
5130881c | 1420 | double arglist_c3 = 0; |
b2c0d4ba | 1421 | MyMinuit_c3.mnexcm("SET NOWarnings",&arglist_c3,1, ierflg_c3); |
5130881c | 1422 | arglist_c3 = -1; |
b2c0d4ba | 1423 | MyMinuit_c3.mnexcm("SET PRint",&arglist_c3,1, ierflg_c3); |
1424 | //arglist_c3=2;// improve Minimization Strategy (1 is default) | |
1425 | //MyMinuit_c3.mnexcm("SET STR",&arglist_c3,1,ierflg_c3); | |
1426 | //arglist_c3 = 0; | |
1427 | //MyMinuit_c3.mnexcm("SCAN", &arglist_c3,1,ierflg_c3); | |
1428 | arglist_c3 = 5000; | |
1429 | MyMinuit_c3.mnexcm("MIGRAD", &arglist_c3 ,1,ierflg_c3); | |
1430 | ||
0e610a08 | 1431 | TF1 *c3Fit1D_Base=new TF1("c3Fit1D_Base","[0]*(1+[1]*exp(-pow([2]*x/0.19733,[3]) * [4]/2.))",0,1); |
1432 | c3Fit1D_Base->FixParameter(3,ExpPower); | |
1433 | TF1 *c3Fit1D_Expan; | |
1434 | if(Gaussian) { | |
1435 | c3Fit1D_Base->FixParameter(4,1.0); | |
1436 | c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)),3))",0,1); | |
1437 | }else { | |
1438 | c3Fit1D_Base->FixParameter(4,sqrt(3.)); | |
1439 | c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,1) * sqrt(3.)/2.) * pow(1 + [3]*([2]*x/sqrt(3.)/0.19733 - 1) + [4]/2*(pow([2]*x/sqrt(3.)/0.19733,2) - 4*[2]*x/sqrt(3.)/0.19733 + 2),3))",0,1); | |
1440 | } | |
9e040a65 | 1441 | double OutputPar_c3[npar_c3]={0}; |
1442 | double OutputPar_c3_e[npar_c3]={0}; | |
1443 | ||
1444 | double par_c3[npar_c3]; // the start values | |
1445 | double stepSize_c3[npar_c3]; // step sizes | |
1446 | double minVal_c3[npar_c3]; // minimum bound on parameter | |
1447 | double maxVal_c3[npar_c3]; // maximum bound on parameter | |
1448 | string parName_c3[npar_c3]; | |
9e040a65 | 1449 | if(SameCharge && !MCcase){ |
b2c0d4ba | 1450 | for(int ft=0; ft<2; ft++){// Gaussian, Edgeworth |
0e610a08 | 1451 | if(ft==1 && !IncludeExpansion) continue; |
1452 | par_c3[0] = 1.0; par_c3[1] = .5; par_c3[2] = 5; par_c3[3] = kappa3; par_c3[4] = kappa4;// was par_c3[3] = 0.; par_c3[4] = 0.; | |
b2c0d4ba | 1453 | stepSize_c3[0] = 0.01; stepSize_c3[1] = 0.1; stepSize_c3[2] = 0.2; stepSize_c3[3] = 0.01; stepSize_c3[4] = 0.01; |
0e610a08 | 1454 | minVal_c3[0] = 0.95; minVal_c3[1] = 0.2; minVal_c3[2] = 0.5; minVal_c3[3] = -1; minVal_c3[4] = -1;// was minVal_c3[3] = -1; minVal_c3[4] = -1; |
1455 | maxVal_c3[0] = 1.1; maxVal_c3[1] = 2.5; maxVal_c3[2] = 15.; maxVal_c3[3] = +1; maxVal_c3[4] = +1;// was minVal_c3[3] = +1; minVal_c3[4] = +1; | |
1456 | parName_c3[0] = "N"; parName_c3[1] = "#lambda_{3}"; parName_c3[2] = "R_{inv}"; parName_c3[3] = "coeff_{3}"; parName_c3[4] = "coeff_{4}"; | |
1457 | c3Fit1D_Base->SetParName(0,"N"); c3Fit1D_Base->SetParName(1,"#lambda_{3}"); c3Fit1D_Base->SetParName(2,"R_{inv}"); | |
b2c0d4ba | 1458 | for (int i=0; i<npar_c3; i++){ |
1459 | MyMinuit_c3.DefineParameter(i, parName_c3[i].c_str(), par_c3[i], stepSize_c3[i], minVal_c3[i], maxVal_c3[i]); | |
9e040a65 | 1460 | } |
0e610a08 | 1461 | // |
1462 | if(!FixExpansionAvg) { | |
1463 | maxVal_c3[1] = 2.0; | |
1464 | double RadiusFix; | |
1465 | if(CollisionType==0) RadiusFix = 10 - 7*Mbin/15.; | |
1466 | else RadiusFix = 2.5 - 1.25*(Mbin-12)/7.; | |
1467 | MyMinuit_c3.DefineParameter(2, parName_c3[2].c_str(), RadiusFix, stepSize_c3[2], minVal_c3[2], maxVal_c3[2]); MyMinuit_c3.FixParameter(2); | |
1468 | MyMinuit_c3.DefineParameter(1, parName_c3[1].c_str(), 2.0, stepSize_c3[1], minVal_c3[1], maxVal_c3[1]); MyMinuit_c3.FixParameter(1); | |
1469 | } | |
1470 | // kappa_3 and kappa_4 | |
b2c0d4ba | 1471 | if(ft==0){// Gaussian |
1472 | MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), 0., stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3); | |
1473 | MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), 0., stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4); | |
1474 | }else{// Edgeworth | |
0e610a08 | 1475 | if(FixExpansionAvg){ |
b2c0d4ba | 1476 | MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), kappa3, stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3); |
1477 | MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), kappa4, stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4); | |
1478 | }else{ | |
1479 | Int_t err=0; | |
1480 | Double_t tmp[1]; | |
1481 | tmp[0] = 3+1; | |
1482 | MyMinuit_c3.mnexcm( "RELease", tmp, 1, err );// kappa3 | |
1483 | tmp[0] = 4+1; | |
1484 | MyMinuit_c3.mnexcm( "RELease", tmp, 1, err );// kappa4 | |
1485 | } | |
1486 | } | |
0e610a08 | 1487 | |
9e040a65 | 1488 | |
b2c0d4ba | 1489 | ///////////////////////////////////////////////////////////// |
1490 | // Do the minimization! | |
1491 | cout<<"Start Three-d fit"<<endl; | |
1492 | MyMinuit_c3.Migrad();// Minuit's best minimization algorithm | |
1493 | cout<<"End Three-d fit"<<endl; | |
1494 | ///////////////////////////////////////////////////////////// | |
1495 | MyMinuit_c3.mnexcm("SHOw PARameters", &arglist_c3, 1, ierflg); | |
1496 | cout<<"c3 Fit: Chi2/NDF = "<<Chi2_c3/(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl; | |
1497 | ||
1498 | for(int i=0; i<npar_c3; i++){ | |
1499 | MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]); | |
1500 | } | |
1501 | if(ft==0){ | |
0e610a08 | 1502 | c3Fit1D_Base->SetParameter(0,OutputPar_c3[0]); c3Fit1D_Base->SetParError(0,OutputPar_c3_e[0]); |
1503 | c3Fit1D_Base->SetParameter(1,OutputPar_c3[1]); c3Fit1D_Base->SetParError(1,OutputPar_c3_e[1]); | |
1504 | c3Fit1D_Base->SetParameter(2,OutputPar_c3[2]); c3Fit1D_Base->SetParError(2,OutputPar_c3_e[2]); | |
1505 | c3Fit1D_Base->SetLineStyle(2); | |
1506 | c3Fit1D_Base->Draw("same"); | |
b2c0d4ba | 1507 | }else{ |
0e610a08 | 1508 | c3Fit1D_Expan->SetParameter(0,OutputPar_c3[0]); c3Fit1D_Expan->SetParError(0,OutputPar_c3_e[0]); |
1509 | c3Fit1D_Expan->SetParameter(1,OutputPar_c3[1]); c3Fit1D_Expan->SetParError(1,OutputPar_c3_e[1]); | |
1510 | c3Fit1D_Expan->SetParameter(2,OutputPar_c3[2]); c3Fit1D_Expan->SetParError(2,OutputPar_c3_e[2]); | |
1511 | c3Fit1D_Expan->SetParameter(3,OutputPar_c3[3]); c3Fit1D_Expan->SetParError(3,OutputPar_c3_e[3]); | |
1512 | c3Fit1D_Expan->SetParameter(4,OutputPar_c3[4]); c3Fit1D_Expan->SetParError(4,OutputPar_c3_e[4]); | |
1513 | c3Fit1D_Expan->SetLineStyle(1); | |
1514 | c3Fit1D_Expan->Draw("same"); | |
1515 | double lambdaStar=0, lambdaStar_e=0; | |
1516 | if(Gaussian){ | |
1517 | lambdaStar = OutputPar_c3[1]*pow(1 + OutputPar_c3[4]/8.,3); | |
1518 | lambdaStar_e = pow(OutputPar_c3_e[1]*pow(1 + OutputPar_c3[4]/8.,3),2); | |
1519 | lambdaStar_e = sqrt(lambdaStar_e); | |
1520 | }else{ | |
1521 | lambdaStar = OutputPar_c3[1]*pow(1 - OutputPar_c3[3] + OutputPar_c3[4],3); | |
1522 | lambdaStar_e = pow(OutputPar_c3_e[1]*pow(1 - OutputPar_c3[3] + OutputPar_c3[4],3),2); | |
1523 | lambdaStar_e = sqrt(lambdaStar_e); | |
1524 | } | |
1525 | cout<<"lambda*_3 = "<<lambdaStar<<" +- "<<lambdaStar_e<<endl; | |
b2c0d4ba | 1526 | } |
0e610a08 | 1527 | //if(ft==0) {MainFitParams[1]=OutputPar_c3[2]; MainFitParams[3]=OutputPar_c3[1];} |
1528 | //else {MainFitParams[5]=OutputPar_c3[2]; MainFitParams[7]=OutputPar_c3[1];} | |
b2c0d4ba | 1529 | }// fit type |
1530 | }// SC and !MC | |
0e610a08 | 1531 | |
1532 | // project 3D EW fit onto 1D histogram | |
1533 | if(SameCharge && !MCcase){ | |
1534 | for(int i=2; i<=BINLIMIT_3; i++){// bin number | |
1535 | double Q12 = BinCenters[i-1];// true center | |
1536 | for(int j=2; j<=BINLIMIT_3; j++){// bin number | |
1537 | double Q13 = BinCenters[j-1];// true center | |
1538 | for(int k=2; k<=BINLIMIT_3; k++){// bin number | |
1539 | double Q23 = BinCenters[k-1];// true center | |
1540 | // | |
1541 | double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2)); | |
1542 | ||
1543 | if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible | |
1544 | if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible | |
1545 | ||
1546 | double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics) | |
1547 | if(TERM5==0) continue; | |
1548 | // | |
1549 | if(AddCC){ | |
1550 | if(CHARGE==-1){// base is (SC,MC,MC) | |
1551 | TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(j,k,i); | |
1552 | }else {// base is (MC,MC,SC) | |
1553 | TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(k,i,j); | |
1554 | } | |
1555 | } | |
1556 | // | |
1557 | int MRC_Q3bin = int(Q3/0.01) + 1; | |
1558 | TERM5 *= (MomRes1d[0][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1; | |
1559 | // | |
1560 | double radius_exp = c3Fit1D_Expan->GetParameter(2)/FmToGeV; | |
1561 | double arg12 = Q12*radius_exp; | |
1562 | double arg13 = Q13*radius_exp; | |
1563 | double arg23 = Q23*radius_exp; | |
1564 | ||
1565 | Float_t EW12 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12); | |
1566 | EW12 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12); | |
1567 | Float_t EW13 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13); | |
1568 | EW13 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12); | |
1569 | Float_t EW23 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23); | |
1570 | EW23 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12); | |
1571 | // | |
1572 | double cumulant_fitvalue=0; | |
1573 | cumulant_fitvalue = c3Fit1D_Expan->GetParameter(0)*(c3Fit1D_Expan->GetParameter(1)*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5); | |
1574 | ||
1575 | GenSignalExpected_num->Fill(Q3, cumulant_fitvalue); | |
1576 | GenSignalExpected_den->Fill(Q3, TERM5); | |
1577 | /////////////////////////////////////////////////////////// | |
1578 | ||
1579 | ||
1580 | ||
1581 | } | |
1582 | } | |
1583 | } | |
1584 | } | |
1585 | GenSignalExpected_num->Sumw2(); | |
1586 | GenSignalExpected_num->Divide(GenSignalExpected_den); | |
5be394f6 | 1587 | |
0e610a08 | 1588 | TSpline3 *c3Fit1D_ExpanSpline = new TSpline3(GenSignalExpected_num); |
1589 | if(CollisionType==0) c3Fit1D_ExpanSpline->SetLineColor(1); | |
1590 | if(CollisionType==1) c3Fit1D_ExpanSpline->SetLineColor(2); | |
1591 | if(CollisionType==2) c3Fit1D_ExpanSpline->SetLineColor(4); | |
1592 | c3Fit1D_ExpanSpline->SetLineWidth(2); | |
1593 | double xpoints[1000], ypoints[1000]; | |
1594 | bool splineOnly=kFALSE; | |
1595 | for(int iii=0; iii<1000; iii++){ | |
1596 | xpoints[iii] = 0 + (iii+0.5)*0.001; | |
1597 | if(!splineOnly){ | |
1598 | ypoints[iii] = c3Fit1D_Expan->Eval(xpoints[iii]); | |
1599 | if(c3Fit1D_Expan->Eval(xpoints[iii])<2 && fabs(c3Fit1D_Expan->Eval(xpoints[iii])-c3Fit1D_ExpanSpline->Eval(xpoints[iii])) < 0.01) splineOnly=kTRUE; | |
1600 | } | |
1601 | else {ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]); splineOnly=kTRUE;} | |
1602 | } | |
1603 | TGraph *gr_c3Spline = new TGraph(1000,xpoints,ypoints); | |
1604 | gr_c3Spline->SetLineWidth(2); | |
1605 | if(CollisionType==0) gr_c3Spline->SetLineColor(1); | |
1606 | if(CollisionType==1) gr_c3Spline->SetLineColor(2); | |
1607 | if(CollisionType==2) gr_c3Spline->SetLineColor(4); | |
1608 | gr_c3Spline->Draw("c same"); | |
1609 | //c3Fit1D_ExpanSpline->Draw("same"); | |
1610 | ||
1611 | ||
5be394f6 | 1612 | // uncomment to display fit box |
0e610a08 | 1613 | //c3_hist->GetListOfFunctions()->Add(c3Fit1D_Base); |
5be394f6 | 1614 | |
1615 | // Fit Range Systematics | |
1616 | //for(int i=0; i<8; i++){cout<<MainFitParams[i]<<", ";} | |
1617 | //cout<<endl; | |
1618 | //for(int i=0; i<4; i++){cout<<int(100*(MainFitParams[i]-RefMainFitParams[i])/MainFitParams[i])<<"$\\%$ & ";}// Gaussian | |
1619 | //cout<<endl; | |
1620 | //for(int i=4; i<8; i++){cout<<int(100*(MainFitParams[i]-RefMainFitParams[i])/MainFitParams[i])<<"$\\%$ & ";}// Edgeworth | |
1621 | //cout<<endl; | |
1622 | ||
b2c0d4ba | 1623 | |
5be394f6 | 1624 | /*TString *lamName=new TString("#lambda_{3}="); |
1625 | TString *RName=new TString("#R_{inv,3}="); | |
0e610a08 | 1626 | *lamName += int(100*c3Fit1D_Base->GetParameter(1))/100.; lamName->Append(" #pm "); *lamName += int(100*c3Fit1D_Base->GetParError(1))/100.; |
1627 | *RName += round(100*c3Fit1D_Base->GetParameter(2))/100; RName->Append(" #pm "); *RName += int(100*c3Fit1D_Base->GetParError(2))/100.; RName->Append(" fm"); | |
5be394f6 | 1628 | TLatex *fitBox1=new TLatex(0.5,0.9,lamName->Data()); |
b2c0d4ba | 1629 | fitBox1->SetNDC(kTRUE); |
1630 | fitBox1->Draw(); | |
5be394f6 | 1631 | TLatex *fitBox2=new TLatex(0.5,0.8,RName->Data()); |
b2c0d4ba | 1632 | fitBox2->SetNDC(kTRUE); |
5be394f6 | 1633 | fitBox2->Draw();*/ |
1634 | //TLatex *fitBox3=new TLatex(0.5,0.5,"<N_{rec,pions}>=103"); | |
1635 | //fitBox3->SetNDC(kTRUE); | |
1636 | //fitBox3->Draw(); | |
1637 | ||
0e610a08 | 1638 | //TPaveStats *c3Stats=(TPaveStats*)c3Fit1D_Base->FindObject("stats"); |
5130881c | 1639 | //c3Stats->SetX1NDC(0.15); |
1640 | //c3Stats->SetX2NDC(0.52); | |
1641 | //c3Stats->SetY1NDC(.2); | |
1642 | //c3Stats->SetY2NDC(.3); | |
1643 | //c3Stats->Draw("same"); | |
9e040a65 | 1644 | |
9e040a65 | 1645 | |
5be394f6 | 1646 | // The below 1D fit method has less well defined bin centers |
b2c0d4ba | 1647 | //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,1); |
9e040a65 | 1648 | //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + (-0.12/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + (0.43/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1); |
1649 | //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1); | |
1650 | //c3Fit1D->SetLineColor(4); | |
b2c0d4ba | 1651 | //c3Fit1D->SetParameter(0,1); c3Fit1D->SetParName(0,"N"); |
1652 | //c3Fit1D->SetParameter(1,0.5); c3Fit1D->SetParName(1,"#lambda_{3}"); | |
1653 | //c3Fit1D->SetParameter(2,5); c3Fit1D->SetParName(2,"R_{inv}"); | |
1654 | //c3Fit1D->SetParLimits(2,0.5,15); | |
1655 | //c3Fit1D->SetLineColor(1); | |
1656 | // | |
1657 | //c3Fit1D->FixParameter(0,OutputPar_c3[0]); | |
1658 | //c3Fit1D->FixParameter(1,OutputPar_c3[1]); | |
1659 | //c3Fit1D->FixParameter(2,OutputPar_c3[2]); | |
0e610a08 | 1660 | //c3Fit1D->SetParName(3,"coeff_{3}"); c3Fit1D->SetParName(4,"coeff_{4}"); |
9e040a65 | 1661 | //c3Fit1D->SetParameter(3,.24); c3Fit1D->SetParameter(3,.16); |
b2c0d4ba | 1662 | //c3_hist->Fit(c3Fit1D,"IMEN","",0.0,Q3Limit); |
1663 | //c3Fit1D->Draw("same"); | |
9e040a65 | 1664 | //cout<<"1d fit: Chi2/NDF = "<<c3Fit1D->GetChisquare()/c3Fit1D->GetNDF()<<endl; |
1665 | //dentimesFit_c3->DrawCopy("l same"); | |
1666 | // | |
1667 | //c3Fit1D->FixParameter(0, fitcopy_c3->GetParameter(0)); | |
1668 | //c3Fit1D->FixParameter(1, fitcopy_c3->GetParameter(1)); | |
1669 | //c3Fit1D->FixParameter(2, fitcopy_c3->GetParameter(2)); | |
5130881c | 1670 | |
b2c0d4ba | 1671 | //TF1 *c3Fit1D_2=new TF1("c3Fit1D_2","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,1); |
1672 | //c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)),3))",0,1); | |
1673 | //c3Fit1D->FixParameter(0,OutputPar_c3[0]); | |
1674 | //c3Fit1D->FixParameter(1,OutputPar_c3[1]); | |
1675 | //c3Fit1D->FixParameter(2,OutputPar_c3[2]); | |
1676 | //c3Fit1D->FixParameter(3,OutputPar_c3[3]); | |
1677 | //c3Fit1D->FixParameter(4,OutputPar_c3[4]); | |
1678 | //c3Fit1D->SetLineColor(4); | |
1679 | //c3Fit1D->Draw("same"); | |
1680 | ||
5130881c | 1681 | |
1682 | //MomRes1d[1][0]->Draw(); | |
1683 | ||
9e040a65 | 1684 | //for(int i=1; i<50; i++) cout<<c3_hist->GetBinContent(i)<<", "; |
1685 | //cout<<endl; | |
1686 | //for(int i=1; i<50; i++) cout<<c3_hist->GetBinError(i)<<", "; | |
1687 | // pp M2, pi- | |
0e610a08 | 1688 | |
9e040a65 | 1689 | |
1690 | Cterm1->Draw("same"); | |
1691 | legend2->Draw("same"); | |
1692 | ||
5be394f6 | 1693 | /*if(SameCharge==kTRUE && MCcase==kFALSE){ |
5130881c | 1694 | TString *savename=new TString("C3c3_SC_"); |
1695 | if(CollisionType==0) savename->Append("PbPb_K"); | |
1696 | if(CollisionType==1) savename->Append("pPb_K"); | |
1697 | if(CollisionType==2) savename->Append("pp_K"); | |
1698 | *savename += EDbin; | |
1699 | savename->Append("_M"); | |
1700 | *savename += Mbin; | |
1701 | savename->Append(".eps"); | |
1702 | can2->SaveAs(savename->Data()); | |
1703 | }*/ | |
9e040a65 | 1704 | |
1705 | //for(int ii=0; ii<10; ii++) cout<<c3_hist->GetBinContent(ii+1)<<", "; | |
1706 | //Coul_GRiverside->Draw(); | |
1707 | //Coul_Riverside->Draw("same"); | |
0e610a08 | 1708 | /*TLegend *legend4 = new TLegend(.3,.8, .5,.95,NULL,"brNDC"); |
5130881c | 1709 | legend4->SetBorderSize(0); |
1710 | legend4->SetTextFont(42);//42 | |
1711 | legend4->SetTextSize(.04);// small .03; large .06 | |
9e040a65 | 1712 | legend4->SetFillColor(0); |
5130881c | 1713 | |
1714 | gPad->SetRightMargin(0.025); gPad->SetLeftMargin(0.1); gPad->SetTopMargin(0.02); | |
0e610a08 | 1715 | TwoFrac=lambda_PM; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5); |
1716 | gPad->SetGridx(0); | |
1717 | gPad->SetGridy(0); | |
1718 | ||
9e040a65 | 1719 | TH1D *c3_Extended = (TH1D*)Three_1d[CB1][CB2][CB3][SCBin][0]->Clone(); |
1720 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], -( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )); | |
1721 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][1], -(1-OneFrac)); | |
1722 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][2], -(1-OneFrac)); | |
1723 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][3], -(1-OneFrac)); | |
1724 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], (1-OneFrac)*3*(1-TwoFrac)); | |
1725 | c3_Extended->Scale(1/ThreeFrac); | |
1726 | // | |
1727 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][1], -1/TwoFrac); | |
1728 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][2], -1/TwoFrac); | |
1729 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][3], -1/TwoFrac); | |
1730 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], 3*(1-TwoFrac)/TwoFrac); | |
1731 | c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], 3); | |
1732 | // | |
1733 | c3_Extended->Divide(Three_1d[CB1][CB2][CB3][SCBin][4]); | |
5130881c | 1734 | c3_Extended->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})"); |
0e610a08 | 1735 | c3_Extended->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}^{#pm#pm#pm}"); |
1736 | c3_Extended->SetMinimum(0.9); c3_Extended->SetMaximum(2.0); | |
9e040a65 | 1737 | c3_Extended->SetMarkerStyle(24); |
1738 | c3_Extended->SetMarkerColor(2); | |
0e610a08 | 1739 | c3_Extended->SetLineColor(2); |
1740 | ||
9e040a65 | 1741 | // |
1742 | TH1D *C3_Extended = (TH1D*)Three_1d[CB1][CB2][CB3][SCBin][0]->Clone(); | |
1743 | C3_Extended->Divide(Three_1d[CB1][CB2][CB3][SCBin][4]); | |
5130881c | 1744 | C3_Extended->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/c)"); |
1745 | C3_Extended->GetYaxis()->SetTitle("#font[12]{C}_{3}"); | |
0e610a08 | 1746 | c3_Extended->GetXaxis()->SetTitleOffset(0.83); |
1747 | c3_Extended->GetYaxis()->SetTitleOffset(0.8); | |
1748 | c3_Extended->GetYaxis()->SetTitleSize(0.05); | |
1749 | c3_Extended->GetXaxis()->SetTitleSize(0.05); | |
9e040a65 | 1750 | C3_Extended->SetMarkerStyle(20); |
1751 | C3_Extended->SetMarkerColor(4); | |
0e610a08 | 1752 | c3_Extended->GetXaxis()->SetRangeUser(0,2); |
1753 | c3_Extended->SetMarkerStyle(20); | |
1754 | c3_Extended->Draw();*/ | |
9e040a65 | 1755 | //legend4->AddEntry(c3_Extended,"c_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p"); |
1756 | //legend4->AddEntry(C3_Extended,"C_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p"); | |
1757 | //legend4->Draw("same"); | |
0e610a08 | 1758 | //Unity->Draw("same"); |
5130881c | 1759 | |
0e610a08 | 1760 | //Specif_System->Draw("same"); |
1761 | //Specif_KT3->Draw("same"); | |
1762 | //Specif_FSI->Draw("same"); | |
1763 | //Specif_Disclaimer->Draw("same"); | |
1764 | ||
1765 | ||
9e040a65 | 1766 | /* |
5130881c | 1767 | gPad->SetGridx(0); gPad->SetGridy(0); |
1768 | int binLow=20; | |
9e040a65 | 1769 | int binHigh=50; |
5130881c | 1770 | //TwoFrac=0.9; OneFrac=pow(TwoFrac,0.5); ThreeFrac=pow(TwoFrac,1.5); |
9e040a65 | 1771 | TH1D *K0sProjection = (TH1D*)(Three_3d[0][0][1][SCBin][0]->ProjectionY("K0sProjection",binLow,binHigh,binLow,binHigh)); |
1772 | TH1D *K0sProjection_term1 = (TH1D*)(Three_3d[0][0][1][SCBin][0]->ProjectionY("K0sProjection_term1",binLow,binHigh,binLow,binHigh)); | |
1773 | TH1D *K0sProjection_term2 = (TH1D*)(Three_3d[0][0][1][SCBin][1]->ProjectionY("K0sProjection_term2",binLow,binHigh,binLow,binHigh)); | |
1774 | TH1D *K0sProjection_term3 = (TH1D*)(Three_3d[0][0][1][SCBin][2]->ProjectionY("K0sProjection_term3",binLow,binHigh,binLow,binHigh)); | |
1775 | TH1D *K0sProjection_term4 = (TH1D*)(Three_3d[0][0][1][SCBin][3]->ProjectionY("K0sProjection_term4",binLow,binHigh,binLow,binHigh)); | |
1776 | TH1D *K0sProjection_term5 = (TH1D*)(Three_3d[0][0][1][SCBin][4]->ProjectionY("K0sProjection_term5",binLow,binHigh,binLow,binHigh)); | |
1777 | K0sProjection_term1->Add(K0sProjection_term5, -( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )); | |
1778 | K0sProjection_term1->Add(K0sProjection_term2, -(1-OneFrac)); | |
1779 | K0sProjection_term1->Add(K0sProjection_term3, -(1-OneFrac)); | |
1780 | K0sProjection_term1->Add(K0sProjection_term4, -(1-OneFrac)); | |
1781 | K0sProjection_term1->Add(K0sProjection_term5, (1-OneFrac)*3*(1-TwoFrac)); | |
1782 | K0sProjection_term1->Scale(1/ThreeFrac); | |
1783 | // | |
1784 | K0sProjection_term1->Add(K0sProjection_term2, -1/TwoFrac); | |
1785 | K0sProjection_term1->Add(K0sProjection_term3, -1/TwoFrac); | |
1786 | K0sProjection_term1->Add(K0sProjection_term4, -1/TwoFrac); | |
1787 | K0sProjection_term1->Add(K0sProjection_term5, 3*(1-TwoFrac)/TwoFrac); | |
1788 | K0sProjection_term1->Add(K0sProjection_term5, 3); | |
1789 | // | |
1790 | K0sProjection->Divide(K0sProjection_term5); | |
1791 | K0sProjection_term1->Divide(K0sProjection_term5); | |
1792 | K0sProjection_term1->GetXaxis()->SetRangeUser(0,.5); | |
1793 | K0sProjection_term1->SetMarkerStyle(24); | |
1794 | K0sProjection_term1->SetMarkerColor(4); | |
1795 | K0sProjection->SetMarkerStyle(20); | |
1796 | K0sProjection->SetMarkerColor(4); | |
5130881c | 1797 | K0sProjection->SetMinimum(0.92); |
1798 | K0sProjection->GetXaxis()->SetTitle("#font[12]{q_{13}^{#pm#mp}} (GeV/#font[12]{c})"); | |
1799 | K0sProjection->GetYaxis()->SetTitle("#font[12]{C_{3}} or #font[12]{#bf{c}_{3}}"); | |
1800 | K0sProjection->GetXaxis()->SetTitleOffset(1.2); K0sProjection->GetYaxis()->SetTitleOffset(1.3); | |
1801 | ||
1802 | K0sProjection->Draw(); | |
5be394f6 | 1803 | // Sys errors |
1804 | TH1D *Sys_K0sProjection=(TH1D*)K0sProjection->Clone(); | |
1805 | TH1D *Sys_K0sProjection_term1=(TH1D*)K0sProjection_term1->Clone(); | |
1806 | Sys_K0sProjection->SetMarkerSize(0); Sys_K0sProjection_term1->SetMarkerSize(0); | |
1807 | Sys_K0sProjection->SetFillColor(kBlue-10); Sys_K0sProjection_term1->SetFillColor(kBlue-10); | |
1808 | Sys_K0sProjection->SetFillStyle(1000); Sys_K0sProjection_term1->SetFillStyle(1000); | |
1809 | for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) { | |
1810 | Sys_K0sProjection->SetBinError(i+1, 0.01 * Sys_K0sProjection->GetBinContent(i+1)); | |
1811 | Sys_K0sProjection_term1->SetBinError(i+1, 0.01 * Sys_K0sProjection_term1->GetBinContent(i+1)); | |
1812 | } | |
1813 | for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection->GetBinContent(i+1)<<", "; | |
1814 | cout<<endl; | |
1815 | for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection->GetBinError(i+1)<<", "; | |
1816 | cout<<endl; | |
1817 | for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection_term1->GetBinContent(i+1)<<", "; | |
1818 | cout<<endl; | |
1819 | for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection_term1->GetBinError(i+1)<<", "; | |
1820 | cout<<endl; | |
1821 | ||
1822 | Sys_K0sProjection->Draw("E2 same"); | |
1823 | Sys_K0sProjection_term1->Draw("E2 same"); | |
1824 | K0sProjection->Draw("same"); | |
9e040a65 | 1825 | K0sProjection_term1->Draw("same"); |
5be394f6 | 1826 | |
5130881c | 1827 | legend4->AddEntry(K0sProjection,"#font[12]{C_{3}^{#pm#pm#mp}} projection","p"); |
1828 | legend4->AddEntry(K0sProjection_term1,"#font[12]{#bf{c}_{3}^{#pm#pm#mp}} projection","p"); | |
9e040a65 | 1829 | legend4->Draw("same"); |
5130881c | 1830 | |
9e040a65 | 1831 | |
9e040a65 | 1832 | */ |
1833 | ||
0e610a08 | 1834 | |
9e040a65 | 1835 | ////////////////////////////////////////////////////////////////////////////////////// |
1836 | ||
1837 | ||
1838 | ||
1839 | ||
1840 | ||
1841 | ||
1842 | if(SaveToFile){ | |
1843 | TString *savefilename = new TString("OutFiles/OutFile"); | |
1844 | if(CollisionType==0) savefilename->Append("PbPb"); | |
1845 | else if(CollisionType==1) savefilename->Append("pPb"); | |
1846 | else savefilename->Append("pp"); | |
1847 | // | |
1848 | if(MCcase) savefilename->Append("MonteCarlo"); | |
1849 | // | |
1850 | if(SameCharge) savefilename->Append("SC"); | |
1851 | else savefilename->Append("MC"); | |
1852 | // | |
1853 | if(CHARGE==1) savefilename->Append("Pos"); | |
1854 | else savefilename->Append("Neg"); | |
1855 | // | |
0e610a08 | 1856 | |
9e040a65 | 1857 | savefilename->Append("Kt3_"); |
1858 | *savefilename += EDbin+1; | |
1859 | ||
1860 | savefilename->Append("_M"); | |
1861 | *savefilename += Mbin; | |
1862 | savefilename->Append(".root"); | |
1863 | TFile *savefile = new TFile(savefilename->Data(),"RECREATE"); | |
1864 | can1->Write("can1"); | |
1865 | C2_ss->Write("C2_ss"); | |
1866 | C2_os->Write("C2_os"); | |
1867 | fitC2ss_h->Write("fitC2ss_h"); | |
0e610a08 | 1868 | fitC2ss_Base->Write("fitC2ss_Base"); |
1869 | fitC2ss_Expan->Write("fitC2ss_Expan"); | |
9e040a65 | 1870 | Cterm1->Write("C3"); |
1871 | c3_hist->Write("c3"); | |
9e040a65 | 1872 | Combinatorics_1d->Write("Combinatorics_1d"); |
0e610a08 | 1873 | c3Fit1D_Base->Write("c3Fit1D_Base"); |
1874 | c3Fit1D_Expan->Write("c3Fit1D_Expan"); | |
1875 | c3Fit1D_ExpanSpline->Write("c3Fit1D_ExpanSpline"); | |
1876 | gr_c3Spline->Write("gr_c3Spline"); | |
9e040a65 | 1877 | MyMinuit_c3.Write("MyMinuit_c3"); |
9e040a65 | 1878 | // |
1879 | savefile->Close(); | |
1880 | } | |
1881 | ||
1882 | ||
1883 | } | |
1884 | void ReadCoulCorrections(int index){ | |
1885 | /////////////////////// | |
1886 | TString *fname2 = new TString("KFile.root"); | |
1887 | ||
1888 | TFile *File=new TFile(fname2->Data(),"READ"); | |
0e610a08 | 1889 | if(index>=12) cout<<"FSI index not acceptable in 2-particle Coulomb read"<<endl; |
9e040a65 | 1890 | TString *nameSS = new TString("K2ss_"); |
1891 | TString *nameOS = new TString("K2os_"); | |
1892 | *nameSS += index; | |
1893 | *nameOS += index; | |
1894 | TString *nameSS_2 = new TString("K2ss_"); | |
1895 | TString *nameOS_2 = new TString("K2os_"); | |
0e610a08 | 1896 | *nameSS_2 += index+1; *nameOS_2 += index+1; |
1897 | /*if(index<9) {*nameSS_2 += index+1; *nameOS_2 += index+1;} | |
1898 | else {*nameSS_2 += index; *nameOS_2 += index;}*/ | |
9e040a65 | 1899 | TH1D *temp_ss = (TH1D*)File->Get(nameSS->Data()); |
1900 | TH1D *temp_os = (TH1D*)File->Get(nameOS->Data()); | |
1901 | TH1D *temp_ss_2 = (TH1D*)File->Get(nameSS_2->Data()); | |
1902 | TH1D *temp_os_2 = (TH1D*)File->Get(nameOS_2->Data()); | |
1903 | ||
1904 | CoulCorr2SS = (TH1D*)temp_ss->Clone(); | |
1905 | CoulCorr2OS = (TH1D*)temp_os->Clone(); | |
1906 | CoulCorr2SS_2 = (TH1D*)temp_ss_2->Clone(); | |
1907 | CoulCorr2OS_2 = (TH1D*)temp_os_2->Clone(); | |
0e610a08 | 1908 | |
9e040a65 | 1909 | // |
1910 | CoulCorr2SS->SetDirectory(0); | |
1911 | CoulCorr2OS->SetDirectory(0); | |
1912 | CoulCorr2SS_2->SetDirectory(0); | |
1913 | CoulCorr2OS_2->SetDirectory(0); | |
1914 | ||
1915 | File->Close(); | |
1916 | } | |
1917 | double CoulCorr2(int chargeproduct, double Q2){// returns K2 | |
1918 | int indexL=0; | |
1919 | int indexH=0; | |
1920 | double slope=0; | |
1921 | double value1=1.0, value2=1.0; | |
1922 | ||
1923 | indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1))); | |
1924 | indexH = indexL+1; | |
0e610a08 | 1925 | |
1926 | ||
1927 | if(indexH >= CoulCorr2SS->GetNbinsX()) { | |
1928 | return (Gamov(chargeproduct, Q2) * CoulCorr2SS->GetBinContent(CoulCorr2SS->GetNbinsX()-1) / Gamov(chargeproduct, CoulCorr2SS->GetXaxis()->GetBinCenter(CoulCorr2SS->GetNbinsX()-1))); | |
1929 | } | |
1930 | // Low index value | |
1931 | if(chargeproduct==1){ | |
1932 | slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1)); | |
1933 | slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1)); | |
1934 | value1 = slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1); | |
1935 | }else { | |
1936 | slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1)); | |
1937 | slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1)); | |
1938 | value1 = slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1); | |
1939 | } | |
1940 | // High index value | |
1941 | if(chargeproduct==1){ | |
1942 | slope = (CoulCorr2SS_2->GetBinContent(indexL+1) - CoulCorr2SS_2->GetBinContent(indexH+1)); | |
1943 | slope /= (CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexH+1)); | |
1944 | value2 = slope*(Q2 - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS_2->GetBinContent(indexL+1); | |
1945 | }else { | |
1946 | slope = (CoulCorr2OS_2->GetBinContent(indexL+1) - CoulCorr2OS_2->GetBinContent(indexH+1)); | |
1947 | slope /= (CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexH+1)); | |
1948 | value2 = slope*(Q2 - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS_2->GetBinContent(indexL+1); | |
1949 | } | |
1950 | ||
1951 | if(value1>0 && value2>0){// interpolation only done between 9 and 12 Mbin. | |
1952 | return (value1 + (value2-value1)/(2.));// always use the half-way point. | |
1953 | }else if(value1>0){ | |
1954 | return value1; | |
1955 | }else if(value2>0){ | |
1956 | return value2; | |
1957 | }else return 1.0; | |
1958 | ||
1959 | ||
1960 | // | |
1961 | // old way with Gaussian sources for low mult systems | |
1962 | /*if(indexH >= CoulCorr2SS->GetNbinsX()) { | |
1963 | return (Gamov(chargeproduct, Q2) * CoulCorr2SS->GetBinContent(CoulCorr2SS->GetNbinsX()-1) / Gamov(chargeproduct, CoulCorr2SS->GetXaxis()->GetBinCenter(CoulCorr2SS->GetNbinsX()-1))); | |
1964 | } | |
9e040a65 | 1965 | if(chargeproduct==1){ |
1966 | slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1)); | |
1967 | slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1)); | |
1968 | value1 = slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1); | |
1969 | }else { | |
1970 | slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1)); | |
1971 | slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1)); | |
1972 | value1 = slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1); | |
1973 | } | |
1974 | if(Mbin_def<=9 || Mbin_def>12) return value1; | |
1975 | else{// interpolate in K factor transition region | |
1976 | if(chargeproduct==1){ | |
1977 | slope = (CoulCorr2SS_2->GetBinContent(indexL+1) - CoulCorr2SS_2->GetBinContent(indexH+1)); | |
1978 | slope /= (CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexH+1)); | |
1979 | value2 = slope*(Q2 - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS_2->GetBinContent(indexL+1); | |
1980 | }else { | |
1981 | slope = (CoulCorr2OS_2->GetBinContent(indexL+1) - CoulCorr2OS_2->GetBinContent(indexH+1)); | |
1982 | slope /= (CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexH+1)); | |
1983 | value2 = slope*(Q2 - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS_2->GetBinContent(indexL+1); | |
1984 | } | |
0e610a08 | 1985 | |
5130881c | 1986 | if(value1>0 && value2>0){// interpolation only done between 9 and 12 Mbin. |
9e040a65 | 1987 | return (value1 + (Mbin_def-9)*(value2-value1)/(3.)); |
1988 | }else if(value1>0){ | |
1989 | return value1; | |
1990 | }else if(value2>0){ | |
1991 | return value2; | |
1992 | }else return 1.0; | |
1993 | ||
1994 | } | |
0e610a08 | 1995 | */ |
9e040a65 | 1996 | // |
1997 | ||
1998 | } | |
1999 | ||
2000 | //---------------------------------------------------------------------- | |
0e610a08 | 2001 | |
9e040a65 | 2002 | |
2003 | //________________________________________________________________________ | |
2004 | void fcnC2_Global(int& npar, double* deriv, double& f, double par[], int flag){ | |
2005 | ||
2006 | double qinvSS=0; | |
2007 | ||
2008 | double Rch=par[3]/FmToGeV; | |
9e040a65 | 2009 | double CSS=0, COS=0; |
2010 | double SumChi2=0; | |
0e610a08 | 2011 | double Expan=0; |
9e040a65 | 2012 | NFitPoints_C2global=0; |
2013 | if(LinkN) par[9]=par[0];// Link N factors | |
2014 | ||
2015 | for(int i=1; i<BINRANGE_C2global; i++){ | |
2016 | ||
2017 | qinvSS = BinCenters[i]; | |
2018 | if(qinvSS > Q2Limit) continue; | |
2019 | ||
9e040a65 | 2020 | // |
2021 | double arg=qinvSS*Rch; | |
0e610a08 | 2022 | if(Gaussian){// Edgeworth expansion |
2023 | Expan = 1 + par[5]/(6.*pow(2,1.5)) * (8.*pow(arg,3) - 12.*pow(arg,1)); | |
2024 | Expan += par[6]/(24.*pow(2,2)) * (16.*pow(arg,4) -48.*pow(arg,2) + 12); | |
2025 | Expan += par[7]/(120.*pow(2,2.5)) * (32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg); | |
2026 | Expan += par[8]/(720.*pow(2,3)) * (64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120); | |
2027 | }else{ | |
2028 | Expan = 1 + par[5] * (arg - 1); | |
2029 | Expan += par[6]/2. * (pow(arg,2) - 4*arg + 2); | |
2030 | Expan += par[7] * 0.; | |
2031 | Expan += par[8] * 0.; | |
2032 | } | |
2033 | //Expan = 1.0; | |
9e040a65 | 2034 | // |
7c6e4906 | 2035 | |
0e610a08 | 2036 | // old way without undilution |
2037 | /*CSS = 1 + exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
9e040a65 | 2038 | CSS *= par[1]*K2SS[i]; |
2039 | CSS += 1-par[1]; | |
0e610a08 | 2040 | CSS *= par[0];*/ |
2041 | // undilution method | |
2042 | //CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
2043 | //CSS *= par[0]; | |
2044 | // undilution method with correct normalization location (was previously applied to C2^QS) | |
2045 | CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
2046 | CSS *= lambda_PM * K2SS[i]; | |
2047 | CSS += 1-lambda_PM; | |
9e040a65 | 2048 | CSS *= par[0]; |
0e610a08 | 2049 | |
9e040a65 | 2050 | // |
2051 | COS = 1; | |
9e040a65 | 2052 | COS *= par[1]*K2OS[i]; |
2053 | COS += 1-par[1]; | |
2054 | COS *= par[9];// different Norm | |
2055 | // | |
2056 | if(C2ssFitting[i] > 0){ | |
2057 | if(IncludeSS) { | |
0e610a08 | 2058 | //double errorSS = sqrt(pow( (CSS/par[0] - (1-par[1]))/K2SS[i] * K2SS_e[i] * par[0],2) + pow(C2ssFitting_e[i],2)); |
2059 | double errorSS = sqrt(pow(C2ssFitting_e[i],2));// new way with undilution already done in C2ssFitting[i] | |
9e040a65 | 2060 | SumChi2 += pow((CSS - C2ssFitting[i])/errorSS,2); |
2061 | NFitPoints_C2global++; | |
2062 | } | |
2063 | } | |
2064 | if(IncludeOS) { | |
2065 | double errorOS = sqrt(pow( (COS/par[9] - (1-par[1]))/K2OS[i] * K2OS_e[i] * par[9],2) + pow(C2osFitting_e[i],2)); | |
9e040a65 | 2066 | SumChi2 += pow((COS - C2osFitting[i])/errorOS,2); |
2067 | NFitPoints_C2global++; | |
2068 | } | |
2069 | } | |
2070 | ||
2071 | ||
2072 | f = SumChi2; | |
2073 | Chi2_C2global = f; | |
2074 | ||
2075 | } | |
2076 | //________________________________________________________________________ | |
2077 | double C2ssFitFunction(Double_t *x, Double_t *par) | |
2078 | { | |
2079 | double Rch=par[3]/FmToGeV; | |
9e040a65 | 2080 | int qbin = int(fabs(x[0]/BinWidthQ2)); |
0e610a08 | 2081 | |
9e040a65 | 2082 | double qinvSS = BinCenters[qbin]; |
9e040a65 | 2083 | |
9e040a65 | 2084 | double arg=qinvSS*Rch; |
0e610a08 | 2085 | double Expan = 1; |
2086 | if(Gaussian){// Edgeworth expansion | |
2087 | Expan += par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1)); | |
2088 | Expan += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12); | |
2089 | Expan += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg); | |
2090 | Expan += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120); | |
2091 | }else{// Laguerre | |
2092 | Expan = 1 + par[5] * (arg - 1); | |
2093 | Expan += par[6]/2. * (pow(arg,2) - 4*arg + 2); | |
2094 | Expan += par[7] * 0.; | |
2095 | Expan += par[8] * 0.; | |
2096 | } | |
2097 | ||
2098 | // old way without undilution | |
2099 | /*double CSS = 1 + exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
9e040a65 | 2100 | double K2=1.0; |
2101 | if(qbin < BINRANGE_C2global) K2=K2SS[qbin]; | |
9e040a65 | 2102 | CSS *= par[1]*K2; |
2103 | CSS += 1-par[1]; | |
0e610a08 | 2104 | CSS *= par[0];*/ |
2105 | // undilution method | |
2106 | //double CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
2107 | //CSS *= par[0]; | |
2108 | // undilution method with correct normalization location (was previously applied to C2^QS) | |
2109 | double CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2); | |
2110 | CSS *= lambda_PM * K2SS[qbin]; | |
2111 | CSS += 1-lambda_PM; | |
9e040a65 | 2112 | CSS *= par[0]; |
0e610a08 | 2113 | |
9e040a65 | 2114 | return CSS; |
2115 | } | |
2116 | //________________________________________________________________________ | |
2117 | double C2osFitFunction(Double_t *x, Double_t *par) | |
2118 | { | |
2119 | if(LinkN) par[9]=par[0];// Link N factors | |
9e040a65 | 2120 | int qbin = int(fabs(x[0]/BinWidthQ2)); |
0e610a08 | 2121 | |
2122 | //double qinvOS = BinCenters[qbin]; | |
9e040a65 | 2123 | |
9e040a65 | 2124 | double COS = 1; |
2125 | double K2=1.0; | |
2126 | if(qbin < BINRANGE_C2global) K2=K2OS[qbin]; | |
9e040a65 | 2127 | COS *= par[1]*K2; |
2128 | COS += 1-par[1]; | |
2129 | COS *= par[9]; | |
2130 | return COS; | |
2131 | } | |
2132 | //__________________________________________________________________________ | |
2133 | ||
2134 | void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){ | |
2135 | ||
2136 | double q12=0, q13=0, q23=0; | |
0e610a08 | 2137 | double Expan12=0, Expan13=0, Expan23=0; |
9e040a65 | 2138 | double C=0; |
2139 | double Rch=par[2]/FmToGeV; | |
2140 | double SumChi2=0; | |
2141 | //double lnL=0, term1=0, term2=0; | |
2142 | NFitPoints_c3=0; | |
2143 | ||
2144 | for(int i=0; i<=BINLIMIT_3; i++){// q12 | |
2145 | for(int j=0; j<=BINLIMIT_3; j++){// q13 | |
2146 | for(int k=0; k<=BINLIMIT_3; k++){// q23 | |
2147 | ||
2148 | if(B_3[i][j][k] == 0) continue; | |
2149 | if(A_3[i][j][k] == 0) continue; | |
2150 | if(A_3_e[i][j][k] == 0) continue; | |
2151 | ||
2152 | q12 = BinCenters[i]; | |
2153 | q13 = BinCenters[j]; | |
2154 | q23 = BinCenters[k]; | |
2155 | double q3 = sqrt(pow(q12,2)+pow(q13,2)+pow(q23,2)); | |
2156 | if(q3 > Q3Limit) continue; | |
2157 | // | |
2158 | double arg12 = q12*Rch; | |
2159 | double arg13 = q13*Rch; | |
2160 | double arg23 = q23*Rch; | |
0e610a08 | 2161 | if(Gaussian){// Edgeworth expansion |
2162 | Expan12 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg12,3) - 12*pow(arg12,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg12,4) -48*pow(arg12,2) + 12); | |
2163 | Expan13 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg13,3) - 12*pow(arg13,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg13,4) -48*pow(arg13,2) + 12); | |
2164 | Expan23 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg23,3) - 12*pow(arg23,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg23,4) -48*pow(arg23,2) + 12); | |
2165 | }else{// Laguerre expansion | |
2166 | Expan12 = 1 + par[3]*(arg12 - 1) + par[4]/2.*(pow(arg12,2) - 4*arg12 + 2); | |
2167 | Expan13 = 1 + par[3]*(arg13 - 1) + par[4]/2.*(pow(arg13,2) - 4*arg13 + 2); | |
2168 | Expan23 = 1 + par[3]*(arg23 - 1) + par[4]/2.*(pow(arg23,2) - 4*arg23 + 2); | |
2169 | } | |
9e040a65 | 2170 | // |
0e610a08 | 2171 | C = 1 + par[1]*exp(-(pow(arg12,ExpPower)+pow(arg13,ExpPower)+pow(arg23,ExpPower))/2.)*Expan12*Expan13*Expan23; |
9e040a65 | 2172 | C *= par[0];// Norm |
2173 | ||
2174 | double error = pow(A_3_e[i][j][k]/B_3[i][j][k],2); | |
b2c0d4ba | 2175 | error += pow(sqrt(B_3[i][j][k])*A_3[i][j][k]/pow(B_3[i][j][k],2),2); |
9e040a65 | 2176 | error = sqrt(error); |
2177 | SumChi2 += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2); | |
2178 | // | |
2179 | /*if(A_3[i][j][k] >= 1) term1 = C*(A_3[i][j][k]+B_3[i][j][k])/(A_3[i][j][k]*(C+1)); | |
2180 | else term1 = 0; | |
2181 | term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1)); | |
2182 | ||
2183 | if(term1 > 0.0 && term2 > 0.0){ | |
2184 | lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2); | |
2185 | }else if(term1==0 && term2 > 0.0){ | |
2186 | lnL += B_3[i][j][k]*log(term2); | |
2187 | }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;} | |
2188 | */ | |
b2c0d4ba | 2189 | float DegenerateCount=1; |
2190 | if(i==j && i==k) DegenerateCount=1; | |
2191 | else if(i==j || i==k || j==k) DegenerateCount=3; | |
2192 | else DegenerateCount=6; | |
2193 | NFitPoints_c3 += 1/DegenerateCount; | |
2194 | ||
9e040a65 | 2195 | } |
2196 | } | |
2197 | } | |
2198 | //f = -2.0*lnL;// log-liklihood minimization | |
2199 | f = SumChi2;// Chi2 minimization | |
2200 | Chi2_c3 = f; | |
0e610a08 | 2201 | |
9e040a65 | 2202 | } |
2203 | //________________________________________________________________________ | |
2204 | double Dfitfunction_c3(Double_t *x, Double_t *par) | |
2205 | { | |
2206 | double Rch = par[2]/FmToGeV; | |
2207 | double arg12 = x[0]*Rch; | |
2208 | double arg13 = x[1]*Rch; | |
2209 | double arg23 = x[2]*Rch; | |
0e610a08 | 2210 | double Expan12 = 1, Expan13 = 1, Expan23 = 1; |
2211 | if(Gaussian){// Edworth Expansion | |
2212 | Expan12 += par[3]/(6.*pow(2,1.5))*(8*pow(arg12,3) - 12*pow(arg12,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg12,4) -48*pow(arg12,2) + 12); | |
2213 | Expan13 += par[3]/(6.*pow(2,1.5))*(8*pow(arg13,3) - 12*pow(arg13,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg13,4) -48*pow(arg13,2) + 12); | |
2214 | Expan23 += par[3]/(6.*pow(2,1.5))*(8*pow(arg23,3) - 12*pow(arg23,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg23,4) -48*pow(arg23,2) + 12); | |
2215 | }else{// Laguerre Expansion | |
2216 | Expan12 = 1 + par[3]*(arg12 - 1) + par[4]/2.*(pow(arg12,2) - 4*arg12 + 2); | |
2217 | Expan13 = 1 + par[3]*(arg13 - 1) + par[4]/2.*(pow(arg13,2) - 4*arg13 + 2); | |
2218 | Expan23 = 1 + par[3]*(arg23 - 1) + par[4]/2.*(pow(arg23,2) - 4*arg23 + 2); | |
2219 | } | |
2220 | ||
9e040a65 | 2221 | // |
0e610a08 | 2222 | double C = 1 + par[1]*exp(-(pow(arg12,ExpPower)+pow(arg13,ExpPower)+pow(arg23,ExpPower))/2.)*Expan12*Expan13*Expan23; |
9e040a65 | 2223 | C *= par[0];// Norm |
2224 | ||
2225 | return C; | |
2226 | } | |
2227 | //________________________________________________________________________ | |
2228 | double CoulCorrGRS(bool SC, double Q_12, double Q_13, double Q_23){ | |
2229 | ||
2230 | int index12L = int(fabs(Q_12 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1))); | |
2231 | int index12H = index12L+1; | |
2232 | int index13L = int(fabs(Q_13 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1))); | |
2233 | int index13H = index13L+1; | |
2234 | int index23L = int(fabs(Q_23 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1))); | |
2235 | int index23H = index23L+1; | |
2236 | ||
2237 | if(Tricubic){ | |
2238 | // Tricubic Interpolation | |
2239 | double arr[4][4][4]={{{0}}}; | |
2240 | for(int x=0; x<4; x++){ | |
2241 | for(int y=0; y<4; y++){ | |
2242 | for(int z=0; z<4; z++){ | |
2243 | if(SC){ | |
2244 | arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2SS->GetBinContent((index23L)+y)*CoulCorr2SS->GetBinContent((index13L)+z); | |
2245 | }else{ | |
2246 | arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2OS->GetBinContent((index23L)+y)*CoulCorr2OS->GetBinContent((index13L)+z); | |
2247 | } | |
2248 | ||
2249 | } | |
2250 | } | |
2251 | } | |
2252 | return tricubicInterpolate(arr, Q_12, Q_23, Q_13); | |
2253 | }else{ | |
2254 | // Trilinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation | |
2255 | // | |
2256 | double value1=1.0, value2=1.0; | |
2257 | // | |
2258 | double xd = (Q_12-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1)); | |
2259 | xd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1)); | |
2260 | double yd = (Q_13-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1)); | |
2261 | yd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index13H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1)); | |
2262 | double zd = (Q_23-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1)); | |
2263 | zd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index23H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1)); | |
2264 | // | |
2265 | if(SC){ | |
2266 | double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd); | |
2267 | c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd; | |
2268 | double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd); | |
2269 | c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd; | |
2270 | double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd); | |
2271 | c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd; | |
2272 | double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd); | |
2273 | c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd; | |
2274 | // | |
2275 | double c0 = c00*(1-yd) + c10*yd; | |
2276 | double c1 = c01*(1-yd) + c11*yd; | |
2277 | value1 = (c0*(1-zd) + c1*zd); | |
2278 | }else{ | |
2279 | double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd); | |
2280 | c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd; | |
2281 | double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd); | |
2282 | c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd; | |
2283 | double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd); | |
2284 | c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd; | |
2285 | double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd); | |
2286 | c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd; | |
2287 | // | |
2288 | double c0 = c00*(1-yd) + c10*yd; | |
2289 | double c1 = c01*(1-yd) + c11*yd; | |
2290 | value1 = (c0*(1-zd) + c1*zd); | |
2291 | } | |
2292 | ||
0e610a08 | 2293 | //if(Mbin_def<=9 || Mbin_def>12) return value1;// old way for v5 |
9e040a65 | 2294 | // |
2295 | // interpolation for K factor tansition | |
2296 | if(SC){ | |
2297 | double c00 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*(1-xd); | |
2298 | c00 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*xd; | |
2299 | double c10 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*(1-xd); | |
2300 | c10 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*xd; | |
2301 | double c01 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*(1-xd); | |
2302 | c01 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*xd; | |
2303 | double c11 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*(1-xd); | |
2304 | c11 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*xd; | |
2305 | // | |
2306 | double c0 = c00*(1-yd) + c10*yd; | |
2307 | double c1 = c01*(1-yd) + c11*yd; | |
2308 | value2 = (c0*(1-zd) + c1*zd); | |
2309 | }else{ | |
2310 | double c00 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*(1-xd); | |
2311 | c00 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*xd; | |
2312 | double c10 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*(1-xd); | |
2313 | c10 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*xd; | |
2314 | double c01 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*(1-xd); | |
2315 | c01 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*xd; | |
2316 | double c11 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*(1-xd); | |
2317 | c11 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*xd; | |
2318 | // | |
2319 | double c0 = c00*(1-yd) + c10*yd; | |
2320 | double c1 = c01*(1-yd) + c11*yd; | |
2321 | value2 = (c0*(1-zd) + c1*zd); | |
2322 | } | |
2323 | ||
2324 | if(value1>0 && value2>0){ | |
0e610a08 | 2325 | return (value1 + (value2-value1)/(2.));// always take the half-way point |
2326 | //return (value1 + (Mbin_def-9)*(value2-value1)/(3.));// old way for v5 | |
9e040a65 | 2327 | }else if(value1>0){ |
2328 | return value1; | |
2329 | }else if(value2>0){ | |
2330 | return value2; | |
2331 | }else return 1.0; | |
2332 | ||
2333 | } | |
2334 | ||
2335 | } | |
2336 | ||
2337 | double Gamov(int chargeProduct, double qinv){ | |
2338 | ||
2339 | double arg = chargeProduct*2.*PI/(BohrR*qinv/2.); | |
2340 | ||
2341 | return arg/(exp(arg)-1); | |
2342 | } | |
2343 | void ReadMomResFile(int Mbin){ | |
2344 | ||
2345 | int MRCindex_L=0, MRCindex_H=0; | |
2346 | float MRCindex_weight=0; | |
2347 | if(Mbin<=2) {MRCindex_L=0; MRCindex_H=1; MRCindex_weight = Mbin/2.;} | |
2348 | else if(Mbin<=6) {MRCindex_L=1; MRCindex_H=2; MRCindex_weight = fabs(Mbin-3)/3.;} | |
2349 | else if(Mbin<=11) {MRCindex_L=2; MRCindex_H=3; MRCindex_weight = fabs(Mbin-7)/4.;} | |
2350 | else {MRCindex_L=3; MRCindex_H=4; MRCindex_weight = fabs(Mbin-12)/7.;} | |
2351 | ||
2352 | ||
2353 | TH1D *temp_L[2][5]; | |
2354 | TH1D *temp_H[2][5]; | |
2355 | TH1D *tempC2_L[2]; | |
2356 | TH1D *tempC2_H[2]; | |
2357 | TString *momresfilenameL = new TString("MomResFile_M"); | |
2358 | *momresfilenameL += MRCindex_L; | |
2359 | momresfilenameL->Append(".root"); | |
2360 | TFile *MomResFileL = new TFile(momresfilenameL->Data(),"READ"); | |
2361 | TString *names1D_L[2][5];// SC/MC, term# | |
2362 | for(int ChProd=0; ChProd<2; ChProd++){ | |
2363 | for(int term=0; term<5; term++){ | |
2364 | // | |
2365 | if(ChProd==0) {names1D_L[ChProd][term] = new TString("MRC3_SC_term");} | |
2366 | else {names1D_L[ChProd][term] = new TString("MRC3_MC_term");} | |
2367 | *names1D_L[ChProd][term] += term+1; | |
2368 | temp_L[ChProd][term] = (TH1D*)MomResFileL->Get(names1D_L[ChProd][term]->Data()); | |
2369 | temp_L[ChProd][term]->SetDirectory(0); | |
2370 | ||
2371 | } | |
2372 | TString *C2MRCname=new TString("MomResHisto_"); | |
2373 | if(ChProd==0) C2MRCname->Append("pp"); | |
2374 | else C2MRCname->Append("mp"); | |
2375 | tempC2_L[ChProd]=(TH1D*)(((TH2D*)(MomResFileL->Get(C2MRCname->Data())))->ProjectionY("C2MRCproL",MRC2index,MRC2index)); | |
2376 | tempC2_L[ChProd]->SetDirectory(0); | |
2377 | } | |
2378 | ||
2379 | // | |
2380 | TString *momresfilenameH = new TString("MomResFile_M"); | |
2381 | *momresfilenameH += MRCindex_H; | |
2382 | momresfilenameH->Append(".root"); | |
2383 | TFile *MomResFileH = new TFile(momresfilenameH->Data(),"READ"); | |
2384 | TString *names1D_H[2][5];// SC/MC, term# | |
2385 | for(int ChProd=0; ChProd<2; ChProd++){ | |
2386 | for(int term=0; term<5; term++){ | |
2387 | // | |
2388 | if(ChProd==0) {names1D_H[ChProd][term] = new TString("MRC3_SC_term");} | |
2389 | else {names1D_H[ChProd][term] = new TString("MRC3_MC_term");} | |
2390 | *names1D_H[ChProd][term] += term+1; | |
2391 | temp_H[ChProd][term] = (TH1D*)MomResFileH->Get(names1D_H[ChProd][term]->Data()); | |
2392 | temp_H[ChProd][term]->SetDirectory(0); | |
2393 | MomRes1d[ChProd][term] = (TH1D*)MomResFileH->Get(names1D_H[ChProd][term]->Data()); | |
2394 | MomRes1d[ChProd][term]->SetDirectory(0); | |
2395 | } | |
2396 | TString *C2MRCname=new TString("MomResHisto_"); | |
2397 | if(ChProd==0) C2MRCname->Append("pp"); | |
2398 | else C2MRCname->Append("mp"); | |
2399 | tempC2_H[ChProd]=(TH1D*)(((TH2D*)(MomResFileH->Get(C2MRCname->Data())))->ProjectionY("C2MRCproH",MRC2index,MRC2index)); | |
2400 | tempC2_H[ChProd]->SetDirectory(0); | |
2401 | TString *name=new TString("MomResC2_"); | |
2402 | *name += ChProd; | |
2403 | MomResC2[ChProd] = (TH1D*)(((TH2D*)(MomResFileH->Get(C2MRCname->Data())))->ProjectionY(name->Data(),MRC2index,MRC2index)); | |
2404 | MomResC2[ChProd]->SetDirectory(0); | |
2405 | } | |
2406 | ||
2407 | for(int ChProd=0; ChProd<2; ChProd++){ | |
2408 | // C3 MRC | |
2409 | for(int term=0; term<5; term++){ | |
2410 | for(int bin=1; bin<=temp_H[ChProd][term]->GetNbinsX(); bin++){ | |
2411 | double value=1; | |
2412 | if(temp_L[ChProd][term]->GetBinContent(bin)>0 && temp_H[ChProd][term]->GetBinContent(bin)>0){// both have entries | |
2413 | value = temp_L[ChProd][term]->GetBinContent(bin) + MRCindex_weight * temp_H[ChProd][term]->GetBinContent(bin); | |
2414 | }else if(temp_L[ChProd][term]->GetBinContent(bin)>0){ | |
2415 | value = temp_L[ChProd][term]->GetBinContent(bin); | |
2416 | }else if(temp_H[ChProd][term]->GetBinContent(bin)>0){ | |
2417 | value = temp_H[ChProd][term]->GetBinContent(bin); | |
2418 | }else value=1.0; | |
2419 | ||
2420 | MomRes1d[ChProd][term]->SetBinContent(bin,value); | |
2421 | } | |
2422 | } | |
2423 | // C2 MRC | |
2424 | for(int bin=1; bin<=tempC2_H[ChProd]->GetNbinsX(); bin++){ | |
2425 | double value=1; | |
2426 | if(tempC2_L[ChProd]->GetBinContent(bin)>0 && tempC2_H[ChProd]->GetBinContent(bin)>0){// both have entries | |
2427 | value = tempC2_L[ChProd]->GetBinContent(bin)*(1-MRCindex_weight) + MRCindex_weight * tempC2_H[ChProd]->GetBinContent(bin); | |
2428 | }else if(tempC2_L[ChProd]->GetBinContent(bin)>0){ | |
2429 | value = tempC2_L[ChProd]->GetBinContent(bin); | |
2430 | }else if(tempC2_H[ChProd]->GetBinContent(bin)>0){ | |
2431 | value = tempC2_H[ChProd]->GetBinContent(bin); | |
2432 | }else value=1.0; | |
2433 | MomResC2[ChProd]->SetBinContent(bin,value); | |
2434 | } | |
2435 | } | |
2436 | ||
2437 | for(int ChProd=0; ChProd<2; ChProd++){ | |
2438 | for(int term=0; term<5; term++){ | |
2439 | for(int i=0; i<MomRes1d[0][0]->GetNbinsX(); i++){ | |
2440 | if(SourceType==0 && Mbin<10){ | |
2441 | if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.1) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0); | |
2442 | }else if(SourceType==0 && Mbin>=10){ | |
2443 | if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.2) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0); | |
2444 | }else{ | |
2445 | if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.5) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0); | |
2446 | } | |
2447 | ||
2448 | } | |
2449 | } | |
2450 | ||
2451 | } | |
2452 | ||
2453 | MomResFileL->Close(); | |
2454 | MomResFileH->Close(); | |
2455 | ||
2456 | } | |
2457 | double cubicInterpolate (double p[4], double x) { | |
2458 | return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));// Paulinternet | |
2459 | } | |
2460 | double bicubicInterpolate (double p[4][4], double x, double y) { | |
2461 | double arr[4]; | |
2462 | arr[0] = cubicInterpolate(p[0], y); | |
2463 | arr[1] = cubicInterpolate(p[1], y); | |
2464 | arr[2] = cubicInterpolate(p[2], y); | |
2465 | arr[3] = cubicInterpolate(p[3], y); | |
2466 | return cubicInterpolate(arr, x); | |
2467 | } | |
2468 | double tricubicInterpolate (double p[4][4][4], double x, double y, double z) { | |
2469 | double arr[4]; | |
2470 | arr[0] = bicubicInterpolate(p[0], y, z); | |
2471 | arr[1] = bicubicInterpolate(p[1], y, z); | |
2472 | arr[2] = bicubicInterpolate(p[2], y, z); | |
2473 | arr[3] = bicubicInterpolate(p[3], y, z); | |
2474 | return cubicInterpolate(arr, x); | |
2475 | } | |
5130881c | 2476 | //____________________________________________________________________________________________________ |
5be394f6 | 2477 | void ReadMuonCorrections(int ct, int mbin){ |
2478 | TString *name = new TString("MuonCorrection_"); | |
2479 | if(ct==0){// PbPb | |
2480 | if(mbin<=1) *name += 0; | |
2481 | else if(mbin<=3) *name += 1; | |
2482 | else if(mbin<=5) *name += 2; | |
2483 | else if(mbin<=7) *name += 3; | |
2484 | else if(mbin<10) *name += 4; | |
2485 | else *name += 5; | |
2486 | }else *name += 6; | |
2487 | name->Append(".root"); | |
2488 | TFile *muonfile=new TFile(name->Data(),"READ"); | |
2489 | C2muonCorrection_SC = (TH1D*)muonfile->Get("C2muonCorrection_SC"); | |
2490 | C2muonCorrection_MC = (TH1D*)muonfile->Get("C2muonCorrection_MC"); | |
2491 | if(SameCharge_def) C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_SC"); | |
2492 | else C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_MC"); | |
2493 | // | |
2494 | C2muonCorrection_SC->SetDirectory(0); | |
2495 | C2muonCorrection_MC->SetDirectory(0); | |
2496 | C3muonCorrection->SetDirectory(0); | |
2497 | // | |
2498 | muonfile->Close(); | |
2499 | } | |
2500 | //____________________________________________________________________________________________________ | |
5130881c | 2501 | void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2) |
2502 | { | |
2503 | // revision on July 24th or 25th 2012 | |
2504 | // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!) | |
2505 | x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw()); | |
2506 | // Printf("%f %f %f %f", x1, x2, y1, y2); | |
2507 | ||
2508 | TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2); | |
2509 | myPadLogo->SetLeftMargin(0); | |
2510 | myPadLogo->SetTopMargin(0); | |
2511 | myPadLogo->SetRightMargin(0); | |
2512 | myPadLogo->SetBottomMargin(0); | |
2513 | myPadLogo->Draw(); | |
2514 | myPadLogo->cd(); | |
2515 | TASImage *myAliceLogo = new TASImage((prel) ? "~/Pictures/2011-Nov-24-ALICE_PRELIMARY_logo_BLACK_small_usage_design.eps" : "~/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.eps"); | |
2516 | myAliceLogo->Draw(); | |
2517 | } |