17 #include "TProfile2D.h"
27 #include "TPaveStats.h"
29 #define BohrR 1963.6885 // Bohr Radius for pions
30 #define FmToGeV 0.19733 // conversion to fm
32 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
36 int FileSetting=0;// 0(standard), 1(r3 Lambda), 2(TTC), 3(PID), 4(Pos B), 5(Neg B), 6(Gauss), 7(4 kt bins old TTC cuts)
37 bool MCcase_def=kFALSE;// MC data?
38 int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case
39 bool SameCharge_def=kTRUE;// 3-pion same-charge?
40 int EDbin_def=0;// 0 or 1: Kt3 bin
42 bool IncludeG_def=kFALSE;// Include coherence?
43 bool IncludeEW_def=kTRUE;// Include EdgeWorth coefficients?
44 bool IncludeEWfromTherm_def=kFALSE;// Include EdgeWorth coefficients from Therminator?
45 bool GRS_def=kTRUE;// Generalized RiverSide 3-body FSI correction?
46 int Mbin_def=0;// 0-9: centrality bin in widths of 5%
47 int Ktbin_def=10;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
48 double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase
51 bool IncludeMJcorrection=kTRUE;// linear Mini-Jet correction for denominator of r3?
52 bool SaveToFile_def=kFALSE;// Save outputs to file?
53 int SourceType=1;// 0=Gaussian, 1=Therminator, 2=Lorentzian (keep at 1 for default)
54 bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin?
55 bool GofP=kFALSE;// Include momentum dependence of coherent fraction?
56 bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states?
57 bool LinkN=kFALSE;// Make N(++)=N(+-)?
58 bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
59 bool Tricubic=kFALSE;// Tricubic or Trilinear interpolation? Tricubic was found to be unstable.
60 bool IncludeSS=kTRUE;// Include same-charge two-pion correlations in the fit?
61 bool IncludeOS=kTRUE;// Include mixed-charge two-pion correlations in the fit?
66 const int Sbins_2=1;// 2-particle Species bins. 1=Pion-Pion only. max=6 (pi-pi, pi-k, pi-p, k-p, k-k, p-p)
67 const int Sbins_3=1;// 3-particle Species bins. 1=Pion-Pion-Pion only. max=10
68 const int fitlimitLowBin = 3;
69 const int fitlimitHighBin = 14;// max = 20 (100MeV)
70 bool LHC11h=kTRUE;// always kTRUE
71 bool ExplicitLoop=kFALSE;// always kFALSE
72 bool PairCut=kTRUE;// always kTRUE
73 bool PbPbcase=kTRUE;// always kTRUE
77 int bBin;// impact parameter bin
78 int MomResCentBin;// momentum resolution centrality bin
79 int RValue;// Radius setting for Gaussian source profile
80 int bValue;// impact parameter for Therminator source profile (default)
81 float TwoFrac;// Lambda parameter
82 float OneFrac;// Lambda^{1/2}
83 float ThreeFrac;// Lambda^{3/2}
84 float TwoFracMomRes;// Lambda parameter for momentum resolution corrections
85 float RValueMomRes;// Gaussian radius value for momentum resolution corrections
86 float TwoFracr3;// Lambda parameter for r3
87 int TPN_bin;// Two Particle Normalization bin (always 0)
88 double Qcut_pp = 0.6;// 0.6
89 double Qcut_PbPb = 0.1;
90 double NormQcutLow_pp = 1.0;
91 double NormQcutHigh_pp = 1.5;
92 double NormQcutLow_PbPb = .15;//1.06
93 double NormQcutHigh_PbPb = .175;//1.1
95 // single-pion emitter size (for fits with momentum dependence of coherence)
96 double RT = 0.72/FmToGeV;
97 double Temp = 0.135;// 0.135 GeV
100 const int Nlines = 50;
101 TH3D *CoulCorrOmega0SS;
102 TH3D *CoulCorrOmega0OS;
106 //static double Lednicky_qinv[74];
107 //static double Lednicky_CoulStrong[74];
109 void ReadCoulCorrections(int, int, int, int);
110 void ReadCoulCorrections_Omega0();
111 //void ReadLednickyFile(int);
112 void ReadMomResFile(int, double);
113 double CoulCorr2(int, double);
114 double CoulCorrOmega0(bool, double, double, double);
115 double CoulCorrGRS(bool, double, double, double);
116 double OSLfitfunction(Double_t *x, Double_t *par);
117 double D3fitfunction_3(Double_t *x, Double_t *par);
118 double r3_fitfunction(Double_t *x, Double_t *par);
119 double Gamov(int, double);
120 //double LednickyCorr(double);
121 double C2ssFitFunction(Double_t *x, Double_t *par);
122 double C2osFitFunction(Double_t *x, Double_t *par);
123 double cubicInterpolate(double[4], double);
124 double bicubicInterpolate(double[4][4], double, double);
125 double tricubicInterpolate(double[4][4][4], double, double, double);
127 void fcnC2_Global(int&, double*, double&, double[], int);
128 void fcnOSL(int&, double*, double&, double[], int);
130 const int BINRANGE_C2global=20;
131 double C2ssFitting[BINRANGE_C2global];
132 double C2osFitting[BINRANGE_C2global];
133 double C2ssFitting_e[BINRANGE_C2global];
134 double C2osFitting_e[BINRANGE_C2global];
135 double K2SS[BINRANGE_C2global];
136 double K2OS[BINRANGE_C2global];
137 double K2SS_e[BINRANGE_C2global];
138 double K2OS_e[BINRANGE_C2global];
139 double K3SS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
140 double K3OS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
141 double Chi2_C2global;
142 double NFitPoints_C2global;
143 double Chi2_C3global;
144 double NFitPoints_C3global;
146 const int BINS_OSL=40;// out,side,long bins
147 double A[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
148 double B[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
149 double A_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
150 double B_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
151 double avg_q[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
152 double K_OSL[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
153 double K_OSL_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
155 double NFitPoints_OSL;
157 const int BINRANGE_3=50;// q12,q13,q23
158 double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
159 double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
161 double C3_base[BINRANGE_3][BINRANGE_3][BINRANGE_3];
162 double N_term1[BINRANGE_3][BINRANGE_3][BINRANGE_3];
163 double N_term2[BINRANGE_3][BINRANGE_3][BINRANGE_3];
164 double N_term3[BINRANGE_3][BINRANGE_3][BINRANGE_3];
165 double N_term4[BINRANGE_3][BINRANGE_3][BINRANGE_3];
166 double N_term5[BINRANGE_3][BINRANGE_3][BINRANGE_3];
168 double MomRes_C2_pp[40];// Qinv bins
169 double MomRes_C2_mp[40];//
170 double MomRes_term1_pp[40];
171 double MomRes_term2_pp[40];
172 double MomRes_term1_mp[40];
173 double MomRes_term2_mp[40];
175 TH3D *MomRes3d[2][5];// SC/MC, term#
176 TH1D *MomRes1d[2][5];// SC/MC, term#
177 TH3D *MomRes3d_FVP[2][5];// sum#, term#
178 TH3D *AvgK3[2];// SC/MC
179 double AvgP[6][20];// kt bin, qinv bin
182 double MomResDev_C2_pp[40];// Qinv bins
183 double MomResDev_C2_mp[40];
184 double MomResDev_term1_pp[40];
185 double MomResDev_term2_pp[40];
188 double r3_3d_arr[20][20][20];
189 double r3_3d_arr_e[20][20][20];
190 double AvgQinvSS[30];
191 double AvgQinvOS[30];
195 void Plot_PDCumulants(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool IncludeEWfromTherm=IncludeEWfromTherm_def, bool SameCharge=SameCharge_def, bool IncludeG=IncludeG_def, bool IncludeEW=IncludeEW_def, bool GRS=GRS_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
198 SaveToFile_def=SaveToFile;
200 IncludeEWfromTherm_def=IncludeEWfromTherm;
202 IncludeG_def=IncludeG;
203 IncludeEW_def=IncludeEW;
204 SameCharge_def=SameCharge;// 3-pion same-charge
207 if(Ktbin_def==10) Ktbin_GofP=2;
208 else Ktbin_GofP=Ktbin_def;
210 Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
211 Kthighbin=(Ktbin)*2+4;
214 if(ConstantFSI) KtbinFSI=10;
218 TwoFracr3 = 0.7;// standard lambda is 0.7 for r3
219 if(FileSetting==1) TwoFracr3 = 0.6;
221 if(Mbin == 0) {// 0-5%
225 RValueMomRes=10;// for C2
228 MomResCentBin=1;// for C3
229 }else if(Mbin == 1) {// 5-10%
237 }else if(Mbin <= 3) {// 10-20%
245 }else if(Mbin <= 5) {// 20-30%
253 }else if(Mbin <= 7){// 30-40%
271 OneFrac = sqrt(TwoFrac);
272 ThreeFrac = pow(TwoFrac,3/2.);
275 if(SourceType==0 && RValue > 8) {cout<<"Radius value too large!!!"<<endl; return;}
277 cout<<"Mbin = "<<Mbin<<" Kt = "<<Ktbin<<" R input = "<<RValue<<" lambda input = "<<TwoFrac<<endl;
279 // Core-Halo modeling of single-particle and triple-particle core fraction
280 //float AvgN[10]={472.56, 390.824, 319.597, 265.933, 218.308, 178.865, 141.2, 115.88, 87.5556, 69.3235};// (avg total pion mult)/2. 2.0sigma
282 float fc = (1+sqrt(1 + 4*AvgN[Mbin]*TwoFrac*(AvgN[Mbin]-1)))/(2*AvgN[Mbin]);// Two component model (core + non-interacting halo)
283 cout<<"Before: "<<OneFrac<<" "<<ThreeFrac<<endl;
285 ThreeFrac = (fc*AvgN[Mbin]*(fc*AvgN[Mbin]-1)*(fc*AvgN[Mbin]-2))/(AvgN[Mbin]*(AvgN[Mbin]-1)*(AvgN[Mbin]-2));
286 cout<<"After: "<<OneFrac<<" "<<ThreeFrac<<" "<<endl;
289 // avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values. From unsmeared HIJING 0-5% with input signal
290 double AvgQinvSS_temp[30]={0.00421164, 0.00796583, 0.0128473, 0.0178262, 0.0228416, 0.0276507, 0.0325368, 0.0376114, 0.0424707, 0.0476274, 0.0526015, 0.0575645, 0.0625478, 0.0675416, 0.0725359, 0.0775219, 0.0825521, 0.0875211, 0.0925303, 0.0975319, 0.102544, 0.10753, 0.112499, 0.117545, 0.122522, 0.127522, 0.132499, 0.137514, 0.142495, 0.147521};
291 double AvgQinvOS_temp[30]={0.00352789, 0.00797596, 0.0128895, 0.0177198, 0.0226397, 0.0276331, 0.0326309, 0.0376471, 0.0426217, 0.047567, 0.0525659, 0.0575472, 0.0625886, 0.0675261, 0.0725543, 0.077564, 0.0825617, 0.0875465, 0.092539, 0.0975348, 0.102529, 0.107527, 0.112506, 0.117531, 0.122536, 0.1275, 0.132508, 0.137508, 0.14251, 0.147511};
292 for(int i=0; i<30; i++){
293 AvgQinvSS[i] = AvgQinvSS_temp[i];
294 AvgQinvOS[i] = AvgQinvOS_temp[i];
298 // average total momentum in each kt and qinv bin. Used for a test of momentum dependence of fits including coherence
299 double AvgP_kt1[20]={0.25, 0.254421, 0.270286, 0.27346, 0.274017, 0.274252, 0.274438, 0.274675, 0.27497, 0.275354, 0.275844, 0.276433, 0.27708, 0.277768, 0.278469, 0.279179, 0.279871, 0.280536, 0.281167, 0.281749};
300 double AvgP_kt2[20]={0.34, 0.34282, 0.370731, 0.38009, 0.381977, 0.382246, 0.382231, 0.38212, 0.381979, 0.381813, 0.381653, 0.381502, 0.381358, 0.381227, 0.381111, 0.381008, 0.380923, 0.380864, 0.380835, 0.380824};
301 double AvgP_kt3[20]={0.47, 0.47, 0.469358, 0.486368, 0.492106, 0.49335, 0.493497, 0.49341, 0.493263, 0.493048, 0.492821, 0.492587, 0.49234, 0.492093, 0.491842, 0.491607, 0.491363, 0.491126, 0.490887, 0.490645};
302 double AvgP_kt4[20]={0.59, 0.59, 0.59, 0.588732, 0.598683, 0.601864, 0.602459, 0.602302, 0.602005, 0.601617, 0.601192, 0.600794, 0.600365, 0.599961, 0.599571, 0.599208, 0.598841, 0.598511, 0.59816, 0.597836};
303 double AvgP_kt5[20]={0.67, 0.67, 0.67, 0.676098, 0.694986, 0.701248, 0.703524, 0.704307, 0.704402, 0.704432, 0.704255, 0.703942, 0.703521, 0.703031, 0.702484, 0.701869, 0.701206, 0.700521, 0.699822, 0.699084};
304 double AvgP_kt6[20]={0.79, 0.79, 0.79, 0.79, 0.790909, 0.799252, 0.801984, 0.801945, 0.800688, 0.799177, 0.797658, 0.796156, 0.794652, 0.793256, 0.791942, 0.790697, 0.789588, 0.788608, 0.787698, 0.786941};
305 for(int ktbin=1; ktbin<=6; ktbin++){
306 for(int qbin=1; qbin<=20; qbin++){
307 if(ktbin==1) AvgP[ktbin-1][qbin-1] = AvgP_kt1[qbin-1];
308 else if(ktbin==2) AvgP[ktbin-1][qbin-1] = AvgP_kt2[qbin-1];
309 else if(ktbin==3) AvgP[ktbin-1][qbin-1] = AvgP_kt3[qbin-1];
310 else if(ktbin==4) AvgP[ktbin-1][qbin-1] = AvgP_kt4[qbin-1];
311 else if(ktbin==5) AvgP[ktbin-1][qbin-1] = AvgP_kt5[qbin-1];
312 else AvgP[ktbin-1][qbin-1] = AvgP_kt6[qbin-1];
318 gStyle->SetOptStat(0);
319 gStyle->SetOptDate(0);
320 gStyle->SetOptFit(0111);
326 myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11.root","READ");
327 //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11_Smeared.root","READ");
329 myfile = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p68R11_C2plots.root","READ");
332 //myfile = new TFile("Results/PDC_11h_v1paper.root","READ");// v1 paper
333 //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// Old standard (bad r3 Interp, 0.035TTC)
334 //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// v3 Standard
336 //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins_FullTPC.root","READ");// FullTPC runlist
337 if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// Standard
338 if(FileSetting==1) myfile = new TFile("Results/PDC_11h_Lam0p7_and_Lam0p6.root","READ");// Lam=0.7 and Lam=0.6 (3ktbins, 0.035TTC)
339 if(FileSetting==2) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// TTC variation
340 if(FileSetting==3) myfile = new TFile("Results/PDC_11h_PID1p5.root","READ");// PID variation (0.035TTC)
341 if(FileSetting==4) myfile = new TFile("Results/PDC_11h_PosB.root","READ");// Positive B field for r3num (0.035TTC)
342 if(FileSetting==5) myfile = new TFile("Results/PDC_11h_NegB.root","READ");// Negative B field for r3num (0.035TTC)
343 if(FileSetting==6) myfile = new TFile("Results/PDC_11h_3sigmaTTCDenextrap_GaussDenextrap.root","READ");// Gaussian, 3sigmaTTC)
344 if(FileSetting==7) myfile = new TFile("Results/PDC_11h_4ktbins.root","READ");// 4 kt bins (0.035TTC)
345 if(FileSetting==8) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// 2 Kt3 bins
348 cout<<"pp not currently supported"<<endl;
352 ReadCoulCorrections(SourceType, RValue, bValue, KtbinFSI);
353 //ReadLednickyFile(RValue);
354 ReadMomResFile(RValueMomRes, TwoFracMomRes);
355 ReadCoulCorrections_Omega0();
357 /////////////////////////////////////////////////////////
363 NormQcutLow = NormQcutLow_PbPb;
364 NormQcutHigh = NormQcutHigh_PbPb;
366 NormQcutLow = NormQcutLow_pp;
367 NormQcutHigh = NormQcutHigh_pp;
373 TDirectoryFile *tdir = (TDirectoryFile*)myfile->Get("PWGCF.outputChaoticityAnalysis.root");
374 if(FileSetting!=1 && FileSetting !=6) MyList=(TList*)tdir->Get("ChaoticityOutput_1");
375 else MyList=(TList*)tdir->Get("ChaoticityOutput_2");
377 MyList=(TList*)myfile->Get("MyList");
382 TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
385 cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
388 TH1D *ChiSquaredNDF = new TH1D("ChiSquaredNDF","",2,0.5,2.5);// Chi^2/NDF records
390 // Explicit Loop Histos
391 TH2D *Two_ex_2d[2][2][Sbins_2][2];
392 TH1D *Two_ex[2][2][Sbins_2][2];
395 double NormH_pc[5]={0};
397 //TH3D *PionNorm_e[2];
399 double norm_ex_2[6][2]={{0}};
401 // 3d method histograms
402 TH3D *Three_3d[2][2][2][Sbins_3][5];
403 // 4-vect Product Sum
405 ///////////////////////////////////
409 for(int c1=0; c1<2; c1++){
410 for(int c2=0; c2<2; c2++){
411 for(int sc=0; sc<Sbins_2; sc++){
412 for(int term=0; term<2; term++){
414 TString *name2_ex = new TString("Explicit2_Charge1_");
416 name2_ex->Append("_Charge2_");
418 name2_ex->Append("_SC_");
420 name2_ex->Append("_M_");
422 name2_ex->Append("_ED_");
423 *name2_ex += 0;// EDbin
424 name2_ex->Append("_Term_");
427 if(sc==0 || sc==3 || sc==5){
428 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
433 Two_ex_2d[c1][c2][sc][term] = (TH2D*)MyList->FindObject(name2_ex->Data());
434 Two_ex_2d[c1][c2][sc][term]->Sumw2();
435 TString *proname = new TString(name2_ex->Data());
436 proname->Append("_pro");
438 if(Ktbin==10) {Ktlowbin=1; Kthighbin=Two_ex_2d[c1][c2][sc][term]->GetNbinsX();}
439 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)
441 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));
442 Two_ex[c1][c2][sc][term]->Scale(norm_ex_2[sc][0]/norm_ex_2[sc][term]);
444 Two_ex[c1][c2][sc][term]->SetMarkerStyle(20);
445 Two_ex[c1][c2][sc][term]->GetXaxis()->SetTitle("q_{inv}");
446 Two_ex[c1][c2][sc][term]->GetYaxis()->SetTitle("C_{2}");
447 Two_ex[c1][c2][sc][term]->SetTitle("");
453 for(int c3=0; c3<2; c3++){
454 for(int sc=0; sc<Sbins_3; sc++){
455 for(int term=0; term<5; term++){
457 TString *name3_ex = new TString("Explicit3_Charge1_");
459 name3_ex->Append("_Charge2_");
461 name3_ex->Append("_Charge3_");
463 name3_ex->Append("_SC_");
465 name3_ex->Append("_M_");
467 name3_ex->Append("_ED_");
469 name3_ex->Append("_Term_");
473 TString *name3_pc = new TString("PairCut3_Charge1_");
475 name3_pc->Append("_Charge2_");
477 name3_pc->Append("_Charge3_");
479 name3_pc->Append("_SC_");
481 name3_pc->Append("_M_");
483 name3_pc->Append("_ED_");
485 name3_pc->Append("_Term_");
488 ///////////////////////////////////////
489 // skip degenerate histograms
490 if(sc==0 || sc==6 || sc==9){// Identical species
491 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
492 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
494 if( (c1+c2)==1) {if(c1!=0) continue;}
495 }else {}// do nothing for pi-k-p case
497 /////////////////////////////////////////
502 Three_ex[c1][c2][c3][sc][term] = (TH1D*)MyList->FindObject(name3_ex->Data());
504 Three_ex[c1][c2][c3][sc][term]->SetMarkerStyle(20);
505 Three_ex[c1][c2][c3][sc][term]->SetLineColor(2);
506 Three_ex[c1][c2][c3][sc][term]->SetMarkerColor(2);
507 Three_ex[c1][c2][c3][sc][term]->GetXaxis()->SetTitle("Q_{3}");
508 Three_ex[c1][c2][c3][sc][term]->GetYaxis()->SetTitle("C_{3}");
509 Three_ex[c1][c2][c3][sc][term]->SetTitle("");
511 TString *name = new TString(name3_ex->Data());
512 name->Append("_Norm");
513 NormH_ex[term] = ((TH1D*)MyList->FindObject(name->Data()))->Integral();
515 if(NormH_ex[term] > 0){
516 Three_ex[c1][c2][c3][sc][term]->Scale(NormH_ex[0]/NormH_ex[term]);
517 }else{ cout<<"normalization = 0. Skipping this SC."<<endl;}
522 TString *nameNorm=new TString("PairCut3_Charge1_");
524 nameNorm->Append("_Charge2_");
526 nameNorm->Append("_Charge3_");
528 nameNorm->Append("_SC_");
530 nameNorm->Append("_M_");
532 nameNorm->Append("_ED_");
533 *nameNorm += 0;// EDbin
534 nameNorm->Append("_Term_");
536 nameNorm->Append("_Norm");
537 NormH_pc[term] = ((TH1D*)MyList->FindObject(nameNorm->Data()))->Integral();
539 if(NormH_pc[term] > 0){
543 TString *name_3d = new TString(name3_pc->Data());
544 name_3d->Append("_3D");
545 Three_3d[c1][c2][c3][sc][term] = (TH3D*)MyList->FindObject(name_3d->Data());
546 Three_3d[c1][c2][c3][sc][term]->Sumw2();
547 Three_3d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]);
548 //cout<<"Scale factor = "<<NormH_pc[0]/NormH_pc[term]<<endl;
550 // 4-vect Product Sum
551 if(c1==c2 && c1==c3 && sc==0 && term==4){
553 TString *nameDenType=new TString("PairCut3_Charge1_");
555 nameDenType->Append("_Charge2_");
557 nameDenType->Append("_Charge3_");
559 nameDenType->Append("_SC_");
561 nameDenType->Append("_M_");
562 *nameDenType += Mbin;
563 nameDenType->Append("_ED_");
564 *nameDenType += EDbin;
565 nameDenType->Append("_TPN_");
568 PionNorm[c1] = (TH3D*)MyList->FindObject(nameDenType->Data());
569 PionNorm[c1]->Sumw2();
570 PionNorm[c1]->Scale(NormH_pc[0]/NormH_pc[term]);
572 TPNRejects = (TH3D*)MyList->FindObject("fTPNRejects4");
573 TPNRejects->Scale(NormH_pc[0]/NormH_pc[term]);
580 cout<<"normalization = 0. Skipping this SC. SC="<<sc<<" c1="<<c1<<" c2="<<c2<<" c3="<<c3<<endl;
596 cout<<"Done getting Histograms"<<endl;
599 TCanvas *can1 = new TCanvas("can1", "can1",11,53,700,500);
600 can1->SetHighLightColor(2);
601 can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
602 gStyle->SetOptFit(0111);
603 can1->SetFillColor(10);
604 can1->SetBorderMode(0);
605 can1->SetBorderSize(2);
608 can1->SetFrameFillColor(0);
609 can1->SetFrameBorderMode(0);
610 can1->SetFrameBorderMode(0);
612 TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
613 legend1->SetBorderSize(1);
614 legend1->SetTextSize(.04);// small .03; large .036
615 legend1->SetFillColor(0);
617 /////////////////////////////////////////////////////////
618 /////////////////////////////////////////////////////////
619 // This part for plotting track splitting/merging effects in MC data only
621 TH3F *Merge3d_num=(TH3F*)MyList->FindObject("fPairsDetaDPhiNum");
622 TH3F *Merge3d_den=(TH3F*)MyList->FindObject("fPairsDetaDPhiDen");
623 //TH3F *Merge3d_num=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_num");
624 //TH3F *Merge3d_den=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_den");
626 TH1D *Merge1d_num[10];
627 TH1D *Merge1d_den[10];
628 TString *newnamenum[10];
629 TString *newnameden[10];
630 TF1 *MergedGaus=new TF1("MergedGaus","1-[0]*exp(-pow(x/[1],2))",-0.1, 0.1);
631 MergedGaus->SetParName(0,"Amplitude");
632 MergedGaus->SetParName(1,"width");
633 MergedGaus->SetParameter(0,0.06);
634 MergedGaus->SetParameter(1,0.01);
635 MergedGaus->SetParLimits(0,0.001,0.5);
636 MergedGaus->SetParLimits(1,0.001,0.1);
638 for(int i=2; i<10; i++){
639 if(i!=5 && i!=8) continue;// 5 and 8
640 newnamenum[i]=new TString("namenum_");
642 newnameden[i]=new TString("nameden_");
645 Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionZ(newnamenum[i]->Data(),i+1,i+1,90,110);//90,110 (phi)
646 Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionZ(newnameden[i]->Data(),i+1,i+1,90,110);// (phi)
647 //Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionY(newnamenum[i]->Data(),i+1,i+1,190,410);// 190,410 (eta)
648 //Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionY(newnameden[i]->Data(),i+1,i+1,190,410);// (eta)
649 //Merge1d_num[i]->Rebin(2);
650 //Merge1d_den[i]->Rebin(2);
651 Merge1d_num[i]->Sumw2();
652 Merge1d_den[i]->Sumw2();
653 Merge1d_num[i]->SetMarkerStyle(20);
655 if(Merge1d_den[i]->Integral(1,100)<=0) continue;
656 double SF_merge = Merge1d_num[i]->Integral(1,100)/Merge1d_den[i]->Integral(1,100);// Z projection (phi)
657 //double SF_merge = Merge1d_num[i]->Integral(1,50)/Merge1d_den[i]->Integral(1,50);// Y projection (eta)
658 Merge1d_den[i]->Scale(SF_merge);
659 Merge1d_num[i]->Divide(Merge1d_den[i]);
663 Merge1d_num[i]->SetLineColor(i+1);
664 Merge1d_num[i]->SetMarkerColor(i+1);
666 Merge1d_num[i]->SetLineColor(11);
667 Merge1d_num[i]->SetMarkerColor(11);
670 Merge1d_num[i]->SetLineColor(2);
671 Merge1d_num[i]->SetMarkerColor(2);
674 Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#phi*");
675 //Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#eta");
676 Merge1d_num[i]->GetYaxis()->SetTitle("C_{2}^{HIJING}");
677 Merge1d_num[i]->GetXaxis()->SetRangeUser(-.1,.1);
678 Merge1d_num[i]->SetMinimum(.91);
679 Merge1d_num[i]->SetMaximum(1.1);
680 Merge1d_num[i]->DrawCopy();
682 //Merge1d_num[i]->Fit(MergedGaus,"IME","",-0.1,0.1);
684 Merge1d_num[i]->DrawCopy("same");
687 TString *Dname=new TString("D=0.2*");
690 legend1->AddEntry(newnamenum[i]->Data(),Dname->Data(),"lpf");
692 legend1->Draw("same");
693 gStyle->SetOptFit(111);
694 Merge1d_num[8]->Fit(MergedGaus,"IME","",-0.1,0.1);
695 MergedGaus->Draw("same");
697 /*TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsPadRowNum");// kt, shfrac, qinv
698 TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsPadRowDen");// kt, shfrac, qinv
699 PadRowDen3->Scale(PadRowNum3->Integral(1,20,1,159, 35,40)/PadRowDen3->Integral(1,20,1,159, 35,40));
700 PadRowNum3->GetYaxis()->SetRangeUser(0,0.01);
701 PadRowDen3->GetYaxis()->SetRangeUser(0,0.01);
702 TH1D *PadRowNum=(TH1D*)PadRowNum3->Project3D("z");
703 TH1D *PadRowDen=(TH1D*)PadRowDen3->Project3D("z");
704 PadRowNum->Divide(PadRowDen);
707 TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiNum");// r, shfrac, deltaphi
708 TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiDen");// r, shfrac, deltaphi
709 PadRowDen3->Scale(PadRowNum3->Integral(1,10,1,159, 90,100)/PadRowDen3->Integral(1,10,1,159, 90,100));
710 PadRowNum3->GetXaxis()->SetRange(5,5);
711 PadRowDen3->GetXaxis()->SetRange(5,5);
712 TH2D *PadRowNum=(TH2D*)PadRowNum3->Project3D("zy");
713 TH2D *PadRowDen=(TH2D*)PadRowDen3->Project3D("zy");
714 PadRowNum->Divide(PadRowDen);
715 PadRowNum->Draw("lego");
717 /////////////////////////
725 /////////////////////////
726 TH1D *Two_pipi_mp = (TH1D*)Two_ex[0][1][0][0]->Clone();
727 Two_pipi_mp->Divide(Two_ex[0][1][0][1]);
729 // Normalization range counting. Just to evaluate required normalization width in qinv.
730 //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;
731 //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;
732 //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;
738 TH1D *Two_ex_clone_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
739 Two_ex_clone_mm->Divide(Two_ex[0][0][SCOI_2][1]);
740 TH1D *Two_ex_clone_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
741 Two_ex_clone_mp->Divide(Two_ex[0][1][SCOI_2][1]);
742 TH1D *Two_ex_clone_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
743 Two_ex_clone_pp->Divide(Two_ex[1][1][SCOI_2][1]);
745 // Mini-jet ++ background linear estimation.
746 TF1 *MJ_bkg_ss=new TF1("MJ_bkg_ss","pol1",0,1);
747 Two_ex_clone_mm->Fit(MJ_bkg_ss,"IMENQ","",0.2,0.4);
748 cout<<"Non-femto bkg C2(q=0.01) = "<<MJ_bkg_ss->Eval(0.01)<<endl;
750 Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}");
751 Two_ex_clone_mm->SetTitle("");
752 Two_ex_clone_mp->GetYaxis()->SetTitle("C_{2}");
753 Two_ex_clone_mm->SetMarkerColor(2);
754 Two_ex_clone_mm->SetLineColor(2);
755 Two_ex_clone_mp->SetMarkerColor(1);
756 Two_ex_clone_pp->SetMarkerColor(4);
757 Two_ex_clone_pp->SetLineColor(4);
758 Two_ex_clone_mm->GetXaxis()->SetRangeUser(0,0.1);
759 Two_ex_clone_mm->SetMinimum(0.95);
760 Two_ex_clone_mm->SetMaximum(1.4);
763 Two_ex_clone_mp->SetMarkerColor(4);
764 Two_ex_clone_mp->SetLineColor(4);
765 Two_ex_clone_mm->Add(Two_ex_clone_pp);
766 Two_ex_clone_mm->Scale(0.5);
767 Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}^{HIJING}");
768 Two_ex_clone_mm->GetYaxis()->SetTitleOffset(1.3);
769 Two_ex_clone_mm->SetMinimum(0.97);
770 Two_ex_clone_mm->SetMaximum(1.02);
771 Two_ex_clone_mm->DrawCopy();
772 //Two_ex_clone_pp->DrawCopy("same");
773 Two_ex_clone_mp->DrawCopy("same");
774 legend1->AddEntry(Two_ex_clone_mm,"same-charge","p");
775 //legend1->AddEntry(Two_ex_clone_pp,"++","p");
776 legend1->AddEntry(Two_ex_clone_mp,"mixed-charge","p");
777 legend1->Draw("same");
782 ////////////////////////////////////////////
783 // Get Therminator EdgeWorth coefficients
784 TString *ThermName = new TString("../ThermData/Therm_FSI_b");
785 *ThermName += bValue;
786 ThermName->Append(".root");
787 TFile *Thermfile = new TFile(ThermName->Data(),"READ");
788 TH2D *thermNum2d_ss=(TH2D*)Thermfile->Get("Num_CosFSI_ss");
789 TH2D *thermNumSq2d_ss=(TH2D*)Thermfile->Get("NumSq_CosFSI_ss");
790 TH2D *thermDen2d_ss=(TH2D*)Thermfile->Get("Den_ss");
791 TH2D *thermLargeR2D_ss=(TH2D*)Thermfile->Get("LargeRpairs_ss");
792 TH1D *thermNum_ss=(TH1D*)thermNum2d_ss->ProjectionY("thermNum_ss",Ktlowbin,Kthighbin);
793 TH1D *thermNumSq_ss=(TH1D*)thermNumSq2d_ss->ProjectionY("thermNumSq_ss",Ktlowbin,Kthighbin);
794 TH1D *thermDen_ss=(TH1D*)thermDen2d_ss->ProjectionY("thermDen_ss",Ktlowbin,Kthighbin);
795 TH1D *thermLargeR_ss=(TH1D*)thermLargeR2D_ss->ProjectionY("thermLargeR_ss",Ktlowbin,Kthighbin);
797 TH2D *thermNum2d_os=(TH2D*)Thermfile->Get("Num_CosFSI_os");
798 TH2D *thermNumSq2d_os=(TH2D*)Thermfile->Get("NumSq_CosFSI_os");
799 TH2D *thermDen2d_os=(TH2D*)Thermfile->Get("Den_os");
800 TH2D *thermLargeR2D_os=(TH2D*)Thermfile->Get("LargeRpairs_os");
801 TH1D *thermNum_os=(TH1D*)thermNum2d_os->ProjectionY("thermNum_os",Ktlowbin,Kthighbin);
802 TH1D *thermNumSq_os=(TH1D*)thermNumSq2d_os->ProjectionY("thermNumSq_os",Ktlowbin,Kthighbin);
803 TH1D *thermDen_os=(TH1D*)thermDen2d_os->ProjectionY("thermDen_os",Ktlowbin,Kthighbin);
804 TH1D *thermLargeR_os=(TH1D*)thermLargeR2D_os->ProjectionY("thermLargeR_os",Ktlowbin,Kthighbin);
805 TH1D *C2Therm_ss = (TH1D*)thermNum_ss->Clone();
806 TH1D *C2Therm_os = (TH1D*)thermNum_os->Clone();
807 TH1D *C2Den_ss = (TH1D*)thermDen_ss->Clone();
808 TH1D *C2Den_os = (TH1D*)thermDen_os->Clone();
809 C2Therm_ss->Add(thermLargeR_ss);
810 C2Den_ss->Add(thermLargeR_ss);
811 C2Therm_os->Add(thermLargeR_os);
812 C2Den_os->Add(thermLargeR_os);
816 C2Therm_ss->Divide(C2Den_ss);
817 C2Therm_os->Divide(C2Den_os);
820 for(int i=1; i<=thermDen_ss->GetNbinsX(); i++){
821 if(thermDen_ss->GetBinContent(i) > 0){
822 double err = thermNumSq_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i);
823 err -= pow(thermNum_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i),2);
824 err /= C2Den_ss->GetBinContent(i);
825 err = err + pow(sqrt(thermLargeR_ss->GetBinContent(i))/C2Den_ss->GetBinContent(i),2);
826 err += pow(sqrt(C2Den_ss->GetBinContent(i))*thermLargeR_ss->GetBinContent(i)/pow(C2Den_ss->GetBinContent(i),2),2);
828 C2Therm_ss->SetBinError(i, err);
830 if(thermDen_os->GetBinContent(i) > 0){
831 double err = thermNumSq_os->GetBinContent(i)/C2Den_os->GetBinContent(i);
832 err -= pow(thermNum_os->GetBinContent(i)/C2Den_os->GetBinContent(i),2);
833 err /= C2Den_os->GetBinContent(i);
834 err = err + pow(sqrt(thermLargeR_os->GetBinContent(i))/C2Den_os->GetBinContent(i),2);
835 err += pow(sqrt(C2Den_os->GetBinContent(i))*thermLargeR_os->GetBinContent(i)/pow(C2Den_os->GetBinContent(i),2),2);
837 C2Therm_os->SetBinError(i, err);
841 C2Therm_ss->SetMarkerStyle(20);
842 C2Therm_os->SetMarkerStyle(25);
843 C2Therm_ss->SetMarkerColor(2);
844 C2Therm_os->SetMarkerColor(4);
845 C2Therm_ss->SetLineColor(2);
846 C2Therm_os->SetLineColor(4);
847 C2Therm_ss->SetMarkerSize(1.5);
848 C2Therm_os->SetMarkerSize(1.5);
849 C2Therm_ss->SetDirectory(0);
850 C2Therm_os->SetDirectory(0);
851 C2Therm_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
852 C2Therm_ss->GetYaxis()->SetTitle("1+<cos(qx)>");
853 C2Therm_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
854 C2Therm_os->GetYaxis()->SetTitle("C_{2}^{+-}");
855 C2Therm_os->SetMinimum(0.98);
856 C2Therm_os->SetMaximum(1.25);
857 //C2Therm_ss->DrawCopy();
858 //C2Therm_os->DrawCopy("same");
859 C2Therm_ss->SetDirectory(0);
860 C2Therm_os->SetDirectory(0);
861 //C2Therm->Fit(GaussEW,"IME","",0.005,0.1);
866 /////////////////////////////////////////////////////
867 // Global fitting C2os and C2ss
868 double ThermEW[4]={0};
869 double C2ssSys_e[BINRANGE_C2global]={0};
870 double C2osSys_e[BINRANGE_C2global]={0};
873 TF1 *fitC2ssThermGaus;
874 TF1 *fitC2osThermGaus;
876 const int npar=11;// 10
877 TMinuit MyMinuit(npar);
878 MyMinuit.SetFCN(fcnC2_Global);
879 double OutputPar[npar]={0};
880 double OutputPar_e[npar]={0};
882 double par[npar]; // the start values
883 double stepSize[npar]; // step sizes
884 double minVal[npar]; // minimum bound on parameter
885 double maxVal[npar]; // maximum bound on parameter
886 string parName[npar];
888 TH1D *temp_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
889 temp_mm->Divide(Two_ex[0][0][SCOI_2][1]);
890 TH1D *temp_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
891 temp_mp->Divide(Two_ex[0][1][SCOI_2][1]);
892 TH1D *temp_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
893 temp_pp->Divide(Two_ex[1][1][SCOI_2][1]);
894 TH1D *C2ssRaw=(TH1D*)temp_mm->Clone();// MRC uncorrected
895 TH1D *C2osRaw=(TH1D*)temp_mp->Clone();// MRC uncorrected
896 C2ssRaw->SetMarkerStyle(24);
897 C2osRaw->SetMarkerStyle(21);//21
898 C2ssRaw->SetMarkerColor(2);
899 C2osRaw->SetMarkerColor(4);
903 while(iteration<3){// 0=pure Gaussian Therminator fits, 1=EW Therminator fits, 2=real data fits with Therm EW
904 if(!IncludeEWfromTherm) iteration=2;// skip Therminator
906 for(int i=0; i<BINRANGE_C2global; i++){
908 C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.;
909 if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= MomRes_C2_pp[i];
910 C2ssFitting_e[i] = pow(MomRes_C2_pp[i]*sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2);
911 C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.);
912 C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2));
913 C2ssFitting_e[i] += pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
914 C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]);
915 C2osFitting[i] = temp_mp->GetBinContent(i+1);
916 if(!GeneratedSignal && !MCcase) C2osFitting[i] *= MomRes_C2_mp[i];
917 C2osFitting_e[i] = pow(MomRes_C2_mp[i]*temp_mp->GetBinError(i+1),2);
918 C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1));
919 C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1));
920 C2osFitting_e[i] += pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
921 C2osFitting_e[i] = sqrt(C2osFitting_e[i]);
923 C2ssSys_e[i] = pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
924 C2ssSys_e[i] = sqrt(C2ssSys_e[i]);
925 C2osSys_e[i] = pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
926 C2osSys_e[i] = sqrt(C2osSys_e[i]);
929 K2SS[i] = CoulCorr2(+1, AvgQinvSS[i]);
930 K2OS[i] = CoulCorr2(-1, AvgQinvOS[i]);
934 K2SS_e[i] = sqrt(pow((K2SS[i]-1)*0.02,2) + pow((K2SS[i]-Gamov(+1, AvgQinvSS[i]))*0.02,2));
935 K2OS_e[i] = sqrt(pow((K2OS[i]-1)*0.02,2) + pow((K2OS[i]-Gamov(-1, AvgQinvOS[i]))*0.02,2));
941 C2ssFitting[i] = C2Therm_ss->GetBinContent(i+1);// Therminator fit
942 C2ssFitting_e[i] = C2Therm_ss->GetBinError(i+1);// Therminator fit
943 C2osFitting[i] = C2Therm_os->GetBinContent(i+1);// Therminator fit
944 C2osFitting_e[i] = C2Therm_os->GetBinError(i+1);// Therminator fit
953 C2ssFitting[0]=-1000; C2osFitting[0]=-1000;
954 C2ssFitting_e[0]=1000; C2osFitting_e[0]=1000;
955 K2SS_e[0]=1000; K2OS_e[0]=1000;
959 par[0] = 1.0; par[1]=0.5; par[2]=0.5; par[3]=9.2; par[4] = .1; par[5] = .2; par[6] = .2; par[7] = 0.; par[8] = 0.; par[9] = 0.; par[10] = 0.01;
960 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; stepSize[10]=0.01;
961 minVal[0] = 0.995; minVal[1] = 0.2; 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[10] = 0;
962 maxVal[0] = 1.1; maxVal[1] = 1.0; 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[10]=1.;
963 parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh";
964 parName[5] = "kappa_3"; parName[6] = "kappa_4"; parName[7] = "kappa_5"; parName[8]="kappa_6"; parName[9]="Norm_2"; parName[10]="P_{coh}";
966 MyMinuit.DefineParameter(10, parName[10].c_str(), 0., stepSize[10], minVal[10], maxVal[10]); MyMinuit.FixParameter(10);// extra parameter
968 for (int i=0; i<npar; i++){
969 MyMinuit.DefineParameter(i, parName[i].c_str(), par[i], stepSize[i], minVal[i], maxVal[i]);
971 //MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
972 //MyMinuit.DefineParameter(1, parName[1].c_str(), .68, stepSize[1], minVal[1], maxVal[1]); MyMinuit.FixParameter(1);// lambda
973 //MyMinuit.DefineParameter(2, parName[2].c_str(), 0.9999, stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
975 if(IncludeG==kFALSE || iteration<2) {
976 MyMinuit.DefineParameter(2, parName[2].c_str(), 0., stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
977 MyMinuit.DefineParameter(4, parName[4].c_str(), 0., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
982 MyMinuit.mnexcm( "RELease", tmp, 1, err );// G
984 MyMinuit.mnexcm( "RELease", tmp, 1, err );// Rcoh
987 //MyMinuit.DefineParameter(3, parName[3].c_str(), 10.88, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
988 //MyMinuit.DefineParameter(4, parName[4].c_str(), 5., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
989 MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
990 MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
993 if(!IncludeEWfromTherm){
994 // kappa3=0.24, kappa4=0.16 for testing
995 MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
996 MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
997 MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
998 MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1001 MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1002 MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1003 }else if(iteration==1){
1007 MyMinuit.mnexcm( "RELease", tmp, 1, err );// k3
1009 MyMinuit.mnexcm( "RELease", tmp, 1, err );// k4
1010 }else{// iteration=2
1011 MyMinuit.DefineParameter(5, parName[5].c_str(), ThermEW[0], stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1012 MyMinuit.DefineParameter(6, parName[6].c_str(), ThermEW[1], stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1013 MyMinuit.DefineParameter(7, parName[7].c_str(), ThermEW[2], stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
1014 MyMinuit.DefineParameter(8, parName[8].c_str(), ThermEW[3], stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1016 }// IncludeEWfromTherm
1019 if(IncludeSS==kFALSE){
1020 MyMinuit.DefineParameter(3, parName[3].c_str(), 9.1, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
1021 MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
1022 MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1023 MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1025 if(IncludeOS==kFALSE){
1026 MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2
1032 //arglist[0]=2;// improve Minimization Strategy (1 is default)
1033 //MyMinuit.mnexcm("SET STR",arglist,1,ierflg);
1035 //MyMinuit.mnexcm("SCAN", arglist,1,ierflg);
1037 MyMinuit.mnexcm("MIGRAD", arglist ,1,ierflg);
1038 // Do the minimization!
1039 cout<<"Start C2 Global fit"<<endl;
1040 MyMinuit.Migrad();// generally the best minimization algorithm
1041 cout<<"End C2 Global fit"<<endl;
1043 for (int i=0; i<npar; i++){
1044 MyMinuit.GetParameter(i,OutputPar[i],OutputPar_e[i]);
1047 cout<<"Chi2/NDF = "<<Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars())<<" Chi^2="<<Chi2_C2global<<endl;
1049 ChiSquaredNDF->Fill(1, Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars()));
1051 ThermEW[0]=OutputPar[5]; ThermEW[1]=OutputPar[6]; ThermEW[2]=OutputPar[7]; ThermEW[3]=OutputPar[8];
1052 fitC2ssTherm = new TF1("fitC2ssTherm",C2ssFitFunction, 0.005,0.2, npar);
1053 for(int i=0; i<npar; i++) {
1054 fitC2ssTherm->FixParameter(i,OutputPar[i]);
1055 fitC2ssTherm->SetParError(i,OutputPar_e[i]);
1057 fitC2osTherm = new TF1("fitC2osTherm",C2osFitFunction, 0.005,0.2, npar);
1058 for(int i=0; i<npar; i++) {
1059 fitC2osTherm->FixParameter(i,OutputPar[i]);
1060 fitC2osTherm->SetParError(i,OutputPar_e[i]);
1062 fitC2osTherm->SetLineColor(4);
1065 fitC2ssThermGaus = new TF1("fitC2ssThermGaus",C2ssFitFunction, 0.005,0.2, npar);
1066 for(int i=0; i<npar; i++) {
1067 fitC2ssThermGaus->FixParameter(i,OutputPar[i]);
1068 fitC2ssThermGaus->SetParError(i,OutputPar_e[i]);
1070 fitC2osThermGaus = new TF1("fitC2osThermGaus",C2osFitFunction, 0.005,0.2, npar);
1071 for(int i=0; i<npar; i++) {
1072 fitC2osThermGaus->FixParameter(i,OutputPar[i]);
1073 fitC2osThermGaus->SetParError(i,OutputPar_e[i]);
1075 fitC2osThermGaus->SetLineColor(4);
1083 TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone();
1084 TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone();
1085 TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone();
1086 TH1D *C2_os_Momsys=(TH1D*)Two_ex_clone_mp->Clone();
1087 TH1D *C2_ss_Ksys=(TH1D*)Two_ex_clone_mm->Clone();
1088 TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone();
1089 TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone();
1090 TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone();
1092 for(int i=0; i<BINRANGE_C2global; i++){
1093 C2_ss->SetBinContent(i+1, C2ssFitting[i]);
1094 C2_os->SetBinContent(i+1, C2osFitting[i]);
1095 C2_ss_Momsys->SetBinContent(i+1, C2ssFitting[i]);
1096 C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]);
1097 C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]);
1098 C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]);
1099 double staterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
1100 double staterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
1101 C2_ss->SetBinError(i+1, staterror_ss);
1102 C2_os->SetBinError(i+1, staterror_os);
1103 C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]);
1104 C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]);
1105 C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]);
1106 C2_os_Ksys->SetBinError(i+1, OutputPar[1]*K2OS_e[i]);
1108 K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0);
1109 K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0);
1112 C2_ss_Momsys->SetMarkerSize(0);
1113 C2_ss_Momsys->SetFillStyle(1000);
1114 C2_ss_Momsys->SetFillColor(kRed-10);
1115 C2_os_Momsys->SetMarkerSize(0);
1116 C2_os_Momsys->SetFillStyle(1000);
1117 C2_os_Momsys->SetFillColor(kBlue-10);
1118 C2_ss_Ksys->SetMarkerSize(0);
1119 C2_ss_Ksys->SetFillStyle(3004);
1120 C2_ss_Ksys->SetFillColor(kRed);
1121 C2_os_Ksys->SetMarkerSize(0);
1122 C2_os_Ksys->SetFillStyle(3004);
1123 C2_os_Ksys->SetFillColor(kBlue);
1127 C2_ss->GetXaxis()->SetRangeUser(0,0.098);
1128 C2_ss->GetYaxis()->SetRangeUser(0.98,1.3);
1129 C2_ss->GetXaxis()->SetTitleOffset(.8);
1130 C2_ss->GetYaxis()->SetTitleOffset(.8);
1131 C2_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1132 C2_ss->GetXaxis()->SetTitleSize(0.05);
1133 C2_ss->GetYaxis()->SetTitleSize(0.05);
1134 C2_os->GetXaxis()->SetRangeUser(0,0.098);
1135 C2_os->GetYaxis()->SetRangeUser(0.98,1.3);
1136 C2_os->GetXaxis()->SetTitleOffset(.8);
1137 C2_os->GetYaxis()->SetTitleOffset(.8);
1138 C2_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1139 C2_os->GetXaxis()->SetTitleSize(0.05);
1140 C2_os->GetYaxis()->SetTitleSize(0.05);
1142 C2_ss->SetMarkerSize(1.5);
1143 C2_os->SetMarkerSize(1.5);
1144 C2_os->SetMarkerStyle(25);
1145 C2_os->SetMarkerColor(4);
1147 TF1 *fitC2ss = new TF1("fitC2ss",C2ssFitFunction, 0.005,0.2, npar);
1148 TF1 *fitC2os = new TF1("fitC2os",C2osFitFunction, 0.005,0.2, npar);
1149 for(int i=0; i<npar; i++) {
1150 fitC2ss->FixParameter(i,OutputPar[i]);
1151 fitC2ss->SetParError(i,OutputPar_e[i]);
1153 for(int i=0; i<npar; i++) {
1154 fitC2os->FixParameter(i,OutputPar[i]);
1155 fitC2os->SetParError(i,OutputPar_e[i]);
1160 legend1->AddEntry(C2_ss,"same-charge","p");
1161 C2_os->DrawCopy("same");
1162 legend1->AddEntry(C2_os,"opp-charge","p");
1163 C2_ss_Momsys->DrawCopy("E2 same");
1164 C2_os_Momsys->DrawCopy("E2 same");
1165 C2_ss_Ksys->DrawCopy("E2 same");
1166 C2_os_Ksys->DrawCopy("E2 same");
1167 C2_ss->DrawCopy("same");
1168 C2_os->DrawCopy("same");
1169 fitC2ss->DrawCopy("same");
1170 fitC2os->SetLineColor(4);
1171 fitC2os->DrawCopy("same");
1172 MJ_bkg_ss->SetLineColor(1);
1173 MJ_bkg_ss->Draw("same");
1174 legend1->AddEntry(MJ_bkg_ss,"Non-femto bkg","l");
1175 legend1->Draw("same");
1179 ////////// C2 systematics
1180 // C2 -- base, M0, Pos B (0.035 TTC for all)
1181 //double C2ss_base[20]={0, 1.17121, 1.17272, 1.13947, 1.09819, 1.06512, 1.04115, 1.02563, 1.01558, 1.00937, 1.00551, 1.00309, 1.00157, 1.00065, 1.00023, 0.999911, 0.99976, 0.999733, 0.999593, 0.999641};
1182 //double C2ss_base_e[20]={0, 0.00135764, 0.000430769, 0.000242626, 0.000167895, 0.00012913, 0.000105724, 9.01743e-05, 7.90587e-05, 7.07036e-05, 6.41813e-05, 5.89418e-05, 5.46406e-05, 5.10458e-05, 4.80009e-05, 4.53857e-05, 4.31177e-05, 4.11375e-05, 3.93869e-05, 3.78413e-05};
1184 double C2ss_base[20]={0, 1.16926, 1.17306, 1.13932, 1.09834, 1.06475, 1.04096, 1.02546, 1.01547, 1.00931, 1.00542, 1.00311, 1.00155, 1.00071, 1.00018, 0.999786, 0.999742, 0.999612, 0.999548, 0.999631};
1185 double C2ss_base_e[20]={0, 0.00128558, 0.000410643, 0.000231512, 0.000160337, 0.000123314, 0.000100987, 8.61441e-05, 7.55245e-05, 6.75439e-05, 6.13117e-05, 5.63119e-05, 5.21976e-05, 4.87656e-05, 4.58515e-05, 4.33492e-05, 4.11861e-05, 3.92891e-05, 3.76189e-05, 3.61428e-05};
1187 //double C2ss_base[20]={1.97333, 1.17447, 1.09412, 1.05561, 1.03613, 1.02507, 1.01809, 1.01363, 1.01048, 1.00816, 1.00645, 1.00526, 1.00427, 1.0035, 1.00289, 1.00241, 1.00209, 1.0018, 1.00148, 1.00127};
1188 //double C2ss_base_e[20]={1.61122, 0.00030054, 0.000174753, 0.000123051, 9.53599e-05, 7.81537e-05, 6.64323e-05, 5.79583e-05, 5.15436e-05, 4.65277e-05, 4.25025e-05, 3.92115e-05, 3.64672e-05, 3.4149e-05, 3.21673e-05, 3.04575e-05, 2.89698e-05, 2.76642e-05, 2.65087e-05, 2.54847e-05};
1190 //double C2ss_base[20]={0, 0, 0, 1.15387, 1.14361, 1.10213, 1.06552, 1.04357, 1.02804, 1.01826, 1.01182, 1.00767, 1.00509, 1.00294, 1.00169, 1.0008, 1.00063, 1.00022, 1.00001, 0.999986};
1191 //double C2ss_base_e[20]={0, 0, 0, 0.00229526, 0.000813871, 0.000518539, 0.00039803, 0.000328676, 0.000281681, 0.00024771, 0.000221807, 0.000201394, 0.000184824, 0.000171048, 0.000159468, 0.000149554, 0.000141058, 0.000133583, 0.00012702, 0.000121218};
1193 //double C2ss_base[20]={0, 0, 0, 1.11133, 1.13824, 1.11806, 1.08288, 1.05648, 1.03774, 1.02516, 1.01728, 1.01166, 1.00841, 1.00553, 1.00357, 1.00194, 1.00132, 1.00104, 1.00058, 1.00039};
1194 //double C2ss_base_e[20]={0, 0, 0, 0.0147639, 0.00188068, 0.000910766, 0.000637117, 0.000510753, 0.000432397, 0.000377855, 0.000337624, 0.000306333, 0.000281495, 0.00026099, 0.000243881, 0.000229331, 0.000216961, 0.000206243, 0.000196785, 0.000188487};
1196 //double C2ss_base[20]={0, 0, 0, 0, 1.08106, 1.12109, 1.09838, 1.07193, 1.05087, 1.03534, 1.02433, 1.01785, 1.01339, 1.00928, 1.00662, 1.00441, 1.00326, 1.00192, 1.00175, 1.00143};
1197 //double C2ss_base_e[20]={0, 0, 0, 0, 0.00722887, 0.00211091, 0.00123591, 0.000924715, 0.000765679, 0.000662277, 0.000588573, 0.000533669, 0.000490831, 0.000456049, 0.000427965, 0.000404612, 0.000385285, 0.000368671, 0.00035474, 0.000342708};
1199 cout<<"C2 values:"<<endl;
1200 for(int ii=0; ii<20; ii++){
1201 cout<<Two_ex_clone_mm->GetBinContent(ii+1)<<", ";
1204 cout<<"C2 errors:"<<endl;
1205 for(int ii=0; ii<20; ii++){
1206 cout<<Two_ex_clone_mm->GetBinError(ii+1)<<", ";
1209 TH1D *C2ss_Sys=(TH1D*)Two_ex_clone_mm->Clone();
1210 for(int ii=1; ii<20; ii++){
1211 if(C2ss_base[ii] ==0) {
1212 C2ss_Sys->SetBinContent(ii+1, 100.);
1213 C2ss_Sys->SetBinError(ii+1, 100.);
1216 C2ss_Sys->SetBinContent(ii+1, (C2ss_base[ii]-C2ss_Sys->GetBinContent(ii+1))/C2ss_base[ii]);
1217 C2ss_Sys->SetBinError(ii+1, sqrt(fabs(pow(C2ss_Sys->GetBinError(ii+1),2) - pow(C2ss_base_e[ii],2))));
1219 C2ss_Sys->GetXaxis()->SetRangeUser(0,0.099);
1220 C2ss_Sys->GetYaxis()->SetTitle("(C_{2}^{s}-C_{2}^{v})/C_{2}^{s}");
1221 C2ss_Sys->GetYaxis()->SetTitleOffset(1.3);
1222 C2ss_Sys->SetMinimum(-0.01);
1223 C2ss_Sys->SetMaximum(0.01);
1225 TF1 *C2lineSys=new TF1("C2lineSys","pol0",0,0.099);
1226 //C2ss_Sys->Fit(C2lineSys,"MEQ","",0,0.099);
1229 // Momentum resolution uncorrected C2
1230 /*C2ssRaw->Draw("same");
1231 C2osRaw->Draw("same");
1232 legend1->AddEntry(C2ssRaw,"same-charge, Raw","p");
1233 legend1->AddEntry(C2osRaw,"opp-charge, Raw","p");
1234 legend1->Draw("same");
1237 // FSI corrected C2+-
1238 /*TH1D *C2_os_FSIcorrected=(TH1D*)C2_os->Clone();
1239 for(int ii=2; ii<20; ii++){
1240 double value = (C2_os_FSIcorrected->GetBinContent(ii) - (1-OutputPar[1]))/(OutputPar[1]*K2OS[ii-1]);
1241 C2_os_FSIcorrected->SetBinContent(ii,value);
1243 C2_os_FSIcorrected->SetMinimum(0.95);
1244 C2_os_FSIcorrected->SetMarkerStyle(24);
1245 C2_os_FSIcorrected->Draw();
1249 /*Two_ex_clone_mm->SetMarkerStyle(25);
1250 Two_ex_clone_mp->SetMarkerStyle(25);
1251 Two_ex_clone_mm->SetMarkerColor(2);
1252 Two_ex_clone_mp->SetMarkerColor(4);
1253 Two_ex_clone_mm->Draw("same");
1254 Two_ex_clone_mp->Draw("same");
1255 legend1->AddEntry(C2_ss,"same-charge; Momentum Resolution Corrected","p");
1256 legend1->AddEntry(C2_os,"opp-charge; Momentum Resolution Corrected","p");
1257 legend1->AddEntry(Two_ex_clone_mm,"same-charge; Raw","p");
1258 legend1->AddEntry(Two_ex_clone_mp,"opp-charge; Raw","p");
1259 legend1->Draw("same");
1261 //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1262 //TGraph *gr2 = (TGraph*)MyMinuit.Contour(15,1,2);
1263 //gr2->GetXaxis()->SetTitle("lambda");
1264 //gr2->GetYaxis()->SetTitle("G");
1265 //gr2->SetTitle("");
1266 //gr2->SetFillColor(kYellow);
1268 //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1269 //MyMinuit.SetErrorDef(1);
1270 //TGraph *gr1 = (TGraph*)MyMinuit.Contour(15,1,2);
1271 //gr1->SetFillColor(kGreen);//38
1276 ///////////////////////////
1278 /* double C2_star_mm[50]={0, 0, 0.566654, 1.11195, 1.09964, 1.12566, 1.1479, 1.15596, 1.15658, 1.14583, 1.13182, 1.11544, 1.09851, 1.08326, 1.0701, 1.05686, 1.048, 1.03924, 1.03122, 1.02544, 1.02019, 1.01617, 1.01239, 1.00942, 1.00736, 1.00485, 1.00296, 1.00152, 1.00039, 0.999461, 0.998703, 0.997563, 0.997294, 0.996726, 0.996739, 0.996355, 0.996426, 0.996581, 0.996336, 0.996157, 0.996312, 0.996413, 0.996501, 0.996204, 0.996816, 0.996654, 0.99695, 0.997494, 0.997135, 0.997188};
1279 double C2_star_mm_e[50]={0, 0, 0.454334, 0.0216601, 0.00590084, 0.00297232, 0.001939, 0.00142608, 0.00112851, 0.000929702, 0.000792097, 0.000690265, 0.000613054, 0.000552755, 0.000504494, 0.000464031, 0.000431371, 0.000403169, 0.000378672, 0.000357685, 0.000339133, 0.00032286, 0.000308251, 0.000295149, 0.000283433, 0.000272615, 0.000262784, 0.000253889, 0.000245693, 0.00023813, 0.000231135, 0.00022455, 0.000218593, 0.000212961, 0.000207822, 0.000202906, 0.000198339, 0.000194087, 0.000189971, 0.000186181, 0.000182575, 0.000179202, 0.000175989, 0.000172894, 0.000170072, 0.000167293, 0.000164706, 0.000162286, 0.000159867, 0.000157606};
1280 TH1D *Two_star_mm=(TH1D*)Two_ex_clone_mm->Clone();
1281 Two_star_mm->Add(Two_star_mm,-1);
1282 for(int i=0; i<50; i++) {Two_star_mm->SetBinContent(i+1, C2_star_mm[i]); Two_star_mm->SetBinError(i+1, C2_star_mm_e[i]);}
1283 Two_star_mm->SetMarkerColor(4);
1284 Two_star_mm->SetLineColor(4);
1285 //Two_star_mm->Draw("same");
1286 //legend1->AddEntry(Two_star_mm,"--, 200 GeV","p");
1288 //Two_ex_clone_mp->Multiply(Two_ex_clone_pp);
1289 //Two_ex_clone_mp->Draw();
1290 //Two_ex_clone_mm->Multiply(Two_ex_clone_mp);
1291 //Two_ex_clone_mm->Draw();
1292 //Two_ex_clone_pp->Draw("same");
1294 //legend1->Draw("same");
1301 /////////////////////////////////////////////////////////////////////////
1302 // 3-d fitting (out,side,long). Not used for paper.
1304 TString *name1 = new TString("Explicit2_Charge1_0_Charge2_0_SC_0_M_");
1306 TString *name2 = new TString(name1->Data());
1307 TString *name3 = new TString(name1->Data());
1308 name1->Append("_ED_0_Term_1_osl_b2");// b1 (0.2<kt<0.3). b2 (0.6<kt<0.7)
1309 name2->Append("_ED_0_Term_2_osl_b2");
1310 name3->Append("_ED_0_Term_1_osl_b2_QW");
1313 TH3D *num_osl = (TH3D*)MyList->FindObject(name1->Data());
1314 TH3D *den_osl = (TH3D*)MyList->FindObject(name2->Data());
1315 den_osl->Scale(num_osl->Integral(28,40, 28,40, 28,40)/den_osl->Integral(28,40, 28,40, 28,40));
1316 TH3D *num_osl_QW = (TH3D*)MyList->FindObject(name3->Data());
1318 for(int i=0; i<BINS_OSL; i++){
1319 for(int j=0; j<BINS_OSL; j++){
1320 for(int k=0; k<BINS_OSL; k++){
1321 if(num_osl->GetBinContent(i+1,j+1,k+1) < 1 || den_osl->GetBinContent(i+1,j+1,k+1) < 1) continue;
1323 avg_q[i][j][k] = num_osl_QW->GetBinContent(i+1,j+1,k+1)/num_osl->GetBinContent(i+1,j+1,k+1);
1324 int QinvBin = int((avg_q[i][j][k])/0.005);
1325 if(QinvBin >=20) QinvBin=19;
1326 if(MomRes_term1_pp[QinvBin] ==0) continue;
1327 if(MomRes_term2_pp[QinvBin] ==0) continue;
1329 num_osl->SetBinContent(i+1,j+1,k+1, num_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term1_pp[QinvBin]);
1330 den_osl->SetBinContent(i+1,j+1,k+1, den_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term2_pp[QinvBin]);
1332 A[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1333 B[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1335 A_e[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1336 A_e[i][j][k] += pow(num_osl->GetBinContent(i+1,j+1,k+1)/MomRes_term1_pp[QinvBin] *fabs(MomRes_term1_pp[QinvBin]-1)*0.1,2);
1337 A_e[i][j][k] = sqrt(A_e[i][j][k]);
1338 B_e[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1339 B_e[i][j][k] += pow(den_osl->GetBinContent(i+1,j+1,k+1)/MomRes_term2_pp[QinvBin] *fabs(MomRes_term2_pp[QinvBin]-1)*0.1,2);
1340 B_e[i][j][k] = sqrt(B_e[i][j][k]);
1342 K_OSL[i][j][k] = CoulCorr2(+1, avg_q[i][j][k]);
1343 K_OSL_e[i][j][k] = sqrt(pow((K_OSL[i][j][k]-1)*0.02,2) + pow((K_OSL[i][j][k]-Gamov(+1,avg_q[i][j][k]))*0.01,2));
1344 //K_OSL_e[i][j][k] = 0;
1350 const int npar_osl=6;//5
1351 TMinuit MyMinuit_osl(npar_osl);
1352 MyMinuit_osl.SetFCN(fcnOSL);
1353 double OutputPar_osl[npar_osl]={0};
1354 double OutputPar_osl_e[npar_osl]={0};
1356 double par_osl[npar_osl]; // the start values
1357 double stepSize_osl[npar_osl]; // step sizes
1358 double minVal_osl[npar_osl]; // minimum bound on parameter
1359 double maxVal_osl[npar_osl]; // maximum bound on parameter
1360 string parName_osl[npar_osl];
1362 par_osl[0] = 1.0; par_osl[1]=0.5; par_osl[2]=0.; par_osl[3]=5.; par_osl[4] = 5.; par_osl[5]=5.;
1363 stepSize_osl[0] = 0.01; stepSize_osl[1] = 0.01; stepSize_osl[2] = 0.02; stepSize_osl[3] = 0.2; stepSize_osl[4] = 0.2; stepSize_osl[5] = 0.2;
1364 minVal_osl[0] = 0.8; minVal_osl[1] = 0.01; minVal_osl[2] = 0.0; minVal_osl[3] = 1.; minVal_osl[4] = 1.; minVal_osl[5] = 1.;
1365 maxVal_osl[0] = 1.2; maxVal_osl[1] = 1.0; maxVal_osl[2] = 0.99; maxVal_osl[3] = 20.; maxVal_osl[4] = 20.; maxVal_osl[5] = 20.;
1366 parName_osl[0] = "Norm"; parName_osl[1] = "Lambda"; parName_osl[2] = "G"; parName_osl[3] = "R_out"; parName_osl[4] = "R_side"; parName_osl[5] = "R_long";
1368 for (int i=0; i<npar_osl; i++){
1369 MyMinuit_osl.DefineParameter(i, parName_osl[i].c_str(), par_osl[i], stepSize_osl[i], minVal_osl[i], maxVal_osl[i]);
1371 MyMinuit_osl.FixParameter(2);// G
1372 //MyMinuit.FixParameter(1);// lambda
1373 cout<<"here!!"<<endl;
1375 // Do the minimization!
1376 cout<<"Start Three-d fit"<<endl;
1377 MyMinuit_osl.Migrad();// Minuit's best minimization algorithm
1378 cout<<"End Three-d fit"<<endl;
1379 cout<<"Chi2/NDF = "<<Chi2_OSL/(NFitPoints_OSL-MyMinuit_osl.GetNumFreePars())<<endl;
1381 for (int i=0; i<npar_osl; i++){
1382 MyMinuit_osl.GetParameter(i,OutputPar_osl[i],OutputPar_osl_e[i]);
1386 TF3 *fit3d = new TF3("fit3d",OSLfitfunction, 0,0.2, 0,0.2, 0,0.2, npar_osl);
1387 for(int i=0; i<npar_osl; i++) {fit3d->FixParameter(i,OutputPar_osl[i]);}
1388 Int_t BinsOSL = num_osl->GetNbinsX();
1389 double LimitOSL = num_osl->GetXaxis()->GetBinUpEdge(BinsOSL);
1390 TH3D *den_timesFit = new TH3D("den_timesFit","",BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL);
1391 for(int i=0; i<BinsOSL; i++){
1392 for(int j=0; j<BinsOSL; j++){
1393 for(int k=0; k<BinsOSL; k++){
1395 double qout=(i+.5)*0.005;
1396 double qside=(j+.5)*0.005;
1397 double qlong=(k+.5)*0.005;
1398 den_timesFit->SetBinContent(i+1,j+1,k+1, fit3d->Eval(qout,qside,qlong)*den_osl->GetBinContent(i+1,j+1,k+1));
1399 den_timesFit->SetBinError(i+1,j+1,k+1, 0);
1405 TH1D *num_pro=(TH1D*)num_osl->ProjectionX("num_pro",binL,binH, binL,binH);
1406 TH1D *den_pro=(TH1D*)den_osl->ProjectionX("den_pro",binL,binH, binL,binH);
1407 TH1D *dentimesFit_pro=(TH1D*)den_timesFit->ProjectionX("dentimesFit_pro",binL,binH, binL,binH);
1410 num_pro->Divide(den_pro);
1411 num_pro->SetMarkerStyle(20);
1412 num_pro->SetTitle("");
1413 num_pro->GetXaxis()->SetTitle("q_{out}");
1414 num_pro->GetYaxis()->SetTitle("C_{2}");
1415 num_pro->SetMinimum(0.97);
1416 num_pro->SetMaximum(1.48);
1417 num_pro->DrawCopy();
1419 dentimesFit_pro->Divide(den_pro);
1420 dentimesFit_pro->SetLineColor(2);
1421 dentimesFit_pro->DrawCopy("same");
1424 //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1425 //TGraph *gr2 = (TGraph*)MyMinuit.Contour(10,1,2);
1426 //gr2->SetFillColor(kYellow);
1428 //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1429 //MyMinuit.SetErrorDef(1);
1430 //TGraph *gr1 = (TGraph*)MyMinuit.Contour(10,1,2);
1431 //gr1->SetFillColor(kGreen);//38
1436 //////////////////////////////////////////////////////////////////////////////////////
1439 // To visualize the Qcut and Norm Regions
1440 //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
1441 //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);
1442 //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);
1443 //for(int bin=1; bin<=20; bin++) QcutRegion->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1444 //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1445 //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1446 //QcutRegion->SetFillColor(4);
1447 //NormRegion1->SetFillColor(2);
1448 //NormRegion2->SetFillColor(3);
1449 //Two_ex[0][0][0][0]->Draw();
1450 //QcutRegion->Draw("same");
1451 //NormRegion1->Draw("same");
1452 //NormRegion2->Draw("same");
1457 TCanvas *can2 = new TCanvas("can2", "can2",800,0,900,900);//800,0,600,900
1458 can2->SetHighLightColor(2);
1459 gStyle->SetOptFit(0111);
1460 can2->SetFillColor(10);
1461 can2->SetBorderMode(0);
1462 can2->SetBorderSize(2);
1465 can2->SetFrameFillColor(0);
1466 can2->SetFrameBorderMode(0);
1467 can2->SetFrameBorderMode(0);
1469 gPad->SetLeftMargin(0.14);
1470 gPad->SetRightMargin(0.02);
1472 TLegend *legend2 = new TLegend(.4,.67,1.,.87,NULL,"brNDC");
1473 legend2->SetBorderSize(1);
1474 legend2->SetTextSize(.06);// small .03; large .06
1475 legend2->SetFillColor(0);
1477 /////////////////////////////////////////////
1478 TH3D *C3QS_3d = new TH3D("C3QS_3d","",20,0,.98, 20,0,.1, 20,0,.1);
1479 TH3D *Combinatorics_3d = new TH3D("Combinatorics_3d","",20,0,.1, 20,0,.1, 20,0,.1);
1481 const int Q3BINS = 20;
1482 TH1D *num_withRS = new TH1D("num_withRS","",Q3BINS,0,0.2);
1483 TH1D *num_withGRS = new TH1D("num_withGRS","",Q3BINS,0,0.2);
1484 TH1D *num_withTM = new TH1D("num_withTM","",Q3BINS,0,0.2);
1485 TH1D *Cterm1 = new TH1D("Cterm1","",Q3BINS,0,0.2);
1486 TH1D *Cterm1_noMRC = new TH1D("Cterm1_noMRC","",Q3BINS,0,0.2);
1487 TH1D *Cterm2 = new TH1D("Cterm2","",Q3BINS,0,0.2);
1488 TH1D *Cterm3 = new TH1D("Cterm3","",Q3BINS,0,0.2);
1489 TH1D *Cterm4 = new TH1D("Cterm4","",Q3BINS,0,0.2);
1490 TH1D *num_QS = new TH1D("num_QS","",Q3BINS,0,0.2);
1491 TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,0.2);
1492 TH1D *Combinatorics_1d_noMRC = new TH1D("Combinatorics_1d_noMRC","",Q3BINS,0,0.2);
1493 TH1D *Coul_Riverside = new TH1D("Coul_Riverside","",Q3BINS,0,0.2);
1494 TH1D *Coul_GRiverside = new TH1D("Coul_GRiverside","",Q3BINS,0,0.2);
1495 TH1D *Coul_Omega0 = new TH1D("Coul_Omega0","",Q3BINS,0,0.2);
1496 TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,0.2);
1497 TH1D *c3_hist_STAR = new TH1D("c3_hist_STAR","",Q3BINS,0,0.2);
1498 TProfile *MomRes_2d = new TProfile("MomRes_2d","",Q3BINS,0,0.2, 0,20.,"");
1499 TProfile *MomRes_3d = new TProfile("MomRes_3d","",Q3BINS,0,0.2, 0,20.,"");
1500 TH1D *r3_num_Q3 = new TH1D("r3_num_Q3","",Q3BINS,0,0.2);
1501 TH1D *r3_den_Q3 = new TH1D("r3_den_Q3","",Q3BINS,0,0.2);
1502 r3_num_Q3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1503 r3_num_Q3->GetYaxis()->SetTitle("r_{3}");
1504 r3_num_Q3->GetYaxis()->SetTitleOffset(1.3);
1505 r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
1506 r3_num_Q3->SetMinimum(0);
1507 r3_num_Q3->SetMaximum(2.4);
1508 r3_den_Q3->SetLineColor(2);
1509 r3_num_Q3->SetMarkerStyle(20);
1510 Coul_Omega0->GetXaxis()->SetRangeUser(0,0.099);
1511 Coul_Omega0->GetXaxis()->SetLabelSize(0.04);
1512 Coul_Omega0->GetYaxis()->SetLabelSize(0.04);
1513 c3_hist_STAR->GetXaxis()->SetRangeUser(0,0.099);
1514 c3_hist_STAR->SetMinimum(0.8); c3_hist_STAR->SetMaximum(1.02);
1515 c3_hist_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1516 c3_hist_STAR->GetYaxis()->SetTitle("c_{3}*^{++-}");
1517 c3_hist_STAR->GetYaxis()->SetTitleOffset(1.6);
1518 if(SameCharge) {Cterm1_noMRC->Sumw2(); Combinatorics_1d_noMRC->Sumw2();}
1520 double num_QS_e[Q3BINS]={0};
1521 double c3_e[Q3BINS]={0};
1522 double r3_e[Q3BINS]={0};
1524 // CB=Charge Bin. 0 means pi-, 1 means pi+
1525 int CB1=0, CB2=0, CB3=0;
1526 int CP12=1, CP13=1, CP23=1;
1528 if(SameCharge) {CB1=0; CB2=0; CB3=0;}
1529 else {CB1=0; CB2=0; CB3=1; CP12=+1, CP13=-1, CP23=-1;}
1531 if(SameCharge) {CB1=1; CB2=1; CB3=1;}
1532 else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;}
1535 // SC = species combination. Always 0 meaning pi-pi-pi.
1539 ReadCoulCorrections(SourceType, RValue, bValue, 10);// switch to full kt range, 10.
1540 //ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation (STAR method testing)
1541 TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",20,0,0.2);
1542 TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",20,0,0.2);
1544 double value_num[2]={0};
1545 double value_num_e[2]={0};
1546 double value_den[2]={0};
1547 double N3_QS[2]={0};
1548 double N3_QS_e[2]={0};
1549 double OutTriplets=0, InTriplets=0;
1550 for(int i=2; i<=20; i++){// bin number
1551 //double Q12_m = (i-0.5)*(0.005);// geometric center
1552 double Q12 = AvgQinvSS[i-1]; if(!SameCharge && CHARGE==+1) {Q12 = AvgQinvOS[i-1];}// true center
1553 int Q12bin = int(Q12/0.005)+1;
1555 for(int j=2; j<=20; j++){// bin number
1556 //double Q13_m = (j-0.5)*(0.005);// geometric center
1557 double Q13 = AvgQinvSS[j-1]; if(!SameCharge) {Q13 = AvgQinvOS[j-1];}// true center
1558 int Q13bin = int(Q13/0.005)+1;
1560 for(int k=2; k<=20; k++){// bin number
1561 //double Q23_m = (k-0.5)*(0.005);// geometric center
1562 double Q23 = AvgQinvSS[k-1]; if(!SameCharge && CHARGE==-1) {Q23 = AvgQinvOS[k-1];}// true center
1563 int Q23bin = int(Q23/0.005)+1;
1565 //if(fabs(i-j)>3 || fabs(i-k)>3 || fabs(j-k)>3) continue;// testing
1567 double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
1568 int Q3bin = Cterm1->GetXaxis()->FindBin(Q3);
1570 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
1571 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
1575 double G_12 = Gamov(CP12, Q12);// point-source Coulomb correlation
1576 double G_13 = Gamov(CP13, Q13);
1577 double G_23 = Gamov(CP23, Q23);
1578 double K2_12 = CoulCorr2(CP12, Q12);// finite-source Coulomb+Strong correlation from Therminator.
1579 double K2_13 = CoulCorr2(CP13, Q13);
1580 double K2_23 = CoulCorr2(CP23, Q23);
1581 double K3 = 1.;// 3-body Coulomb+Strong correlation
1582 if(GRS) {// GRS approach
1583 K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
1584 if(CHARGE==+1 && !SameCharge) K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
1586 K3 = CoulCorrOmega0(SameCharge, Q12, Q13, Q23);
1587 if(CHARGE==+1 && !SameCharge) K3 = CoulCorrOmega0(SameCharge, Q23, Q12, Q13);
1589 if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;}
1592 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
1593 double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
1594 double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
1595 double TERM4=Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event
1596 double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
1597 if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
1598 if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
1600 if(Q3>0.08 && Q3<0.09){// just for testing
1601 if(Q12>0.08 || Q13>0.08 || Q23>0.08) OutTriplets++;
1605 // apply momentum resolution correction
1606 if(!MCcase && !GeneratedSignal){
1608 // 3d momentum resolution corrections
1609 TERM1 *= (MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1610 TERM2 *= (MomRes3d[0][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1611 TERM3 *= (MomRes3d[0][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1612 TERM4 *= (MomRes3d[0][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1613 TERM5 *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1614 // Triple produce of 1-d momentum resolution corrections (less accurate).
1615 /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1];
1616 TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];
1617 TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term2_pp[k-1];
1618 TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term1_pp[k-1];
1619 TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];*/
1620 //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1]);
1621 //MomRes_3d->Fill(Q3, MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
1623 if(CHARGE==-1){// pi-pi-pi+
1624 TERM1 *= (MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1625 TERM2 *= (MomRes3d[1][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1626 TERM3 *= (MomRes3d[1][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1627 TERM4 *= (MomRes3d[1][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1628 TERM5 *= (MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1629 /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1];
1630 TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];
1631 TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_mp[k-1];
1632 TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_mp[k-1];
1633 TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];*/
1635 TERM1 *= (MomRes3d[1][0]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1636 TERM2 *= (MomRes3d[1][3]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1637 TERM3 *= (MomRes3d[1][2]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1638 TERM4 *= (MomRes3d[1][1]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1639 TERM5 *= (MomRes3d[1][4]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1640 /*TERM1 *= MomRes_term1_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_pp[k-1];
1641 TERM2 *= MomRes_term1_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];
1642 TERM3 *= MomRes_term2_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_pp[k-1];
1643 TERM4 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_pp[k-1];
1644 TERM5 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];*/
1646 //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1]);// always treats 12 as ss pair
1647 //MomRes_3d->Fill(Q3, MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
1652 for(int LamType=0; LamType<2; LamType++){
1654 if(LamType==1){TwoFrac=TwoFracr3; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);}// r3 assignment
1656 TwoFrac = OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);// Assign to C2 global fit values found previously
1659 // Purify. Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
1660 N3_QS[LamType] = TERM1;
1661 N3_QS[LamType] -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1662 N3_QS[LamType] -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
1663 N3_QS[LamType] /= ThreeFrac;
1664 N3_QS[LamType] /= K3;
1666 if(LamType==0) num_QS->Fill(Q3, N3_QS[LamType]);
1668 // Isolate 3-pion cumulant
1669 value_num[LamType] = N3_QS[LamType];
1670 value_num[LamType] -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac;
1671 value_num[LamType] -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac;
1672 value_num[LamType] -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac;
1673 value_num[LamType] += 2*TERM5;
1677 if(LamType==1 && SameCharge) {
1678 value_den[LamType] = PionNorm[CB1]->GetBinContent(i,j,k);// Raw r3 denominator
1680 // Momentum Resolution Correction Systematic variation. Only important when MRCShift != 1.0.
1681 double denMRC1 = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1682 denMRC1 *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1683 denMRC1 *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1684 double denMRC2 = (C2ssRaw->GetBinContent(i)*((MomRes_C2_pp[i-1]-1)*MRCShift+1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1685 denMRC2 *= (C2ssRaw->GetBinContent(j)*((MomRes_C2_pp[j-1]-1)*MRCShift+1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1686 denMRC2 *= (C2ssRaw->GetBinContent(k)*((MomRes_C2_pp[k-1]-1)*MRCShift+1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1687 // Non-femto background correction (Mini-Jet).
1688 double denMJ = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - (MJ_bkg_ss->Eval(Q12)-1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1689 denMJ *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - (MJ_bkg_ss->Eval(Q13)-1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1690 denMJ *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - (MJ_bkg_ss->Eval(Q23)-1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1692 double den_ratio = sqrt(fabs(denMRC2))/sqrt(fabs(denMRC1));
1693 value_den[LamType] *= den_ratio;// apply Momentum Resolution correction systematic variation if any.
1694 double den_ratioMJ = sqrt(fabs(denMJ))/sqrt(fabs(denMRC1));
1695 if(IncludeMJcorrection) value_den[LamType] *= den_ratioMJ;// Non-femto bkg correction
1697 //value_den[LamType] -= TPNRejects->GetBinContent(i,j,k);// Testing for 0-5% only to estimate the effect when C2^QS < 1.0
1698 value_den[LamType] *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;// additional momentum resolution correction
1703 N3_QS_e[LamType] = TERM1;
1704 N3_QS_e[LamType] += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2);
1705 N3_QS_e[LamType] += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2));
1706 N3_QS_e[LamType] /= pow(ThreeFrac,2);
1707 N3_QS_e[LamType] /= pow(K3,2);
1708 if(LamType==0) num_QS_e[Q3bin-1] += N3_QS_e[LamType];
1710 value_num_e[LamType] = N3_QS_e[LamType];
1711 value_num_e[LamType] += (pow(1/K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(TERM5),2));
1712 value_num_e[LamType] += (pow(1/K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(TERM5),2));
1713 value_num_e[LamType] += (pow(1/K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(TERM5),2));
1714 value_num_e[LamType] += pow(2*sqrt(TERM5),2);
1715 if(LamType==0) c3_e[Q3bin-1] += value_num_e[LamType] + TERM5;// add baseline stat error
1716 else r3_e[Q3bin-1] += value_num_e[LamType];
1720 c3_hist->Fill(Q3, value_num[0] + TERM5);// for cumulant correlation function
1721 C3QS_3d->SetBinContent(i,j,k, N3_QS[0]);
1722 C3QS_3d->SetBinError(i,j,k, N3_QS_e[0]);
1724 Coul_Riverside->Fill(Q3, G_12*G_13*G_23*TERM5);
1725 Coul_GRiverside->Fill(Q3, TERM5*CoulCorrGRS(SameCharge, Q12, Q13, Q23));
1726 Coul_Omega0->Fill(Q3, K3*TERM5);
1727 num_withGRS->Fill(Q3, N3_QS[0]);
1728 Cterm1->Fill(Q3, TERM1);
1729 Cterm2->Fill(Q3, TERM2);
1730 Cterm3->Fill(Q3, TERM3);
1731 Cterm4->Fill(Q3, TERM4);
1732 Combinatorics_1d->Fill(Q3, TERM5);
1733 Combinatorics_3d->Fill(Q12,Q13,Q23, TERM5);
1734 r3_num_Q3->Fill(Q3, value_num[1]);
1735 r3_den_Q3->Fill(Q3, value_den[1]);
1736 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));
1737 Cterm1_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k));
1738 Combinatorics_1d_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1739 double cumulant_STAR = Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k)/(K2_12*K2_13*K2_23);
1740 cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k)/(K2_12);
1741 cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k)/(K2_13);
1742 cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k)/(K2_23);
1743 c3_hist_STAR->Fill(Q3_m, cumulant_STAR + 3*Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1745 ///////////////////////////////////////////////////////////
1746 // Edgeworth 3-pion expection. Not important for r3.
1747 //double radius_exp = (11-(MomResCentBin-1))/FmToGeV;
1748 //TwoFrac=0.68; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1749 double radius_exp = (OutputPar[3])/FmToGeV;
1750 TwoFrac=OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1751 double arg12 = Q12*radius_exp;
1752 double arg13 = Q13*radius_exp;
1753 double arg23 = Q23*radius_exp;
1754 /*Float_t EW12 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1755 EW12 += 0.16/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1756 Float_t EW13 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1757 EW13 += 0.16/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1758 Float_t EW23 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1759 EW23 += 0.16/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);*/
1761 Float_t EW12 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1762 EW12 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1763 Float_t EW13 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1764 EW13 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1765 Float_t EW23 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1766 EW23 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
1768 double cumulant_exp=0, term1_exp=0;
1770 cumulant_exp = (exp(-pow(radius_exp*Q12,2))*pow(EW12,2)+exp(-pow(radius_exp*Q13,2))*pow(EW13,2)+exp(-pow(radius_exp*Q23,2))*pow(EW23,2))*TERM5;
1771 cumulant_exp += 2*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5;
1772 term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1773 term1_exp += TwoFrac*(1-OneFrac)*(K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2))+K2_13*(1+exp(-pow(radius_exp*Q13,2))*pow(EW13,2))+K2_23*(1+exp(-pow(radius_exp*Q23,2))*pow(EW23,2)))*TERM5;
1774 term1_exp += ThreeFrac*K3*cumulant_exp;
1775 //term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1777 cumulant_exp = (exp(-pow(radius_exp*Q12,2))*pow(EW12,2))*TERM5 + TERM5;
1778 term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1779 term1_exp += TwoFrac*(1-OneFrac)*(K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)) + K2_13 + K2_23)*TERM5;
1780 term1_exp += ThreeFrac*K3*cumulant_exp;
1781 term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1782 //term1_exp = ((1-TwoFrac) + TwoFrac*K2_13)*TERM5;
1785 GenSignalExpected_num->Fill(Q3, term1_exp);
1786 GenSignalExpected_den->Fill(Q3, TERM5);
1787 ///////////////////////////////////////////////////////////
1793 cout<<"OutTriplets: "<<OutTriplets<<" InTriplets: "<<InTriplets<<endl;
1794 ////////////////////////////
1796 // Intermediate steps
1797 num_withRS->Sumw2();
1798 num_withGRS->Sumw2();
1799 num_withTM->Sumw2();
1804 Combinatorics_1d->Sumw2();
1805 Combinatorics_3d->Sumw2();
1808 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]));}
1809 for(int i=0; i<Q3BINS; i++) {c3_hist_STAR->SetBinError(i+1, sqrt(c3_e[i]));}// approximate error
1810 num_withRS->Divide(Combinatorics_1d);
1811 num_withGRS->Divide(Combinatorics_1d);
1812 num_withTM->Divide(Combinatorics_1d);
1813 for(int q3bin=1; q3bin<=r3_num_Q3->GetNbinsX(); q3bin++){
1814 r3_num_Q3->SetBinError(q3bin, sqrt(r3_e[q3bin-1]));
1816 Cterm1->Divide(Combinatorics_1d);
1817 Cterm1_noMRC->Divide(Combinatorics_1d_noMRC);
1818 Cterm2->Divide(Combinatorics_1d);
1819 Cterm3->Divide(Combinatorics_1d);
1820 Cterm4->Divide(Combinatorics_1d);
1821 c3_hist->Divide(Combinatorics_1d);
1822 c3_hist_STAR->Divide(Combinatorics_1d_noMRC);
1823 num_QS->Divide(Combinatorics_1d);
1824 r3_num_Q3->Divide(r3_den_Q3);
1825 GenSignalExpected_num->Sumw2();
1826 GenSignalExpected_num->Divide(GenSignalExpected_den);
1830 Coul_Riverside->Divide(Combinatorics_1d);
1831 Coul_GRiverside->Divide(Combinatorics_1d);
1832 Coul_Omega0->Divide(Combinatorics_1d);
1833 for(int ii=1; ii<=Coul_Omega0->GetNbinsX(); ii++){
1834 Coul_Riverside->SetBinError(ii,0.000001);
1835 Coul_GRiverside->SetBinError(ii,0.000001);
1836 Coul_Omega0->SetBinError(ii,0.000001);
1838 ////////////////////////////
1842 ///////////////////////////////////////////////////////////////////////////////////////////////////
1844 Coul_Riverside->SetLineColor(1);
1845 Coul_GRiverside->SetLineColor(2);
1846 Coul_Omega0->SetLineColor(4);
1848 //Coul_Riverside->Draw();
1849 //legend1->AddEntry(Coul_Riverside,"Riverside","l");
1850 //Coul_GRiverside->Draw("same");
1851 //legend1->AddEntry(Coul_GRiverside,"Generalized Riverside","l");
1852 //Coul_Omega0->Draw("same");
1853 //legend1->AddEntry(Coul_Omega0,"Omega0","l");
1856 Cterm1->SetMarkerStyle(20);
1857 Cterm1->SetMarkerColor(4);
1858 Cterm1->SetLineColor(4);
1859 Cterm1->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1860 Cterm1->GetYaxis()->SetTitle("C_{3}");
1861 //Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}");
1862 Cterm1->SetMinimum(0.97);
1863 Cterm1->SetMaximum(1.7);// 6.1 for same-charge
1864 Cterm1->GetXaxis()->SetRangeUser(0,.095);
1865 Cterm1->GetYaxis()->SetTitleOffset(1.4);
1867 legend2->AddEntry(Cterm1,"C_{3} raw","p");
1869 Cterm2->SetMarkerStyle(20);
1870 Cterm2->SetMarkerColor(7);
1872 Cterm3->SetMarkerStyle(25);
1873 Cterm3->SetMarkerColor(8);
1874 Cterm4->SetMarkerStyle(26);
1875 Cterm4->SetMarkerColor(9);
1876 //Cterm2->Draw("same");
1877 //Cterm3->Draw("same");
1878 //Cterm4->Draw("same");
1879 //legend2->AddEntry(Cterm1,"++-","p");
1883 double C3points_HIJING_mmm[10]={0, 0.961834, 1.01827, 0.999387, 1.00202, 1.00081, 1.00082, 1.00037, 0.999074, 0.999099};
1884 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};
1885 double C3points_HIJING_ppp[10]={0, 1.13015, 1.00623, 1.00536, 1.00293, 0.999964, 1.00116, 1.0007, 1.00037, 1.00105};
1886 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};
1887 TH1D *C3_HIJING_mmm=(TH1D*)Cterm1->Clone();
1888 TH1D *C3_HIJING_ppp=(TH1D*)Cterm1->Clone();
1889 for(int i=0; i<10; i++){
1890 C3_HIJING_mmm->SetBinContent(i+1, C3points_HIJING_mmm[i]);
1891 C3_HIJING_mmm->SetBinError(i+1, C3points_HIJING_mmm_e[i]);
1892 C3_HIJING_ppp->SetBinContent(i+1, C3points_HIJING_ppp[i]);
1893 C3_HIJING_ppp->SetBinError(i+1, C3points_HIJING_ppp_e[i]);
1895 C3_HIJING_mmm->SetMarkerColor(2);
1896 C3_HIJING_mmm->SetLineColor(2);
1897 C3_HIJING_ppp->SetMarkerColor(4);
1898 C3_HIJING_ppp->SetLineColor(4);
1899 //C3_HIJING_mmm->Draw("same");
1900 //C3_HIJING_ppp->Draw("same");
1901 //legend2->AddEntry(C3_HIJING_mmm,"---","p");
1902 //legend2->AddEntry(C3_HIJING_ppp,"+++","p");
1905 num_QS->SetMarkerStyle(24);
1906 num_QS->SetMarkerColor(4);
1907 num_QS->SetLineColor(4);
1908 num_QS->GetXaxis()->SetTitle("Q_{3}");
1909 num_QS->GetYaxis()->SetTitle("C_{3}^{QS}");
1910 num_QS->GetXaxis()->SetRangeUser(0,.12);
1911 num_QS->SetMaximum(6);
1912 //num_QS->Draw("same");
1913 //legend2->AddEntry(num_QS,"C_{3}^{QS}","p");
1915 c3_hist->GetXaxis()->SetTitle("Q_{3}");
1916 c3_hist->GetYaxis()->SetTitle("c_{3}^{++-}");
1917 c3_hist->GetXaxis()->SetRangeUser(0,.15);
1918 c3_hist->SetMarkerStyle(25);
1919 c3_hist->SetMarkerColor(2);
1920 c3_hist->SetLineColor(2);
1921 c3_hist->SetMaximum(3);
1922 c3_hist->SetMinimum(.9);
1923 if(!MCcase) c3_hist->Draw("same");
1924 //legend2->AddEntry(c3_hist,"#font[12]{c}_{3} (cumulant correlation)","p");
1927 GenSignalExpected_num->SetMarkerStyle(20);
1928 //GenSignalExpected_num->Draw("same");
1929 //legend2->AddEntry(GenSignalExpected_num,"#kappa_{3}=0.24, #kappa_{4}=0.16, #lambda=0.68, R=6 fm","p");
1930 //legend2->AddEntry(GenSignalExpected_num,"Edeworth expectation (fully chaotic)","p");
1932 //MomRes_2d->SetMarkerStyle(20);
1933 //MomRes_3d->SetMarkerStyle(20);
1934 //MomRes_3d->SetMarkerColor(4);
1935 //MomRes_2d->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1936 //MomRes_2d->GetYaxis()->SetTitle("< MRC >");
1937 //MomRes_2d->SetTitle("");
1938 //MomRes_2d->Draw();
1939 //legend2->AddEntry(MomRes_2d,"2D: Triple pair product","p");
1940 //MomRes_3d->Draw("same");
1941 //legend2->AddEntry(MomRes_3d,"3D","p");
1943 //legend2->Draw("same");
1947 ////////// C3 systematics
1948 // C3 --- base, M0, (0.035 TTC for all)
1949 //double C3_base[10]={0, 1.63072, 1.6096, 1.43731, 1.28205, 1.17777, 1.11973, 1.07932, 1.05459, 1.04029};
1950 //double C3_base_e[10]={0, 0.022528, 0.00387504, 0.00115667, 0.000423799, 0.000238973, 0.000143993, 9.71502e-05, 7.02007e-05, 5.31267e-05};
1951 // C3 --- base, M0, Pos B (0.035 TTC for all)
1952 //double C3_base[10]={0, 1.62564, 1.60461, 1.438, 1.28234, 1.17827, 1.12009, 1.07953, 1.05474, 1.04037};
1953 //double C3_base_e[10]={0, 0.0239233, 0.00409002, 0.0012215, 0.000446701, 0.000251769, 0.000151651, 0.000102284, 7.38993e-05, 5.59212e-05};
1954 // C3 --+ base, M0, (0.035 TTC for all)
1955 //double C3_base[10]={0, 1.66016, 1.41961, 1.25229, 1.16459, 1.11254, 1.08012, 1.05768, 1.04265, 1.0332};
1956 //double C3_base_e[10]={0, 0.00539779, 0.00111398, 0.000387926, 0.000192906, 0.00011428, 7.24765e-05, 5.09126e-05, 3.76421e-05, 2.87533e-05};
1958 // C3 --- base, M0, (New Standard: 0.04 TTC )
1959 //double C3_base[10]={0, 1.63903, 1.60254, 1.43381, 1.2794, 1.17603, 1.11806, 1.07806, 1.05345, 1.03936};
1960 //double C3_base_e[10]={0, 0.0322796, 0.00438361, 0.00122249, 0.000424557, 0.000233965, 0.000139058, 9.28136e-05, 6.66159e-05, 5.01816e-05};
1961 // C3 --- base, M1, (New Standard: 0.04 TTC )
1962 //double C3_base[10]={0, 1.6127, 1.65561, 1.49508, 1.33093, 1.21458, 1.14708, 1.099, 1.06864, 1.05064};
1963 //double C3_base_e[10]={0, 0.0425447, 0.0061176, 0.00172948, 0.000600699, 0.000329342, 0.000194832, 0.000129427, 9.25599e-05, 6.95395e-05};
1964 // C3 --- base, M2, (New Standard: 0.04 TTC )
1965 //double C3_base[10]={0, 1.6509, 1.71863, 1.54652, 1.38092, 1.25226, 1.17549, 1.12068, 1.08408, 1.06236};
1966 //double C3_base_e[10]={0, 0.0981473, 0.0143699, 0.00404612, 0.00141189, 0.000770764, 0.000453944, 0.000300452, 0.000214068, 0.000160448};
1967 // C3 --- base, M9, (New Standard: 0.04 TTC )
1968 //double C3_base[10]={0, 2.41982, 2.18303, 1.93205, 1.80399, 1.60955, 1.49305, 1.38565, 1.29873, 1.23626};
1969 //double C3_base_e[10]={0, 1.60227, 0.177274, 0.0501455, 0.018456, 0.00998147, 0.00583719, 0.00379465, 0.00264116, 0.0019391};
1971 // C3 --+ base, M0, (New Standard: 0.04 TTC )
1972 //double C3_base[10]={0, 1.66087, 1.41943, 1.25081, 1.16313, 1.11143, 1.07917, 1.05701, 1.04215, 1.0328};
1973 //double C3_base_e[10]={0, 0.00584743, 0.00111278, 0.000374009, 0.00018315, 0.000107523, 6.78669e-05, 4.75116e-05, 3.50489e-05, 2.67328e-05};
1975 // HIJING C3 --- base, M0
1976 //double C3_base[10]={0, 0.970005, 1.00655, 1.00352, 1.00291, 1.00071, 1.0002, 0.999524, 0.999404, 0.999397};
1977 //double C3_base_e[10]={0, 0.050845, 0.0099602, 0.00334862, 0.00138008, 0.000841743, 0.000531776, 0.000371712, 0.000274716, 0.00021};
1978 // HIJING C3 --- base, M1
1979 //double C3_base[10]={0, 1.03657, 1.00199, 0.997984, 1.00067, 1.0006, 0.999901, 0.999967, 0.999792, 0.999777};
1980 //double C3_base_e[10]={0, 0.0634232, 0.0117204, 0.0039446, 0.00163131, 0.000996638, 0.000629662, 0.000440266, 0.00032534, 0.000249};
1981 // HIJING C3 --- base, M2
1982 //double C3_base[10]={0, 1.34345, 1.04226, 1.0278, 0.99582, 1.00554, 1.00296, 1.00057, 1.00271, 1.00152};
1983 //double C3_base_e[10]={0, 0.363559, 0.0503531, 0.0170117, 0.00679185, 0.00419035, 0.00264603, 0.00184587, 0.00136663, 0.00104772};
1984 // HIJING C3 --- base, M3
1985 double C3_base[10]={0, 0.890897, 1.05222, 1.00461, 1.01364, 0.998981, 1.00225, 1.00305, 1.00235, 1.00043};
1986 double C3_base_e[10]={0, 0.297124, 0.0604397, 0.0195066, 0.00812906, 0.00490835, 0.00310751, 0.00217408, 0.00160575, 0.00123065};
1987 TH1D *C3baseHisto=(TH1D*)Cterm1->Clone();
1988 for(int ii=0; ii<10; ii++){
1989 C3baseHisto->SetBinContent(ii+1, C3_base[ii]);
1990 C3baseHisto->SetBinError(ii+1, C3_base_e[ii]);
1993 cout<<"C3 values:"<<endl;
1994 for(int ii=0; ii<10; ii++){
1995 cout<<Cterm1->GetBinContent(ii+1)<<", ";
1998 cout<<"C3 errors:"<<endl;
1999 for(int ii=0; ii<10; ii++){
2000 cout<<Cterm1->GetBinError(ii+1)<<", ";
2003 TH1D *C3_Sys=(TH1D*)Cterm1->Clone();
2004 for(int ii=1; ii<10; ii++){
2005 if(C3_base[ii] ==0) {
2006 C3_Sys->SetBinContent(ii+1, 100.);
2007 C3_Sys->SetBinError(ii+1, 100.);
2010 C3_Sys->SetBinContent(ii+1, (C3_base[ii]-C3_Sys->GetBinContent(ii+1))/C3_base[ii]);
2011 C3_Sys->SetBinError(ii+1, sqrt(fabs(pow(C3_Sys->GetBinError(ii+1),2) - pow(C3_base_e[ii],2))));
2013 gStyle->SetOptFit(111);
2014 C3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2015 C3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2016 C3_Sys->GetYaxis()->SetTitleOffset(2);
2017 C3_Sys->SetMinimum(-0.01);
2018 C3_Sys->SetMaximum(0.01);
2020 TF1 *C3lineSys=new TF1("C3lineSys","pol3",0,0.099);
2021 //C3_Sys->Fit(C3lineSys,"MEQ","",0,0.099);
2024 // Plotting +++ added with --- for HIJING
2025 TLegend *legendC3Hijing = new TLegend(.5,.15,.9,.3,NULL,"brNDC");
2026 legendC3Hijing->SetBorderSize(1);
2027 legendC3Hijing->SetTextSize(.03);// small .03; large .06
2028 legendC3Hijing->SetFillColor(0);
2030 Cterm1->Add(C3baseHisto);
2032 Cterm1->SetMinimum(0.9);
2033 Cterm1->SetMaximum(1.1);
2035 legendC3Hijing->AddEntry(Cterm1,"same-charge, HIJING","p");
2036 legendC3Hijing->Draw("same");
2040 ////////// c3 systematics
2041 // c3 --- base, M0, (New Standard: 0.04 TTC )
2042 double c3_base[10]={0, 1.86014, 1.47533, 1.23733, 1.09944, 1.04145, 1.01693, 1.00715, 1.00253, 1.00111};
2043 double c3_base_e[10]={0, 0.104645, 0.0120917, 0.00333303, 0.00118126, 0.0006692, 0.000405246, 0.000274163, 0.000198507, 0.000150258};
2045 cout<<"c3 values:"<<endl;
2046 for(int ii=0; ii<10; ii++){
2047 cout<<c3_hist->GetBinContent(ii+1)<<", ";
2050 cout<<"c3 errors:"<<endl;
2051 for(int ii=0; ii<10; ii++){
2052 cout<<c3_hist->GetBinError(ii+1)<<", ";
2055 TH1D *c3_Sys=(TH1D*)c3_hist->Clone();
2056 for(int ii=1; ii<10; ii++){
2057 if(c3_base[ii] ==0) {
2058 c3_Sys->SetBinContent(ii+1, 100.);
2059 c3_Sys->SetBinError(ii+1, 100.);
2062 c3_Sys->SetBinContent(ii+1, (c3_base[ii]-c3_Sys->GetBinContent(ii+1))/c3_base[ii]);
2063 c3_Sys->SetBinError(ii+1, sqrt(fabs(pow(c3_Sys->GetBinError(ii+1),2) - pow(c3_base_e[ii],2))));
2065 gStyle->SetOptFit(111);
2066 c3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2067 c3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2068 c3_Sys->GetYaxis()->SetTitleOffset(2);
2069 c3_Sys->SetMinimum(-0.01);
2070 c3_Sys->SetMaximum(0.01);
2072 TF1 *c3lineSys=new TF1("c3lineSys","pol3",0,0.099);
2073 c3_Sys->Fit(c3lineSys,"MEQ","",0,0.099);
2078 TPad *pad1 = new TPad("pad1","pad1",0.0,0.0,1.,1.);
2083 pad1->SetTopMargin(0.02);//0.05
2084 //pad1->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2085 pad1->SetRightMargin(0.01);//1e-2
2086 pad1->SetBottomMargin(0.06);//0.12
2087 pad1->Divide(1,2,0,0);
2090 gPad->SetLeftMargin(0.14);
2091 gPad->SetRightMargin(0.02);
2093 Coul_Omega0->SetMinimum(0.48);
2094 Coul_Omega0->SetMaximum(1.01);
2096 Coul_Omega0->SetMinimum(0.99);
2097 Coul_Omega0->SetMaximum(1.32);
2099 Coul_Omega0->SetMarkerStyle(20);
2100 Coul_Omega0->SetMarkerColor(4);
2101 Coul_GRiverside->SetMarkerStyle(25);
2102 Coul_GRiverside->SetMarkerColor(2);
2103 Coul_Riverside->SetMarkerStyle(23);
2104 Coul_Omega0->Draw();
2105 //Coul_Riverside->Draw("same");
2106 Coul_GRiverside->Draw("same");
2107 legend2->AddEntry(Coul_Omega0,"#Omega_{0}","p");
2108 legend2->AddEntry(Coul_GRiverside,"Generalized Riverside","p");
2109 //legend2->AddEntry(Coul_Riverside,"Riverside","p");
2110 legend2->Draw("same");
2111 TLatex *K3Label = new TLatex(-0.011,0.92,"K_{3}");// -0.011,0.92 (ss), 1.26 (os)
2112 K3Label->SetTextSize(0.08);
2113 K3Label->SetTextAngle(90);
2117 gPad->SetLeftMargin(0.14);
2118 gPad->SetRightMargin(0.02);
2119 gPad->SetBottomMargin(0.13);
2120 TH1D *K3_compOmega0 = (TH1D*)Coul_Omega0->Clone();
2121 TH1D *K3_compGRS = (TH1D*)Coul_GRiverside->Clone();
2122 for(int ii=0; ii<K3_compOmega0->GetNbinsX(); ii++){
2123 K3_compOmega0->SetBinContent(ii+1, K3_compOmega0->GetBinContent(ii+1)-1.0);
2124 K3_compGRS->SetBinContent(ii+1, K3_compGRS->GetBinContent(ii+1)-1.0);
2126 K3_compOmega0->Add(K3_compGRS,-1);
2127 K3_compOmega0->Divide(K3_compGRS);
2128 K3_compOmega0->SetMarkerStyle(20);
2129 K3_compOmega0->SetMarkerColor(1);
2130 K3_compOmega0->SetLineColor(1);
2131 K3_compOmega0->SetMinimum(-0.12);// -0.021
2132 K3_compOmega0->SetMaximum(0.12);// 0.021
2133 K3_compOmega0->SetBinContent(1,-100);
2134 K3_compOmega0->Draw();
2135 TLatex *RatioLabel = new TLatex(-.011,-0.05,"#Delta (K_{3}-1)/(K_{3}(#Omega_{0})-1)");// -0.011,0.04
2136 RatioLabel->SetTextSize(0.08);
2137 RatioLabel->SetTextAngle(90);
2139 TLatex *Q3Label = new TLatex(.065,-0.147,"Q_{3} (GeV/c)");//.065,-0.147
2140 Q3Label->SetTextSize(0.08);
2145 /////////////////////////////
2148 for(int i=0; i<BINRANGE_C2global; i++){
2149 for(int j=0; j<BINRANGE_C2global; j++){
2150 for(int k=0; k<BINRANGE_C2global; k++){
2151 C3ssFitting[i][j][k] = C3QS_3d->GetBinContent(i+1,j+1,k+1);
2152 C3ssFitting_e[i][j][k] = C3QS_3d->GetBinError(i+1,j+1,k+1);
2157 TF3 *C3ssFitResult=new TF3("C3ssFitResult",C3ssFitFunction,0,0.2,0,0.2,0,0.2,npar);
2158 for(int i=0; i<npar; i++) {
2159 if(i<=5) C3ssFitResult->FixParameter(i, OutputPar[i]);
2160 else C3ssFitResult->FixParameter(i, OutputPar[i]-OutputPar_e[i]);
2164 const int NbinsC3_1d=Q3BINS;
2165 TH1D *C3fit_1d=new TH1D("C3fit_1d","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2166 TH1D *C3fit_1d_den=new TH1D("C3fit_1d_den","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2167 double C3_err[NbinsC3_1d]={0};
2170 for(int i=0; i<C3QS_3d->GetNbinsX(); i++){// bin number
2171 for(int j=0; j<C3QS_3d->GetNbinsY(); j++){// bin number
2172 for(int k=0; k<C3QS_3d->GetNbinsZ(); k++){// bin number
2174 if(i==0 || j==0 || k==0) continue;
2175 //if(i==1 || j==1 || k==1) continue;
2176 double Q3 = sqrt(pow(AvgQinvSS[i],2) + pow(AvgQinvSS[j],2) + pow(AvgQinvSS[k],2));
2177 C3fit_1d->Fill(Q3, C3ssFitResult->Eval(AvgQinvSS[i],AvgQinvSS[j],AvgQinvSS[k])*Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2178 C3fit_1d_den->Fill(Q3, Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2183 C3fit_1d->Divide(C3fit_1d_den);
2184 for(int q3bin=1; q3bin<=num_QS->GetNbinsX(); q3bin++) C3fit_1d->SetBinError(q3bin, 0.001);// non-zero to display marker
2185 C3fit_1d->SetMarkerStyle(22);
2186 C3fit_1d->SetLineColor(4);
2187 //C3fit_1d->Draw("same");
2189 TF1 *lineC3 = new TF1("lineC3","6",0,10);
2190 lineC3->SetLineColor(4);
2191 lineC3->SetLineStyle(2);
2192 //lineC3->Draw("same");
2193 //legend2->AddEntry(lineC3,"Chaotic Limit C_{3}^{QS}","l");
2195 TF1 *linec3 = new TF1("linec3","3",0,10);
2196 linec3->SetLineColor(2);
2197 linec3->SetLineStyle(2);
2198 //linec3->Draw("same");
2199 //legend2->AddEntry(linec3,"Chaotic Limit #font[12]{c}_{3}","l");
2202 //legend2->Draw("same");
2205 /////////////////////////////////////////////////////////////////
2206 /////////////////////////////////////////////////////////////////
2208 TCanvas *can3 = new TCanvas("can3", "can3",1600,0,700,700);
2209 can3->SetHighLightColor(2);
2210 can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
2211 //gStyle->SetOptFit(0111);
2212 can3->SetFillColor(10);
2213 can3->SetBorderMode(0);
2214 can3->SetBorderSize(2);
2217 can3->SetFrameFillColor(0);
2218 can3->SetFrameBorderMode(0);
2219 can3->SetFrameBorderMode(0);
2221 TPad *pad3 = new TPad("pad3","pad3",0.0,0.0,1.,1.);
2226 pad3->SetTopMargin(0.02);//0.05
2227 //pad3->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2228 pad3->SetRightMargin(0.01);//1e-2
2229 pad3->SetBottomMargin(0.06);//0.12
2232 gPad->SetLeftMargin(0.14);
2233 gPad->SetRightMargin(0.02);
2235 TF1 *QuarticFit = new TF1("QuarticFit","[0]*(1-[1]*pow(x,4))",0,.1);
2236 QuarticFit->SetParameter(0,2); QuarticFit->SetParameter(1,0);
2237 QuarticFit->SetLineColor(4);
2238 QuarticFit->SetParName(0,"I(Q^{4})"); QuarticFit->SetParName(1,"a(Q^{4})");
2239 TF1 *QuadraticFit = new TF1("QuadraticFit","[0]*(1-[1]*pow(x,2))",0,.1);
2240 QuadraticFit->SetParameter(0,2); QuadraticFit->SetParameter(1,0);
2241 QuadraticFit->SetParName(0,"I(Q^{2})"); QuadraticFit->SetParName(1,"a(Q^{2})");
2246 r3_num_Q3->SetMinimum(0); r3_num_Q3->SetMaximum(2.5);
2247 r3_num_Q3->GetXaxis()->SetRangeUser(0.01,0.09);
2250 r3_num_Q3->Fit(QuarticFit,"IME","",0,0.1);
2252 TPaveStats *Quartic_stats=(TPaveStats*)r3_num_Q3->FindObject("stats");
2253 Quartic_stats->SetX1NDC(0.2);
2254 Quartic_stats->SetX2NDC(0.5);
2255 Quartic_stats->SetY1NDC(.8);
2256 Quartic_stats->SetY2NDC(.9);
2258 TH1D *r3_clone=(TH1D*)r3_num_Q3->Clone();
2259 r3_clone->Fit(QuadraticFit,"IME","",0,0.1);
2261 TPaveStats *Quadratic_stats=(TPaveStats*)r3_clone->FindObject("stats");
2262 Quadratic_stats->SetX1NDC(0.55);
2263 Quadratic_stats->SetX2NDC(0.85);
2264 Quadratic_stats->SetY1NDC(.8);
2265 Quadratic_stats->SetY2NDC(.9);
2267 QuarticFit->Draw("same");
2268 QuadraticFit->Draw("same");
2269 Quartic_stats->Draw("same");
2270 Quadratic_stats->Draw("same");
2275 cout<<"r3 values:"<<endl;
2276 for(int ii=0; ii<10; ii++){
2277 cout<<r3_num_Q3->GetBinContent(ii+1)<<", ";
2280 cout<<"r3 errors:"<<endl;
2281 for(int ii=0; ii<10; ii++){
2282 cout<<r3_num_Q3->GetBinError(ii+1)<<", ";
2286 // M0,pi- base values with TTC=0.035 (old standard)
2287 //double r3base[10]={0, 2.02584, 1.82659, 1.85384, 1.84294, 1.82431, 1.72205, 1.5982, 1.17379, 0.881535};
2288 //double r3base_e[10]={0, 0.15115, 0.038311, 0.0231936, 0.0209217, 0.0291552, 0.0406093, 0.0624645, 0.0933104, 0.122683};
2289 // M0,pi+ base values with TTC=0.035 (old standard)
2290 //double r3base[10]={0, 1.84738, 1.81805, 1.80057, 1.82563, 1.76585, 1.6749, 1.37454, 1.37558, 0.74562};
2291 //double r3base_e[10]={0, 0.147285, 0.0381688, 0.0231643, 0.0209571, 0.0292255, 0.0407515, 0.0627228, 0.0937758, 0.123392};
2292 // M1,pi- base values with TTC=0.035 (old standard)
2293 //double r3base[10]={0, 1.72641, 1.82547, 1.79747, 1.84141, 1.83557, 1.7994, 1.67511, 1.44238, 1.03004};
2294 //double r3base_e[10]={0, 0.192385, 0.0471262, 0.0270385, 0.0229849, 0.0302041, 0.0397806, 0.0595701, 0.0902909, 0.123976};
2295 // M2,pi- base values with TTC=0.035 (old standard)
2296 //double r3base[10]={0, 1.61726, 1.58973, 1.7834, 1.8914, 1.75445, 1.76511, 1.81367, 1.55617, 1.17797};
2297 //double r3base_e[10]={0, 0.417924, 0.0973277, 0.0527296, 0.0421998, 0.0522792, 0.065144, 0.0923891, 0.135867, 0.190967};
2298 // M9,pi- base values with TTC=0.035 (old standard)
2299 //double r3base[10]={0, 2.51235, 2.68042, 1.85333, 1.66794, 1.39832, 1.717, 1.5337, 1.80686, 1.42286};
2300 //double r3base_e[10]={0, 3.09403, 0.650924, 0.272827, 0.160562, 0.145494, 0.137767, 0.146212, 0.164592, 0.189126};
2302 // M0 PosB, pi- base values with TTC=0.035 (old standard)
2303 //double r3base[10]={0, 1.81361, 1.84779, 1.81693, 1.8252, 1.78366, 1.71822, 1.41243, 1.31656, 0.719842};
2304 //double r3base_e[10]={0, 0.158656, 0.0405, 0.0244075, 0.0219925, 0.0306193, 0.0426494, 0.0655837, 0.0980259, 0.129011};
2305 // M1 PosB, pi- base values with TTC=0.035 (old standard)
2306 //double r3base[10]={0, 1.7208, 1.79016, 1.76858, 1.85008, 1.80416, 1.67274, 1.6501, 1.45098, 1.1052};
2307 //double r3base_e[10]={0, 0.202391, 0.0496228, 0.0282103, 0.0239023, 0.0313499, 0.041234, 0.061679, 0.0934039, 0.128249};
2309 // M1 PID1p5 special case (M0 MomRes used),pi- base values with TTC=0.035 (old standard)
2310 //double r3base[10]={0, 1.77035, 1.82471, 1.77934, 1.80615, 1.77819, 1.71105, 1.60505, 1.44563, 1.09473};
2311 //double r3base_e[10]={0, 0.193326, 0.0471058, 0.0269928, 0.0229581, 0.0301777, 0.0397553, 0.0595472, 0.0902768, 0.123968};
2313 // M0,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2314 //double r3base[10]={0, 1.93236, 1.77265, 1.71324, 1.63426, 1.56167, 1.38348, 1.25064, 0.897615, 0.691928};
2315 //double r3base_e[10]={0, 0.221991, 0.0433103, 0.023234, 0.0187678, 0.0243522, 0.0318668, 0.0459759, 0.066614, 0.0873847};
2316 // M1,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2317 //double r3base[10]={0, 1.46815, 1.76361, 1.68282, 1.6334, 1.56039, 1.44628, 1.27343, 1.01443, 0.691937};
2318 //double r3base_e[10]={0, 0.277164, 0.053134, 0.0271048, 0.0207043, 0.0253691, 0.0315642, 0.0443283, 0.0641466, 0.0858159};
2319 // M2,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2320 //double r3base[10]={0, 1.43331, 1.72759, 1.71155, 1.68957, 1.53846, 1.45209, 1.40698, 1.16202, 1.04672};
2321 //double r3base_e[10]={0, 0.594776, 0.109832, 0.0531013, 0.0383514, 0.0444564, 0.0525963, 0.0705969, 0.0998679, 0.135742};
2322 // M9,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2323 //double r3base[10]={0, 6.98181, 2.42845, 1.71962, 1.57641, 1.23102, 1.49543, 1.48063, 1.50938, 1.27714};
2324 //double r3base_e[10]={0, 7.09254, 0.731841, 0.277336, 0.152199, 0.130775, 0.118981, 0.121864, 0.132875, 0.148409};
2326 // M0,pi- base values with 3sigma
2327 //double r3base[10]={0, 1.71114, 1.72142, 1.74749, 1.69154, 1.60411, 1.47805, 1.27561, 1.02355, 0.769853};
2328 //double r3base_e[10]={0, 0.251731, 0.0398122, 0.0196354, 0.0153904, 0.0197609, 0.0258499, 0.0377931, 0.0560878, 0.0774703};
2329 // M1,pi- base values with 3sigma
2330 //double r3base[10]={0, 2.02471, 1.78737, 1.68681, 1.69639, 1.62634, 1.57337, 1.42229, 1.28352, 0.958121};
2331 //double r3base_e[10]={0, 0.314475, 0.048507, 0.022855, 0.0169939, 0.0207322, 0.0260424, 0.0375715, 0.0569927, 0.0830874};
2332 // M2,pi- base values with 3sigma
2333 //double r3base[10]={0, 1.50807, 1.75138, 1.71876, 1.76974, 1.6322, 1.5883, 1.5079, 1.4551, 1.49369};
2334 //double r3base_e[10]={0, 0.658863, 0.0979719, 0.0436736, 0.0308217, 0.035654, 0.0425953, 0.0585648, 0.0865817, 0.128056};
2335 // M9,pi- base values with 3sigma
2336 //double r3base[10]={0, 1.17509, 1.5117, 1.98348, 1.63032, 1.57195, 1.48828, 1.48218, 1.53689, 1.39533};
2337 //double r3base_e[10]={0, 4.89181, 0.60212, 0.230712, 0.11983, 0.102686, 0.0937373, 0.0971281, 0.107949, 0.123923};
2339 // M0,pi-,EDbin=0, Femto Plus+Minus
2340 //double r3base[10]={0, 1.77055, 1.75926, 1.7711, 1.67159, 1.57851, 1.41187, 1.11004, 0.870851, 0.442019};
2341 //double r3base_e[10]={0, 0.0688956, 0.0232182, 0.0156323, 0.0155991, 0.0225749, 0.0324584, 0.0516114, 0.0797203, 0.112197};
2342 // M1,pi-,EDbin=0, Femto Plus+Minus
2343 //double r3base[10]={0, 1.83169, 1.79042, 1.7275, 1.67116, 1.6032, 1.5423, 1.30426, 1.13238, 0.656902};
2344 //double r3base_e[10]={0, 0.0891247, 0.0286553, 0.0182713, 0.0171956, 0.0235918, 0.0324304, 0.0510863, 0.0823637, 0.125496};
2345 // M2,pi-,EDbin=0, Femto Plus+Minus
2346 //double r3base[10]={0, 1.60769, 1.81565, 1.75909, 1.74026, 1.62865, 1.58101, 1.50583, 1.3503, 1.27367};
2347 //double r3base_e[10]={0, 0.188422, 0.0581482, 0.0348698, 0.0309363, 0.0401393, 0.0519907, 0.0771148, 0.120583, 0.187679};
2348 // M3,pi-,EDbin=0, Femto Plus+Minus
2349 double r3base[10]={0, 1.52778, 1.81442, 1.78982, 1.69978, 1.70934, 1.59516, 1.56605, 1.45676, 1.17468};
2350 double r3base_e[10]={0, 0.288329, 0.0882351, 0.0510576, 0.0432015, 0.0536294, 0.0667149, 0.0950229, 0.143021, 0.217322};
2351 // M8,pi-,EDbin=0, Femto Plus+Minus
2352 //double r3base[10]={0, 2.27921, 2.40794, 1.72409, 1.66853, 1.70966, 1.52681, 1.64425, 1.18747, 1.35988};
2353 //double r3base_e[10]={0, 1.08315, 0.289559, 0.136951, 0.0919134, 0.0900869, 0.0901754, 0.101887, 0.121289, 0.148063};
2354 // M9,pi-,EDbin=0, Femto Plus+Minus
2355 //double r3base[10]={0, 2.1271, 2.12837, 1.95071, 1.4832, 1.50308, 1.41068, 1.29826, 1.10739, 0.755899};
2356 //double r3base_e[10]={0, 1.52325, 0.385705, 0.181328, 0.116165, 0.10896, 0.105705, 0.115163, 0.131994, 0.156299};
2360 double r3base_noMRC[10]={0, 4.16425, 0.429681, 1.56506, 1.64596, 1.73785, 1.57181, 1.51971, 1.59096, 1.54403};
2361 double r3base_noMRC_e[10]={0, 2.00855, 0.437928, 0.187872, 0.1087, 0.097367, 0.0893871, 0.0933358, 0.103603, 0.11839};
2362 // M9,pi+ (MRC normal)
2363 double r3base_MRC[10]={0, 3.84986, 0.558354, 1.60733, 1.67855, 1.75362, 1.59637, 1.54972, 1.62045, 1.57307};
2364 double r3base_MRC_e[10]={0, 1.93428, 0.402333, 0.171664, 0.0996222, 0.0903509, 0.0840615, 0.0889731, 0.100023, 0.115704};
2365 // M9,pi+ (MRC 50% larger)
2366 double r3base_MRC_over[10]={0, 3.70355, 0.614999, 1.62641, 1.69371, 1.76209, 1.60931, 1.56565, 1.63663, 1.58935};
2367 double r3base_MRC_over_e[10]={0, 1.90276, 0.387084, 0.164703, 0.0956857, 0.0872585, 0.0816823, 0.0870033, 0.0984005, 0.114503};
2371 TH1D *r3_noMRC=(TH1D*)r3_num_Q3->Clone();
2372 TH1D *r3_MRC=(TH1D*)r3_num_Q3->Clone();
2373 TH1D *r3_MRC_over=(TH1D*)r3_num_Q3->Clone();
2375 TH1D *r3Sys=(TH1D*)r3_num_Q3->Clone();
2376 for(int ii=1; ii<10; ii++){
2377 double Stat_RelativeError = sqrt(fabs(pow(r3_num_Q3->GetBinError(ii+1),2) - pow(r3base_e[ii],2)));// + for independent stat's
2378 r3Sys->SetBinContent(ii+1, (r3base[ii]-r3_num_Q3->GetBinContent(ii+1))/r3base[ii]);
2379 r3Sys->SetBinError(ii+1, Stat_RelativeError);
2381 r3_noMRC->SetBinContent(ii+1, r3base_noMRC[ii]);
2382 r3_MRC->SetBinContent(ii+1, r3base_MRC[ii]);
2383 r3_MRC_over->SetBinContent(ii+1, r3base_MRC_over[ii]);
2384 r3_noMRC->SetBinError(ii+1, r3base_noMRC_e[ii]);
2385 r3_MRC->SetBinError(ii+1, r3base_MRC_e[ii]);
2386 r3_MRC_over->SetBinError(ii+1, r3base_MRC_over_e[ii]);
2388 r3_noMRC->SetMarkerStyle(20); r3_noMRC->SetMarkerColor(1); r3_noMRC->SetLineColor(1);
2389 r3_MRC->SetMarkerStyle(23); r3_MRC->SetMarkerColor(4); r3_MRC->SetLineColor(4);
2390 r3_MRC_over->SetMarkerStyle(24); r3_MRC_over->SetMarkerColor(2); r3_MRC_over->SetLineColor(2);
2393 //r3_noMRC->Draw("same");
2394 //r3_MRC_over->Draw("same");
2395 TLegend *legendMRC = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2396 legendMRC->SetBorderSize(1);
2397 legendMRC->SetTextSize(.03);// small .03; large .06
2398 legendMRC->SetFillColor(0);
2399 //legendMRC->AddEntry(r3_MRC,"Momentum Resolution Corrected","p");
2400 //legendMRC->AddEntry(r3_noMRC,"No Correction","p");
2401 //legendMRC->AddEntry(r3_MRC_over,"50% Increased Correction","p");
2402 //legendMRC->Draw("same");
2404 r3Sys->GetYaxis()->SetTitle("(r_{3}^{t1}-r_{3}^{t2})/r_{3}^{t1}");
2405 r3Sys->GetYaxis()->SetTitleOffset(1.8);
2406 r3Sys->SetMinimum(-0.1);
2407 r3Sys->SetMaximum(0.1);
2408 r3Sys->GetXaxis()->SetRangeUser(0,0.099);
2410 TF1 *lineFit=new TF1("lineFit","pol0",0,0.02);
2411 gStyle->SetOptFit(111);
2412 r3Sys->Fit(lineFit,"MEQ","",0,0.1);
2417 ////////////////////////////////////////////////////////////////////////
2419 ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation
2421 TLegend *legendSTAR = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2422 legendSTAR->SetBorderSize(1);
2423 legendSTAR->SetTextSize(.03);// small .03; large .06
2424 legendSTAR->SetFillColor(0);
2426 TH1D *r3_STAR = new TH1D("r3_STAR","",Q3BINS,0,0.2);
2427 TH1D *r3_den_STAR = new TH1D("r3_den_STAR","",Q3BINS,0,0.2);
2428 double r3_STAR_num_e[20]={0};
2429 double r3_STAR_den_e[20]={0};
2431 for(int i=2; i<=20; i++){// bin number
2432 double Q12 = (i-0.5)*(0.005);
2433 int Q12bin = int(Q12/0.005)+1;
2435 for(int j=2; j<=20; j++){// bin number
2436 double Q13 = (j-0.5)*(0.005);
2437 int Q13bin = int(Q13/0.005)+1;
2439 for(int k=2; k<=20; k++){// bin number
2440 double Q23 = (k-0.5)*(0.005);
2441 int Q23bin = int(Q23/0.005)+1;
2443 double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
2444 int Q3bin = Cterm1_noMRC->GetXaxis()->FindBin(Q3);
2445 if(Q3bin>19) continue;
2447 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;
2448 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;
2450 double K2_12 = CoulCorr2(CP12, Q12);
2451 double K2_13 = CoulCorr2(CP13, Q13);
2452 double K2_23 = CoulCorr2(CP23, Q23);
2454 double r3num = Cterm1_noMRC->GetBinContent(Q3bin)/(K2_12*K2_13*K2_23) - 1;
2455 r3num -= C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1;
2456 r3num -= C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1;
2457 r3num -= C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1;
2458 double r3den = (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1);
2459 r3den *= (C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1);
2460 r3den *= (C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1);
2461 if(r3den<0) continue;
2462 r3den = sqrt(r3den);
2464 r3_STAR_num_e[Q3bin-1] += pow(Cterm1_noMRC->GetBinError(Q3bin)/(K2_12*K2_13*K2_23),2);
2465 r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q12bin)/(K2_12),2);
2466 r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q13bin)/(K2_13),2);
2467 r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q23bin)/(K2_23),2);
2468 r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q12bin)/(K2_12) * (C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1)*(C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1) /r3den,2);
2469 r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q13bin)/(K2_13) * (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1)*(C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1) /r3den,2);
2470 r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q23bin)/(K2_23) * (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1)*(C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1) /r3den,2);
2472 r3_STAR->Fill(Q3, r3num);
2473 r3_den_STAR->Fill(Q3, r3den);
2479 for(int i=0; i<20; i++){
2480 r3_STAR->SetBinError(i+1, sqrt( r3_STAR_num_e[i]));
2481 r3_den_STAR->SetBinError(i+1, sqrt( r3_STAR_den_e[i]));
2483 r3_STAR->Divide(r3_den_STAR);
2484 r3_STAR->SetMaximum(5);
2485 r3_STAR->SetMinimum(-10);
2486 r3_STAR->GetXaxis()->SetRangeUser(0,0.099);
2487 r3_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
2488 r3_STAR->GetYaxis()->SetTitle("r_{3}*");
2489 r3_STAR->SetLineColor(2); r3_STAR->SetMarkerColor(2);
2490 r3_STAR->SetMarkerStyle(kFullStar);
2491 r3_STAR->SetMarkerSize(1.2);
2493 r3_STAR->Fit(QuarticFit,"IME","",0.01,0.1);
2495 TPaveStats *Quartic_stats_STAR=(TPaveStats*)r3_STAR->FindObject("stats");
2496 Quartic_stats_STAR->SetX1NDC(0.2);
2497 Quartic_stats_STAR->SetX2NDC(0.5);
2498 Quartic_stats_STAR->SetY1NDC(.8);
2499 Quartic_stats_STAR->SetY2NDC(.9);
2501 TH1D *r3_clone_STAR=(TH1D*)r3_STAR->Clone();
2502 r3_STAR->Fit(QuadraticFit,"IME","",0,0.1);
2504 TPaveStats *Quadratic_stats_STAR=(TPaveStats*)r3_clone_STAR->FindObject("stats");
2505 Quadratic_stats_STAR->SetX1NDC(0.55);
2506 Quadratic_stats_STAR->SetX2NDC(0.85);
2507 Quadratic_stats_STAR->SetY1NDC(.8);
2508 Quadratic_stats_STAR->SetY2NDC(.9);
2510 QuarticFit->Draw("same");
2511 QuadraticFit->Draw("same");
2512 Quartic_stats_STAR->Draw("same");
2513 Quadratic_stats_STAR->Draw("same");
2514 //Cterm1_noMRC->Draw();
2516 /////////////////////////////////////////////////////////////////////////
2517 if(SameCharge) r3_num_Q3->Draw("same");
2518 legendSTAR->AddEntry(r3_STAR,"STAR method","p");
2519 legendSTAR->AddEntry(r3_num_Q3,"ALICE method","p");
2520 legendSTAR->Draw("same");
2524 //c3_hist_STAR->Draw();
2525 //Cterm1_noMRC->Draw();
2527 //r3_num_Q3->SetMarkerStyle(20);
2528 //r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
2529 //r3_num_Q3->SetMinimum(-7);
2530 //r3_num_Q3->SetMaximum(2.2);
2534 double C3_star_noCoul[20]={0, 1.59478, 1.49099, 1.4239, 1.31166, 1.21337, 1.14344, 1.09671, 1.06525, 1.04391, 1.0306, 1.02045, 1.01359, 1.0035, 0.992264, 0.989587, 0.989968, 0.996795, 0, 0};
2535 double C3_star_noCoul_e[20]={0, 0.260834, 0.0191839, 0.00524928, 0.0021307, 0.00110448, 0.000661211, 0.000438894, 0.000311247, 0.00023452, 0.000182819, 0.000147496, 0.000126547, 0.000130428, 0.000159494, 0.000240658, 0.000489113, 0.00387954, 0, 0};
2536 TH1D *C3_star=(TH1D*)Cterm1->Clone();
2537 C3_star->Add(C3_star,-1);
2538 for(int i=0; i<20; i++) {C3_star->SetBinContent(i+1, C3_star_noCoul[i]); C3_star->SetBinError(i+1, C3_star_noCoul_e[i]);}
2539 C3_star->SetMarkerColor(4);
2540 C3_star->SetLineColor(4);
2541 //C3_star->Draw("same");
2542 //legend1->AddEntry(C3_star,"---; 200 GeV","p");
2544 //legend1->Draw("same");
2546 //TGraph *gr_MomRes_term1=new TGraph(10,
2553 //////////////////////////////////////////////////////////////////////////////////////
2561 TString *savefilename = new TString("OutFiles/OutFile");
2562 if(SameCharge) savefilename->Append("SC");
2563 else savefilename->Append("MC");
2565 if(CHARGE==1) savefilename->Append("Pos");
2566 else savefilename->Append("Neg");
2568 if(IncludeG) savefilename->Append("YesG");
2569 else savefilename->Append("NoG");
2571 if(IncludeEWfromTherm) savefilename->Append("EWfromTherm");
2572 if(IncludeEW) savefilename->Append("EW");
2573 if(GRS) savefilename->Append("GRS");
2574 else savefilename->Append("Omega0");
2576 if(EDbin==0) savefilename->Append("Kt3_1");
2577 else savefilename->Append("Kt3_2");
2579 savefilename->Append("_Kt");
2580 *savefilename += Ktbin;
2581 savefilename->Append("_M");
2582 *savefilename += Mbin;
2583 savefilename->Append(".root");
2584 TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
2585 can1->Write("can1");
2586 C2_ss->Write("C2_ss");
2587 C2_os->Write("C2_os");
2588 C2_ss_Momsys->Write("C2_ss_Momsys");
2589 C2_os_Momsys->Write("C2_os_Momsys");
2590 C2_ss_Ksys->Write("C2_ss_Ksys");
2591 C2_os_Ksys->Write("C2_os_Ksys");
2592 fitC2ss->Write("fitC2ss");
2593 fitC2os->Write("fitC2os");
2594 if(IncludeEWfromTherm) {
2595 fitC2ssTherm->Write("fitC2ssTherm");
2596 fitC2osTherm->Write("fitC2osTherm");
2597 fitC2ssThermGaus->Write("fitC2ssThermGaus");
2598 fitC2osThermGaus->Write("fitC2osThermGaus");
2600 C2Therm_ss->Write("C2Therm_ss");
2601 C2Therm_os->Write("C2Therm_os");
2602 K2_ss->Write("K2_ss");
2603 K2_os->Write("K2_os");
2604 Cterm1->Write("C3");
2605 Coul_Omega0->Write("Coul_Omega0");
2606 Coul_GRiverside->Write("Coul_GRiverside");
2607 GenSignalExpected_num->Write("C3_EWexpectation");
2608 num_QS->Write("num_QS");
2609 c3_hist->Write("c3");
2610 Combinatorics_1d->Write("Combinatorics_1d");
2611 ChiSquaredNDF->Write("ChiSquaredNDF");
2612 r3_num_Q3->Write("r3_Q3");
2621 void ReadCoulCorrections(int ST, int RVal, int bVal, int kt){
2622 ///////////////////////
2625 if(FileSetting!=6) fname = new TString("KFile.root");
2626 if(FileSetting==6) fname = new TString("KFile_Gauss.root");
2628 TFile *File=new TFile(fname->Data(),"READ");
2629 if(ST==0){// Gaussian
2630 if(RVal < 3 || RVal > 10) cout<<"Coulomb Correlation Gaussian radius outside of range!!!!!!!!!!!!!!!!"<<endl;
2631 TH2D *tempG_ss = (TH2D*)File->Get("K2ssG");
2632 CoulCorr2SS = (TH1D*)tempG_ss->ProjectionY("CoulCorr2SS",RVal-2, RVal-2);
2633 TH2D *tempG_os = (TH2D*)File->Get("K2osG");
2634 CoulCorr2OS = (TH1D*)tempG_os->ProjectionY("CoulCorr2OS",RVal-2, RVal-2);
2636 if(ST==1){//Therminator
2637 if(bVal!=2 && bVal!=3 && bVal!=5 && bVal!=7 && bVal!=8 && bVal!=9) cout<<"Therminator bVal not acceptable in 2-particle Coulomb read"<<endl;
2639 if(kt==10){// kt integrated
2640 TH2D *tempT_ss = (TH2D*)File->Get("K2ssT");
2641 TH2D *tempT_os = (TH2D*)File->Get("K2osT");
2642 CoulCorr2SS = (TH1D*)tempT_ss->ProjectionY("CoulCorr2SS",bBin, bBin);
2643 CoulCorr2OS = (TH1D*)tempT_os->ProjectionY("CoulCorr2OS",bBin, bBin);
2645 if(kt < 1 || kt > 6) cout<<"kt bin out of range in 2-particle Coulomb read"<<endl;
2646 TH3D *tempT3_ss = (TH3D*)File->Get("K2ssT_kt");
2647 TH3D *tempT3_os = (TH3D*)File->Get("K2osT_kt");
2648 CoulCorr2SS = (TH1D*)tempT3_ss->ProjectionZ("CoulCorr2SS",bBin, bBin, kt,kt);
2649 CoulCorr2OS = (TH1D*)tempT3_os->ProjectionZ("CoulCorr2OS",bBin, bBin, kt,kt);
2652 CoulCorr2SS->SetDirectory(0);
2653 CoulCorr2OS->SetDirectory(0);
2656 double CoulCorr2(int chargeproduct, double Q2){// returns K2
2662 indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2665 if(indexH >= Nlines) return 1.0;
2666 if(chargeproduct==1){
2667 slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
2668 slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
2669 return slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1);
2671 slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1));
2672 slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1));
2673 return slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
2677 int Qbin=CoulCorr2SS->GetXaxis()->FindBin(Q2);
2678 if(chargeproduct==1) return CoulCorr2SS->GetBinContent(Qbin);
2679 else return CoulCorr2OS->GetBinContent(Qbin);
2683 //----------------------------------------------------------------------
2686 //________________________________________________________________________
2687 void fcnC2_Global(int& npar, double* deriv, double& f, double par[], int flag){
2691 double Rch=par[3]/FmToGeV;
2692 double Rcoh=par[4]/FmToGeV;
2694 double CSS=0, COS=0;
2697 double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2698 NFitPoints_C2global=0;
2699 if(LinkN) par[9]=par[0];// Link N factors
2701 for(int i=1; i<BINRANGE_C2global; i++){
2703 qinvSS = AvgQinvSS[i];
2705 if(!GofP) Dp=fabs(par[2])/(1-fabs(par[2]));// p independent
2706 else Dp = fabs(par[2])/(1-fabs(par[2]));
2708 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2) - pow(RT,2))*pow(AvgP[Ktbin_GofP-1][i]-qinvSS/2.,2));
2709 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2) - pow(RT,2))*pow(AvgP[Ktbin_GofP-1][i]+qinvSS/2.,2));
2710 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][i]-qinvSS/2.,2) - 2*MT/Temp);
2711 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][i]+qinvSS/2.,2) - 2*MT/Temp);
2712 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]-qinvSS/2.,2) - 2*MT/Temp);
2713 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]+qinvSS/2.,2) - 2*MT/Temp);
2714 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2715 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2716 double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2717 double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2719 if(!GofP) {Dp1=Dp; Dp2=Dp;}
2722 double arg=qinvSS*Rch;
2723 EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2724 EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2725 EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2726 EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2728 double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
2729 double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
2730 CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
2731 if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2732 else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
2733 //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
2735 CSS *= par[1]*K2SS[i];
2740 if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2741 COS *= par[1]*K2OS[i];
2743 COS *= par[9];// different Norm
2745 if(C2ssFitting[i] > 0){
2747 double errorSS = sqrt(pow( (CSS/par[0] - (1-par[1]))/K2SS[i] * K2SS_e[i] * par[0],2) + pow(C2ssFitting_e[i],2));
2748 //double errorSS = C2ssFitting_e[i];
2749 SumChi2 += pow((CSS - C2ssFitting[i])/errorSS,2);
2750 NFitPoints_C2global++;
2754 double errorOS = sqrt(pow( (COS/par[9] - (1-par[1]))/K2OS[i] * K2OS_e[i] * par[9],2) + pow(C2osFitting_e[i],2));
2755 //double errorOS = C2osFitting_e[i];
2756 SumChi2 += pow((COS - C2osFitting[i])/errorOS,2);
2757 NFitPoints_C2global++;
2766 //________________________________________________________________________
2767 double C2ssFitFunction(Double_t *x, Double_t *par)
2769 double Rch=par[3]/FmToGeV;
2770 double Rcoh=par[4]/FmToGeV;
2772 int qbin = int(fabs(x[0]/0.005));
2773 if(qbin >= BINRANGE_C2global) return 1.0;
2774 double qinvSS = AvgQinvSS[qbin];
2775 double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2777 if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
2778 else Dp = fabs(par[2])/(1-fabs(par[2]));
2780 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]-qinvSS/2.,2));
2781 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]+qinvSS/2.,2));
2782 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][qbin]-qinvSS/2.,2) - 2*MT/Temp);
2783 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][qbin]+qinvSS/2.,2) - 2*MT/Temp);
2784 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2785 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2786 double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2787 double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2789 if(!GofP) {Dp1=Dp; Dp2=Dp;}
2790 double arg=qinvSS*Rch;
2791 double EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2792 EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2793 EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2794 EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2795 double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
2796 double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
2797 double CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
2798 if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2799 else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
2800 //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
2801 CSS *= par[1]*K2SS[qbin];
2805 //cout<<qinvSS<<" "<<Dp1/(1+Dp1) * Dp2/(1+Dp2)<<endl;
2809 //________________________________________________________________________
2810 double C2osFitFunction(Double_t *x, Double_t *par)
2812 if(LinkN) par[9]=par[0];// Link N factors
2813 double Rcoh=par[4]/FmToGeV;
2815 int qbin = int(fabs(x[0]/0.005));
2816 if(qbin >= BINRANGE_C2global) return 1.0;
2817 double qinvOS = AvgQinvOS[qbin];
2819 if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
2820 else Dp = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin],2));
2821 //Dp = fabs(par[2])/(1-fabs(par[2]));
2822 double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]-qinvOS/2.,2));
2823 double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]+qinvOS/2.,2));
2824 //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
2825 //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
2827 if(!GofP) {Dp1=Dp; Dp2=Dp;}
2829 if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2830 COS *= par[1]*K2OS[qbin];
2835 //________________________________________________________________________
2836 void fcnC3_Global(int& npar, double* deriv, double& f, double par[], int flag){
2842 //double D=fabs(par[2])/(1-fabs(par[2]));
2843 double Rch=par[3]/FmToGeV;
2844 //double Rcoh=par[4]/FmToGeV;
2847 double EW_12=0, EW_13=0, EW_23=0;
2848 NFitPoints_C3global=0;
2850 for(int i=0; i<BINRANGE_C2global; i++){
2851 qinvSS_12 = AvgQinvSS[i];
2852 double arg_12=qinvSS_12*Rch;
2853 EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
2854 EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
2855 EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
2856 EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
2858 for(int j=0; j<BINRANGE_C2global; j++){
2859 qinvSS_13 = AvgQinvSS[j];
2860 double arg_13=qinvSS_13*Rch;
2861 EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
2862 EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
2863 EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
2864 EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
2866 for(int k=0; k<BINRANGE_C2global; k++){
2867 qinvSS_23 = AvgQinvSS[k];
2868 double arg_23=qinvSS_23*Rch;
2869 EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
2870 EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
2871 EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
2872 EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
2874 if(i==0 || j==0 || k==0) continue;
2875 CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
2876 CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
2877 CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
2878 CSS += 2*EW_12*EW_13*EW_23 * exp(-1/2.*pow(Rch,2)*(pow(qinvSS_12,2)+pow(qinvSS_13,2)+pow(qinvSS_23,2)));
2881 /*double errorSS = sqrt(pow(C3ssFitting_e[i],2));
2883 SumChi2 += pow((CSS - C3ssFitting[i])/errorSS,2);
2884 NFitPoints_C3global++;
2888 //double errorOS = sqrt(pow(C3osFitting_e[i],2));
2889 //SumChi2 += pow((COS - C3osFitting[i])/errorOS,2);
2890 //NFitPoints_C3global++;
2902 //______________________________________________________________________
2903 /*void C3ssFitFunction(Double_t *x, Double_t *par){
2905 double qinvSS_12=0, qinvOS_12=0;
2906 double qinvSS_13=0, qinvOS_13=0;
2907 double qinvSS_23=0, qinvOS_23=0;
2909 double D=fabs(par[2])/(1-fabs(par[2]));
2910 double Rch=par[3]/FmToGeV;
2911 double Rcoh=par[4]/FmToGeV;
2912 double CSS=0, COS=0;
2913 double term1=0, term2=0;
2915 double EW_12=0, EW_13=0, EW_23=0;
2916 NFitPoints_C3global=0;
2919 int qbin_12 = x[0]*1000/5;
2920 qinvSS_12 = AvgQinvSS[qbin_12];
2921 qinvOS_12 = AvgQinvOS[qbin_12];
2922 double arg_12=qinvSS_12*Rch;
2923 EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
2924 EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
2925 EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
2926 EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
2928 int qbin_13 = x[1]*1000/5;
2929 qinvSS_13 = AvgQinvSS[qbin_13];
2930 qinvOS_13 = AvgQinvOS[qbin_13];
2931 double arg_13=qinvSS_13*Rch;
2932 EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
2933 EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
2934 EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
2935 EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
2937 int qbin_23 = x[2]*1000/5;
2938 qinvSS_23 = AvgQinvSS[qbin_23];
2939 qinvOS_23 = AvgQinvOS[qbin_23];
2940 double arg_23=qinvSS_23*Rch;
2941 EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
2942 EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
2943 EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
2944 EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
2946 if(qbin_12==0 || qbin_13==0 || qbin_23==0) return 0;
2947 CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
2948 CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
2949 CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
2950 CSS += 2*EW_12*EW_13*EW_23 * exp(-1/2.*pow(Rch,2)*(pow(qinvSS_12,2)+pow(qinvSS_13,2)+pow(qinvSS_23,2)));
2955 //----------------------------------------------------------------------
2956 void fcnOSL(int& npar, double* deriv, double& f, double par[], int flag){
2958 double qout=0, qside=0, qlong=0;
2965 for(int i=0; i<BINS_OSL-10; i++){// out, 0 to BINS_OSL-10
2966 for(int j=0; j<BINS_OSL-10; j++){// side, 0 to BINS_OSL-10
2967 for(int k=0; k<BINS_OSL-10; k++){// long, 0 to BINS_OSL-10
2969 if(B[i][j][k] < 1. ) continue;
2971 qout = (i+0.5)*0.005;
2972 qside = (j+0.5)*0.005;
2973 qlong = (k+0.5)*0.005;
2975 Gaus_ch = exp(-pow(par[3]*qout/FmToGeV,2) - pow(par[4]*qside/FmToGeV,2) - pow(par[5]*qlong/FmToGeV,2));
2977 //C = par[0]*(1 + par[1]*(K_OSL[i][j][k]-1) + K_OSL[i][j][k]*par[1]*Gaus_ch);
2978 C = par[0]*( (1-par[1]) + par[1]*K_OSL[i][j][k]*(1+Gaus_ch));
2982 double error = pow(A_e[i][j][k]/B[i][j][k],2) + pow(B_e[i][j][k]*A[i][j][k]/pow(B[i][j][k],2),2);
2983 error += pow((C/par[0] - (1-par[1]))/K_OSL[i][j][k] * K_OSL_e[i][j][k]*par[0],2);
2984 error = sqrt(error);
2985 SumChi2 += pow( (C - A[i][j][k]/B[i][j][k])/error,2);
2988 if(A[i][j][k] >= 1) term1 = C*(A[i][j][k]+B[i][j][k])/(A[i][j][k]*(C+1));
2990 term2 = (A[i][j][k]+B[i][j][k])/(B[i][j][k]*(C+1));
2992 if(term1 > 0.0 && term2 > 0.0){
2993 lnL += A[i][j][k]*log(term1) + B[i][j][k]*log(term2);
2994 }else if(term1==0 && term2 > 0.0){
2995 lnL += B[i][j][k]*log(term2);
2996 }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
3007 //________________________________________________________________________
3008 double OSLfitfunction(Double_t *x, Double_t *par)
3010 int bin_out = int(1000*(x[0])/5.);
3011 int bin_side = int(1000*(x[1])/5.);
3012 int bin_long = int(1000*(x[2])/5.);
3014 double K = CoulCorr2(+1, avg_q[bin_out][bin_side][bin_long]);
3015 double Gaus_ch = exp(-pow(par[3]*x[0]/FmToGeV,2) - pow(par[4]*x[1]/FmToGeV,2) - pow(par[5]*x[2]/FmToGeV,2));
3016 return par[0]*(1 + par[1]*(K-1) + K*par[1]*Gaus_ch);
3019 //__________________________________________________________________________
3021 void fcn_3(int& npar, double* deriv, double& f, double par[], int flag){
3023 double q12=0, q13=0, q23=0;
3024 double K12=0, K13=0, K23=0, K123=0;
3025 double G1=0, G2=0, G3=0, G4=0;
3026 double EW12=0, EW13=0, EW23=0;
3028 double term1=0, term2=0;
3031 // start from 2-4 MeV bins
3032 for(int i=1; i<BINRANGE_3-10; i++){// q12
3033 for(int j=1; j<BINRANGE_3-10; j++){// q13
3034 for(int k=1; k<BINRANGE_3-10; k++){// q23
3036 if(B_3[i][j][k] < 1. ) continue;
3038 q12 = (i+0.5)*0.002;
3039 q13 = (j+0.5)*0.002;
3040 q23 = (k+0.5)*0.002;
3042 //K12 = 1/CoulCorr_CsorgoMate(2, q12);
3043 //K13 = 1/CoulCorr_CsorgoMate(2, q13);
3044 //K23 = 1/CoulCorr_CsorgoMate(2, q23);
3045 //K12 = 1/Coul_pipi_data(q12);
3046 //K13 = 1/Coul_pipi_data(q13);
3047 //K23 = 1/Coul_pipi_data(q23);
3048 K12=1.0; K13=1.0; K23=1.0;
3051 //K123 = 1 - pow(par[1],3)*(1-K12*K13*K23);
3052 //K12 = 1 - pow(par[1],3)*(1-K12);
3053 //K13 = 1 - pow(par[1],3)*(1-K13);
3054 //K23 = 1 - pow(par[1],3)*(1-K23);
3056 EW12 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q12*par[2]/FmToGeV,3) - 12*pow(q12*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q12*par[2]/FmToGeV,4) -48*pow(q12*par[2]/FmToGeV,2) + 12);
3057 EW13 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q13*par[2]/FmToGeV,3) - 12*pow(q13*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q13*par[2]/FmToGeV,4) -48*pow(q13*par[2]/FmToGeV,2) + 12);
3058 EW23 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q23*par[2]/FmToGeV,3) - 12*pow(q23*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q23*par[2]/FmToGeV,4) -48*pow(q23*par[2]/FmToGeV,2) + 12);
3059 //EW12=1; EW13=1; EW23=1;
3061 G1 = K12*exp(-pow(par[2]*q12/FmToGeV,2)/2.)*EW12;
3062 G1 += K13*exp(-pow(par[2]*q13/FmToGeV,2)/2.)*EW13;
3063 G1 += K23*exp(-pow(par[2]*q23/FmToGeV,2)/2.)*EW23;
3065 G2 = K12*exp(-pow(par[2]*q12/FmToGeV,2))*pow(EW12,2);
3066 G2 += K13*exp(-pow(par[2]*q13/FmToGeV,2))*pow(EW13,2);
3067 G2 += K23*exp(-pow(par[2]*q23/FmToGeV,2))*pow(EW23,2);
3069 G3 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q23*q23)/2.)*EW12*EW23;
3070 G3 += exp(-pow(par[2]/FmToGeV,2)*(q13*q13 + q23*q23)/2.)*EW13*EW23;
3071 G3 += exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13)/2.)*EW12*EW13;
3074 G4 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13 + q23*q23)/2.);
3075 G4 *= EW12*EW13*EW23;
3078 C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3079 //C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3080 C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3081 C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3082 C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3083 C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3090 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));
3092 term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1));
3094 if(term1 > 0.0 && term2 > 0.0){
3095 lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2);
3096 }else if(term1==0 && term2 > 0.0){
3097 lnL += B_3[i][j][k]*log(term2);
3098 }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
3107 //________________________________________________________________________
3108 double D3fitfunction_3(Double_t *x, Double_t *par)
3110 double K12=CoulCorr2(+1, x[0]);
3111 double K13=CoulCorr2(+1, x[1]);
3112 double K23=CoulCorr2(+1, x[2]);
3113 K12=1.0; K13=1.0; K23=1.0;
3114 double K123=K12*K13*K23;
3116 //K123 = 1 - pow(par[1],3)*(1 - K123);
3117 //K12 = 1 - pow(par[1],3)*(1-K12);
3118 //K13 = 1 - pow(par[1],3)*(1-K13);
3119 //K23 = 1 - pow(par[1],3)*(1-K23);
3121 double EW12 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[0]*par[2]/FmToGeV,3) - 12*pow(x[0]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[0]*par[2]/FmToGeV,4) -48*pow(x[0]*par[2]/FmToGeV,2) + 12);
3122 double EW13 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[1]*par[2]/FmToGeV,3) - 12*pow(x[1]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[1]*par[2]/FmToGeV,4) -48*pow(x[1]*par[2]/FmToGeV,2) + 12);
3123 double EW23 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[2]*par[2]/FmToGeV,3) - 12*pow(x[2]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[2]*par[2]/FmToGeV,4) -48*pow(x[2]*par[2]/FmToGeV,2) + 12);
3124 //EW12=1; EW13=1; EW23=1;
3126 double G1=0, G2=0, G3=0, G4=0;
3128 G1 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2)/2.)*EW12;
3129 G1 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2)/2.)*EW13;
3130 G1 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2)/2.)*EW23;
3132 G2 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2))*pow(EW12,2);
3133 G2 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2))*pow(EW13,2);
3134 G2 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2))*pow(EW23,2);
3136 G3 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[2]*x[2])/2.)*EW12*EW23;
3137 G3 += exp(-pow(par[2]/FmToGeV,2)*(x[1]*x[1] + x[2]*x[2])/2.)*EW13*EW23;
3138 G3 += exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1])/2.)*EW12*EW13;
3141 G4 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/2.);
3142 G4 *= EW12*EW13*EW23;
3145 double C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3146 //double C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3147 C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3148 C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3149 C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3150 C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3156 void ReadCoulCorrections_Omega0(){
3157 // read in 3d 3-particle coulomb correlations = K3
3160 if(FileSetting!=6) coulfile = new TFile("KFile.root","READ");
3161 if(FileSetting==6) coulfile = new TFile("KFile_Gauss.root","READ");
3163 TString *name=new TString("K3ss_");
3165 CoulCorrOmega0SS = (TH3D*)coulfile->Get(name->Data());
3166 CoulCorrOmega0SS->SetDirectory(0);
3167 name=new TString("K3os_");
3169 CoulCorrOmega0OS = (TH3D*)coulfile->Get(name->Data());
3170 CoulCorrOmega0OS->SetDirectory(0);
3174 double CoulCorrGRS(bool SC, double Q_12, double Q_13, double Q_23){
3175 /*int Q12bin = CoulCorr2SS->GetXaxis()->FindBin(Q_12);
3176 int Q13bin = CoulCorr2SS->GetXaxis()->FindBin(Q_13);
3177 int Q23bin = CoulCorr2SS->GetXaxis()->FindBin(Q_23);*/
3178 int index12L = int(fabs(Q_12 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3179 int index12H = index12L+1;
3180 int index13L = int(fabs(Q_13 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3181 int index13H = index13L+1;
3182 int index23L = int(fabs(Q_23 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3183 int index23H = index23L+1;
3186 // Tricubic Interpolation
3187 double arr[4][4][4]={{{0}}};
3188 for(int x=0; x<4; x++){
3189 for(int y=0; y<4; y++){
3190 for(int z=0; z<4; z++){
3192 arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2SS->GetBinContent((index23L)+y)*CoulCorr2SS->GetBinContent((index13L)+z);
3194 arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2OS->GetBinContent((index23L)+y)*CoulCorr2OS->GetBinContent((index13L)+z);
3200 return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3202 // Trilinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3204 double xd = (Q_12-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3205 xd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3206 double yd = (Q_13-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3207 yd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index13H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3208 double zd = (Q_23-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3209 zd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index23H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3212 double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3213 c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3214 double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3215 c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3216 double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3217 c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3218 double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3219 c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3221 double c0 = c00*(1-yd) + c10*yd;
3222 double c1 = c01*(1-yd) + c11*yd;
3223 return (c0*(1-zd) + c1*zd);
3225 double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3226 c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3227 double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3228 c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3229 double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3230 c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3231 double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3232 c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3234 double c0 = c00*(1-yd) + c10*yd;
3235 double c1 = c01*(1-yd) + c11*yd;
3236 return (c0*(1-zd) + c1*zd);
3241 double CoulCorrOmega0(bool SC, double Q_12, double Q_13, double Q_23){
3242 int Q12bin = CoulCorrOmega0SS->GetXaxis()->FindBin(Q_12);
3243 int Q13bin = CoulCorrOmega0SS->GetZaxis()->FindBin(Q_13);
3244 int Q23bin = CoulCorrOmega0SS->GetYaxis()->FindBin(Q_23);
3245 int index12L = int(fabs(Q_12 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3246 int index12H = index12L+1;
3247 int index13L = int(fabs(Q_13 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3248 int index13H = index13L+1;
3249 int index23L = int(fabs(Q_23 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3250 int index23H = index23L+1;
3253 if(CoulCorrOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
3256 //////////////////////////
3257 // Tricubic Interpolation
3258 double arr[4][4][4]={{{0}}};
3259 for(int x=0; x<4; x++){
3260 for(int y=0; y<4; y++){
3261 for(int z=0; z<4; z++){
3262 arr[x][y][z] = CoulCorrOmega0SS->GetBinContent((index12L)+x, (index23L)+y, (index13L)+z);
3266 return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3268 ///////////////////////////
3269 // Trilinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3271 double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3272 xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3273 double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3274 yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3275 double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3276 zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3278 // interpolate along x
3279 double c00=0, c10=0, c01=0, c11=0;
3280 if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
3281 if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
3282 if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
3283 if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
3286 if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
3287 if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
3289 if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
3291 cout<<"Not all Omega0 bins well defined. Use GRS instead"<<endl;
3292 return CoulCorrGRS(kTRUE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
3296 }else {// mixed charge. Q12bin always designated as the same-charge pair
3297 if(CoulCorrOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
3300 //////////////////////////
3301 // Tricubic Interpolation
3302 double arr[4][4][4]={{{0}}};
3303 for(int x=0; x<4; x++){
3304 for(int y=0; y<4; y++){
3305 for(int z=0; z<4; z++){
3306 arr[x][y][z] = CoulCorrOmega0OS->GetBinContent((Q12bin)+x, (Q23bin)+y, (Q13bin)+z);
3310 return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3312 ///////////////////////////
3313 // TriLinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3315 double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3316 xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3317 double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3318 yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3319 double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3320 zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3321 // interpolate along x
3322 double c00=0, c10=0, c01=0, c11=0;
3323 if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
3324 if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
3325 if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
3326 if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
3329 if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
3330 if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
3332 if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
3334 cout<<"Not all Omega0 bins well defined. Use GRS instead"<<endl;
3335 return CoulCorrGRS(kFALSE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
3342 double Gamov(int chargeProduct, double qinv){
3344 double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
3346 return arg/(exp(arg)-1);
3348 /*double LednickyCorr(double qinv){
3350 int indexL = int(fabs(1000*qinv-Lednicky_qinv[0])/(2*2.027027));// starts from 2.027 MeV
3351 int indexH = indexL+1;
3352 if( (indexL >= 73)) return 1.0;
3353 double slope = (Lednicky_CoulStrong[indexL] - Lednicky_CoulStrong[indexH])/(Lednicky_qinv[indexL]-Lednicky_qinv[indexH]);
3354 return (slope*(1000*qinv - Lednicky_qinv[indexL]) + Lednicky_CoulStrong[indexL]);
3357 void ReadLednickyFile(int Rvalue){
3359 TString *filename=new TString("../Strong_FSI/converted_CoulStrong_");
3360 *filename += Rvalue;
3361 filename->Append("fm.root");
3362 TFile *Ledfile = new TFile(filename->Data(),"READ");
3363 TH1F *LedHisto = (TH1F*)Ledfile->Get("h27");
3364 for(int i=0; i<74; i++){
3365 Lednicky_qinv[i] = 2.*(LedHisto->GetXaxis()->GetBinLowEdge(i+1) + LedHisto->GetXaxis()->GetBinUpEdge(i+1))/2.;
3366 Lednicky_CoulStrong[i] = LedHisto->GetBinContent(i+1);
3370 void ReadMomResFile(int rvalue, double lambda){
3372 for(int i=0; i<40; i++){
3373 MomRes_C2_pp[i] = 1.;
3374 MomRes_C2_mp[i] = 1.;
3375 MomRes_term1_pp[i] = 1.;
3376 MomRes_term2_pp[i] = 1.;
3377 MomRes_term1_mp[i] = 1.;
3378 MomRes_term2_mp[i] = 1.;
3380 MomResDev_C2_pp[i] = 1.;
3381 MomResDev_C2_mp[i] = 1.;
3382 MomResDev_term1_pp[i] = 1.;
3383 MomResDev_term2_pp[i] = 1.;
3387 TString *momresfilename = new TString("MomResFile");
3388 momresfilename->Append("_Offline_");
3389 if(FileSetting==3) momresfilename->Append("PID1p5");
3390 else momresfilename->Append("11h");
3391 momresfilename->Append(".root");// no corresponding file for 10h
3393 TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
3395 TH2D *MomResWeights_pp = (TH2D*)MomResFile->Get("MomResHisto_pp");
3396 TH2D *MomResWeights_mp = (TH2D*)MomResFile->Get("MomResHisto_mp");
3397 TH2D *MomResWeights_pp_term1 = (TH2D*)MomResFile->Get("MomResHisto_pp_term1");
3398 TH2D *MomResWeights_pp_term2 = (TH2D*)MomResFile->Get("MomResHisto_pp_term2");
3399 TH2D *MomResWeights_mp_term1 = (TH2D*)MomResFile->Get("MomResHisto_mp_term1");
3400 TH2D *MomResWeights_mp_term2 = (TH2D*)MomResFile->Get("MomResHisto_mp_term2");
3403 TString *names3[2][5];// SC/MC, term#
3404 TString *names1D[2][5];// SC/MC, term#
3405 TString *names3_AvgK3[2];// SC/MC
3406 for(int ChProd=0; ChProd<2; ChProd++){
3407 if(ChProd==0) names3_AvgK3[ChProd] = new TString("AvgK3_SC_M");
3408 else names3_AvgK3[ChProd] = new TString("AvgK3_MC_M");
3409 *names3_AvgK3[ChProd] += MomResCentBin-1;
3410 //AvgK3[ChProd] = (TH3D*)MomResFile->Get(names3_AvgK3[ChProd]->Data());
3411 //AvgK3[ChProd]->SetDirectory(0);
3413 for(int term=0; term<5; term++){
3415 if(ChProd==0) {names3[ChProd][term] = new TString("MomResHisto3_SC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_SC_term");}
3416 else {names3[ChProd][term] = new TString("MomResHisto3_MC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_MC_term");}
3417 *names3[ChProd][term] += term+1;
3418 *names1D[ChProd][term] += term+1;
3419 names3[ChProd][term]->Append("_M");
3420 names1D[ChProd][term]->Append("_M");
3421 *names3[ChProd][term] += MomResCentBin-1;
3422 *names1D[ChProd][term] += MomResCentBin-1;
3423 MomRes3d[ChProd][term] = (TH3D*)MomResFile->Get(names3[ChProd][term]->Data());
3424 MomRes1d[ChProd][term] = (TH1D*)MomResFile->Get(names1D[ChProd][term]->Data());
3425 MomRes3d[ChProd][term]->SetDirectory(0);
3426 MomRes1d[ChProd][term]->SetDirectory(0);
3433 if(rvalue<=5) LambdaRbin = 1 + int(float(lambda-0.5+0.0001)/0.02);
3434 else LambdaRbin = 1 + 16*(rvalue-5) + int(float(lambda-0.5+0.0001)/0.02);
3435 Int_t LambdaRbinDev;
3436 if(rvalue<=5) LambdaRbinDev = 1 + 16*(1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3437 else LambdaRbinDev = 1 + 16*(rvalue-5 - 1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3439 for(int i=0; i<40; i++){
3440 MomRes_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbin, i+1);
3441 MomRes_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbin, i+1);
3442 MomRes_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbin, i+1);
3443 MomRes_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbin, i+1);
3444 MomRes_term1_mp[i] = MomResWeights_mp_term1->GetBinContent(LambdaRbin, i+1);
3445 MomRes_term2_mp[i] = MomResWeights_mp_term2->GetBinContent(LambdaRbin, i+1);
3447 MomResDev_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbinDev, i+1);
3448 MomResDev_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbinDev, i+1);
3449 MomResDev_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbinDev, i+1);
3450 MomResDev_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbinDev, i+1);
3453 MomResFile->Close();
3457 void fcn_r3_3d(int& npar, double* deriv, double& f, double par[], int flag){
3462 for(int i=fitlimitLowBin-1; i<fitlimitHighBin; i++){
3463 double q12=(i+0.5)*0.005;
3464 for(int j=fitlimitLowBin-1; j<fitlimitHighBin; j++){
3465 double q13=(j+0.5)*0.005;
3466 for(int k=fitlimitLowBin-1; k<fitlimitHighBin; k++){
3467 double q23=(k+0.5)*0.005;
3469 //r3 = 1 - par[2]*pow(q12*q13 + q12*q23 + q13*q23,2) + 1/5.*pow(par[0],2)*(3-par[1]*(q12*q12+q13*q13+q23*q23)) - 112/70.*pow(par[0],3);
3470 //r3 /= pow( (1-par[1]*q12*q12 - 4/5.*pow(par[0],2))*(1-par[1]*q13*q13 - 4/5.*pow(par[0],2))*(1-par[1]*q23*q23 - 4/5.*pow(par[0],2)),0.5);
3471 r3 = par[0] - par[1]*(q12*q13 + q12*q23 + q13*q23) - par[2]*pow(q12*q13 + q12*q23 + q13*q23,2);
3472 if(r3_3d_arr_e[i][j][k] > 0) Chi2sum += pow((r3_3d_arr[i][j][k] - r3)/r3_3d_arr_e[i][j][k],2);
3481 double r3_fitfunction(Double_t *x, Double_t *par){
3483 //double r3 = 1 - par[2]*pow(x[0]*x[1] + x[0]*x[2] + x[1]*x[2],2) + 1/5.*pow(par[0],2)*(3-par[1]*(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])) - 112/70.*pow(par[0],3);
3484 //r3 /= pow( (1-par[1]*x[0]*x[0] - 4/5.*pow(par[0],2))*(1-par[1]*x[1]*x[1] - 4/5.*pow(par[0],2))*(1-par[1]*x[2]*x[2] - 4/5.*pow(par[0],2)),0.5);
3485 double r3 = par[0] - par[1]*(x[0]*x[1] + x[0]*x[2] + x[1]*x[2]) - par[2]*pow(x[0]*x[1] + x[0]*x[2] + x[1]*x[2],2);
3491 double cubicInterpolate (double p[4], double x) {
3492 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
3494 double bicubicInterpolate (double p[4][4], double x, double y) {
3496 arr[0] = cubicInterpolate(p[0], y);
3497 arr[1] = cubicInterpolate(p[1], y);
3498 arr[2] = cubicInterpolate(p[2], y);
3499 arr[3] = cubicInterpolate(p[3], y);
3500 return cubicInterpolate(arr, x);
3502 double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
3504 arr[0] = bicubicInterpolate(p[0], y, z);
3505 arr[1] = bicubicInterpolate(p[1], y, z);
3506 arr[2] = bicubicInterpolate(p[2], y, z);
3507 arr[3] = bicubicInterpolate(p[3], y, z);
3508 return cubicInterpolate(arr, x);