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