]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Split: fix refs to AddTaskCentrality.C
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / Plot_FourPion.C
CommitLineData
264b1bb9 1#include <math.h>
2#include <time.h>
3#include <stdio.h>
4#include <stdlib.h>
5#include <Riostream.h>
6#include <assert.h>
7
8#include "TVector2.h"
9#include "TFile.h"
10#include "TString.h"
11#include "TF1.h"
12#include "TF3.h"
13#include "TH1.h"
14#include "TH2.h"
15#include "TH3.h"
16#include "TProfile.h"
17#include "TProfile2D.h"
18#include "TMath.h"
19#include "TText.h"
20#include "TRandom3.h"
21#include "TArray.h"
22#include "TLegend.h"
23#include "TStyle.h"
24#include "TMinuit.h"
25#include "TCanvas.h"
26#include "TLatex.h"
27#include "TPaveStats.h"
28#include "TASImage.h"
29
30#define BohrR 1963.6885 // Bohr Radius for pions
31#define FmToGeV 0.19733 // conversion to fm
32#define PI 3.1415926
33#define masspiC 0.1395702 // pi+ mass (GeV/c^2)
34
35using namespace std;
36// 0(standard), 1(fcSq=0.65), 2(fcSq=0.75), 3(B minus), 4(B plus),
37// 5(Linear instead of Cubic,also fc^2=0.75), 6(10h), 7(MRC 10% increase), 8(MuonCorr 92% pion-pair purity)
38int FileSetting=0;
39//
cbf4f1cb 40bool MCcase_def=0;// MC data?
57e5b35f 41int CollisionType=0;// PbPb, pPb, pp
264b1bb9 42//
43int Mbin_def=0;// 0-9: centrality bin in widths of 5%
cbf4f1cb 44bool SameCharge_def=1;// same-charge?
57e5b35f 45int CHARGE_def=1;// -1 or +1: + or - pions for same-charge case, --+ or -++, ---+ or -+++
cbf4f1cb 46int MixedCharge4pionType_def = 2;// 1(---+) or 2(--++)
264b1bb9 47//
48int EDbin_def=0;// 0 or 1: Kt3 bin
cbf4f1cb 49int TPNbin=0;// TPN bin for r3 and r4
57e5b35f 50int Gbin = int( (0) /2. ) + 5;// +5 (Rcoh=0), +55 (Rcoh=Rch)
51int c3FitType = 1;// EW(1), LG(2)
264b1bb9 52int Ktbin_def=1;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
53//
54bool MRC=1;// Momentum Resolution Corrections?
55bool MuonCorrection=1;// correct for Muons?
cbf4f1cb 56bool FSICorrection=1;// correct For Final-State-Interactions
57bool InterpCorrection=1;// correct for finite bin interpolation
264b1bb9 58//
59int f_choice=0;// 0(Core/Halo), 1(40fm), 2(70fm), 3(100fm)
60//
61//
62bool SaveToFile_def=kFALSE;// Save outputs to file?
63bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
64//
65//
66//
67//
68
264b1bb9 69int fFSIindex=0;
70int Ktlowbin;
71int Kthighbin;
72float TwoFrac;// Lambda parameter
73float OneFrac;// Lambda^{1/2}
74float ThreeFrac;// Lambda^{3/2}
75float FourFrac;// lambda^{4/2}
76double Qcut_pp = 0.6;// 0.6
77double Qcut_PbPb = 0.1;
78double NormQcutLow_pp = 1.0;
79double NormQcutHigh_pp = 1.5;
80double NormQcutLow_PbPb = .15;
81double NormQcutHigh_PbPb = .2;// was .175
82
cbf4f1cb 83int ThreeParticleRebin;
84int FourParticleRebin;
264b1bb9 85
86int RbinMRC;
87
88TH1D *fFSIss[12];
89TH1D *fFSIos[12];
90
91void SetFSICorrelations();
92void SetFSIindex(Float_t);
93Float_t FSICorrelation(Int_t, Int_t, Float_t);
94void SetMuonCorrections();
95void SetMomResCorrections();
96double Gamov(int, double);
97double cubicInterpolate(double[4], double);
98double bicubicInterpolate(double[4][4], double, double);
99double tricubicInterpolate(double[4][4][4], double, double, double);
100//
101float fact(float);
102
103
104TH1D *MRC_SC_2[2];
105TH1D *MRC_MC_2[2];
106TH1D *MRC_SC_3[3];
107TH1D *MRC_MC_3[4];
108TH1D *MRC_SC_4[5];
109TH1D *MRC_MC1_4[6];
110TH1D *MRC_MC2_4[6];
111
112
113double AvgQinvSS[30];
114double AvgQinvOS[30];
115double BinCentersQ4[20];
116
117//
118TH1D *C2muonCorrectionSC;
119TH1D *C2muonCorrectionMC;
120TH1D *WeightmuonCorrection;
121TH1D *C3muonCorrectionSC[2];
122TH1D *C3muonCorrectionMC[3];
123TH1D *C4muonCorrectionSC[4];
124TH1D *C4muonCorrectionMC1[5];
125TH1D *C4muonCorrectionMC2[5];
126
127
128void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool SameCharge=SameCharge_def, int MixedCharge4pionType=MixedCharge4pionType_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
129
130 EDbin_def=EDbin;
131 SaveToFile_def=SaveToFile;
132 MCcase_def=MCcase;
133 CHARGE_def=CHARGE;
134 SameCharge_def=SameCharge;// 3-pion same-charge
135 MixedCharge4pionType_def=MixedCharge4pionType;
136 Mbin_def=Mbin;
137 Ktbin_def=Ktbin;
138 //
139 Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
140 Kthighbin=(Ktbin)*2+4;
141
142 //
143 if(FileSetting==1) TwoFrac=0.65;
144 else if(FileSetting==2 || FileSetting==5) TwoFrac=0.75;
145 else TwoFrac=0.7;
146
147 OneFrac = sqrt(TwoFrac);
148 ThreeFrac = pow(TwoFrac, 1.5);
149 FourFrac = pow(TwoFrac, 2.);
150
cbf4f1cb 151 //// Core/Halo, 40fm, 70fm, 100fm
152 float ThermShift_f33[4]={0., 0.06933, 0.01637, 0.006326};
153 float ThermShift_f32[4]={0., -0.0185, -0.005301, -0.002286};
154 float ThermShift_f31[4]={0., -0.01382, -0.0004682, 0.0005337};
155 float f33prime = ThreeFrac;
156 float f32prime = TwoFrac*(1-OneFrac);
157 float f31prime = pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2);
158 f33prime += ThermShift_f33[f_choice];
159 f32prime += ThermShift_f32[f_choice];
160 f31prime += ThermShift_f31[f_choice];
161 float f33 = f33prime;
162 float f32 = f32prime/TwoFrac;
163 float f31 = f31prime - 3*((1-TwoFrac)*(1-OneFrac) + ThermShift_f32[f_choice]*(1-TwoFrac)/TwoFrac);
164 //cout<<f33 + 3*f32 + f31<<endl;
165
166 //// Core/Halo, 40fm, 70fm, 100fm
167 float ThermShift_f44[4]={0., 0.08741, 0.005284, -0.01064};
168 float ThermShift_f43[4]={0., -0.01126, -0.001486, 0.001686};
169 float ThermShift_f42[4]={0., -0.006466, -7.683e-05, 0.0004572};
170 float ThermShift_f41[4]={0., -0.003556, 0.00112, 0.00115};
171 float f44prime = FourFrac;
172 float f43prime = ThreeFrac*(1-OneFrac);
173 float f42prime = TwoFrac*pow(1-OneFrac,2);
174 float f41prime = pow(1-OneFrac,4) + 4*OneFrac*pow(1-OneFrac,3);
175 f44prime += ThermShift_f44[f_choice];
176 f43prime += ThermShift_f43[f_choice];
177 f42prime += ThermShift_f42[f_choice];
178 f41prime += ThermShift_f41[f_choice];
179 float f44 = f44prime;
180 float f43 = f43prime/f33prime;
181 float f42 = f42prime/TwoFrac;
182 f42 -= 2*f43prime*f32prime/f33prime/TwoFrac;
183 float f41 = f41prime;
184 f41 -= 4*f43prime*f31prime/f33prime;
185 f41 -= 6*f42prime*(1-TwoFrac)/TwoFrac;
186 f41 += 12*f43prime/f33prime*f32prime/TwoFrac*(1-TwoFrac);
187 //cout<<f44 + 4*f43 + 6*f42 + f41<<endl;
188 //cout<<f44<<" "<<f43<<" "<<f42<<" "<<f41<<endl;
189 //
190 // Core/Halo, 40fm, 70fm, 100fm
191 //float TherminatorMod_f33[4]={1., 1.123, 1.022, 1.008};// 1.008
192 //float TherminatorMod_f32[4]={1., 0.844, 0.933, 0.967};// .967
193 //float TherminatorMod_f31[4]={1., 0.828, 0.982, 1.028};// 1.028
194 /*float TherminatorMod_f44[4]={1., 1.188, 1.008, 0.985};// .985
195 float TherminatorMod_f43[4]={1., 0.885, 0.979, 1.026};// 1.026
196 float TherminatorMod_f42[4]={1., 0.687, 0.99, 1.08};// 1.08
197 float TherminatorMod_f41[4]={1., 0.806, 1.33, 1.548};//1.548
198 float f44 = FourFrac;
199 float f43 = (1-OneFrac);
200 float f42 = -pow(1-OneFrac,2);
201 float f41 = -3*pow(1-OneFrac,4) - 8*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2);
202 f44 *= TherminatorMod_f44[f_choice];
203 f43 *= TherminatorMod_f43[f_choice];
204 f42 *= TherminatorMod_f42[f_choice];
205 f41 *= TherminatorMod_f41[f_choice];*/
206 //f44prime *= TherminatorMod_f44[f_choice];
207 //f43prime *= TherminatorMod_f43[f_choice];
208 //f42prime *= TherminatorMod_f42[f_choice];
209 //f41prime *= TherminatorMod_f41[f_choice];
210 //float f44 = f44prime;
211 //float f43 = f43prime/ThreeFrac;
212 //float f42 = f42prime/TwoFrac - 2*pow(1-OneFrac,2);
213 //float f41 = f41prime - 4*pow(1-OneFrac,4) - 12*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2);
214
215 cout<<"Mbin = "<<Mbin<<" KT3 = "<<EDbin<<" SameCharge = "<<SameCharge<<endl;
264b1bb9 216 if(!SameCharge) cout<<"4-pion MixedCharge type = "<<MixedCharge4pionType<<endl;
217
cbf4f1cb 218 if(CollisionType==0){
219 if(Mbin==0) {RbinMRC=10;}
220 else if(Mbin==1) {RbinMRC=9;}
221 else if(Mbin<=3) {RbinMRC=8;}
222 else if(Mbin<=5) {RbinMRC=7;}
223 else {RbinMRC=6;}
224 }else{
225 RbinMRC=2;
226 }
227
228 if(MCcase) MuonCorrection=kFALSE;
264b1bb9 229
cbf4f1cb 230 if(CollisionType==0){
231 ThreeParticleRebin=2;
232 FourParticleRebin=3;
233 }else{
234 ThreeParticleRebin=2;
235 FourParticleRebin=6;
236 }
264b1bb9 237 float Cutoff_FullWeight_Q3[10]={0.065, 0.065, 0.075, 0.075, 0.075, 0.075, 0.085, 0.085, 0.085, 0.085};
238 float Cutoff_FullWeight_Q4[10]={0.115, 0.115, 0.115, 0.130, 0.130, 0.130, 0.145, 0.145, 0.145, 0.145};
239
240
241 //
242 // avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values. From unsmeared HIJING 0-5% with input signal
243 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};
244 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};
245 for(int i=0; i<30; i++){
246 AvgQinvSS[i] = AvgQinvSS_temp[i];
247 AvgQinvOS[i] = AvgQinvOS_temp[i];
248 }
249
250
251
252 TH2D *TwoParticle_2d[2][2][2];// ch1,ch2,term
253 TH1D *TwoParticle[2][2][2];// ch1,ch2,term
57e5b35f 254 TH2D *UnitMult_2d[2][2][2];// ch1,ch2,term
255 TH1D *UnitMult[2][2][2];// ch1,ch2,term
264b1bb9 256 double norm_2[2]={0};
57e5b35f 257 double norm_2_UM[2]={0};
264b1bb9 258 //
259 TH1D *ThreeParticle[2][2][2][5];// ch1,ch2,ch3,term
260 TProfile *K3avg[2][2][2][4];
261 TH2D *TPFullWeight_ThreeParticle_2D[2];// charge
262 TH2D *TPNegFullWeight_ThreeParticle_2D[2];// charge
263 TH1D *TPN_ThreeParticle[2];// charge
264 TH1D *TPFullWeight_ThreeParticle[2];// charge
265 TH1D *TPNegFullWeight_ThreeParticle[2];// charge
266 double norm_3[5]={0};
267 //
268 TH1D *FourParticle[2][2][2][2][13];// ch1,ch2,ch3,ch4,term
269 TProfile *K4avg[2][2][2][2][12];
270 TH2D *TPFullWeight_FourParticle_2D[2];// charge
271 TH2D *TPNegFullWeight_FourParticle_2D[2];// charge
272 TH1D *TPFullWeight_FourParticle[2];// charge
273 TH1D *TPNegFullWeight_FourParticle[2];// charge
274 TH1D *TPN_FourParticle[2];// charge
cbf4f1cb 275 TH3D *FullBuildFromFits_3D;// charge
276 TH1D *FullBuildFromFits;// charge
277 TH3D *PartialBuildFromFits_3D;// charge
278 TH1D *PartialBuildFromFits;// charge
264b1bb9 279 double norm_4[13]={0};
280
281
282 gStyle->SetOptStat(0);
283 gStyle->SetOptDate(0);
284 gStyle->SetOptFit(0111);
285
286 TFile *_file0;
cbf4f1cb 287 if(CollisionType==0){// PbPb
264b1bb9 288 if(MCcase) {
289 if(Mbin<=1){
cbf4f1cb 290 //_file0 = new TFile("Results/PDC_12a17a_GeneratedSignal.root","READ");
291 //_file0 = new TFile("Results/PDC_12a17a_pTSpectrumWeight.root","READ");
292 _file0 = new TFile("Results/PDC_12a17a_Qweights.root","READ");
264b1bb9 293 }
294 }else{
cbf4f1cb 295 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0_FullRunWrongWeightsNoPadRowTTCandMCTTC_RawWeightFile.root","READ");
296 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p16to0p25_0p25to1p0.root","READ");
297 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0.root","READ");
298 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT0p16to0p25_KT0p35.root","READ");
299 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_1percentCentEM.root","READ");
300 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_0p02eta0p045phi_0p03eta0p067phi.root","READ");
57e5b35f 301 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_c3FitBuild.root","READ");
302 if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_LowNorm_HighNorm.root","READ");
cbf4f1cb 303 //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_extendedGweights.root","READ");// Preliminary results
264b1bb9 304 if(FileSetting==5) _file0 = new TFile("Results/PDC_11h_Cubic_Linear.root","READ");
305 if(FileSetting==1 || FileSetting==2) _file0 = new TFile("Results/PDC_11h_Lam0p65_Lam0p75.root","READ");
306 if(FileSetting==3) _file0 = new TFile("Results/PDC_11h_Cubic_Linear_Bminus.root","READ");
307 if(FileSetting==4) _file0 = new TFile("Results/PDC_11h_Cubic_Linear_Bplus.root","READ");
308 if(FileSetting==6) _file0 = new TFile("Results/PDC_10h_Cubic_Linear.root","READ");
309 if(FileSetting==7 || FileSetting==8) _file0 = new TFile("Results/PDC_11h_MRC10percIncrease_Muon92percent.root","READ");
310 }
cbf4f1cb 311 }else if(CollisionType==1){// pPb
57e5b35f 312 _file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm_HighNorm.root","READ");
264b1bb9 313 }else{// pp
57e5b35f 314 _file0 = new TFile("Results/PDC_10bcde_kMB_LowNorm_HighNorm.root","READ");
264b1bb9 315 }
316
317 SetFSIindex(10.);
318 SetFSICorrelations();
319 SetMomResCorrections();
320 SetMuonCorrections();
321 //
322 /////////////////////////////////////////////////////////
323
324
325 double NormQcutLow;
326 double NormQcutHigh;
cbf4f1cb 327 if(CollisionType==0) {
264b1bb9 328 NormQcutLow = NormQcutLow_PbPb;
329 NormQcutHigh = NormQcutHigh_PbPb;
330 }else {
331 NormQcutLow = NormQcutLow_pp;
332 NormQcutHigh = NormQcutHigh_pp;
333 }
334
335
336 TList *MyList;
337 if(!MCcase){
338 TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root");
57e5b35f 339 if(FileSetting==2 || FileSetting==5 || FileSetting==8 || CollisionType!=0) MyList=(TList*)tdir->Get("FourPionOutput_2");
264b1bb9 340 else MyList=(TList*)tdir->Get("FourPionOutput_1");
341 //MyList=(TList*)_file0->Get("MyList");
342 }else{
343 MyList=(TList*)_file0->Get("MyList");
344 }
345 _file0->Close();
346
347
348 TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
349 //
350
351 cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
352
353
354
355 ///////////////////////////////////
356 // Get Histograms
357
358 // 2-particle
359 for(int c1=0; c1<2; c1++){
360 for(int c2=0; c2<2; c2++){
361 for(int term=0; term<2; term++){
362
363 TString *name2 = new TString("TwoParticle_Charge1_");
364 *name2 += c1;
365 name2->Append("_Charge2_");
366 *name2 += c2;
367 name2->Append("_M_");
368 *name2 += Mbin;
369 name2->Append("_ED_");
370 *name2 += 0;// EDbin
371 name2->Append("_Term_");
372 *name2 += term+1;
373
374 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
375
376 TwoParticle_2d[c1][c2][term] = (TH2D*)MyList->FindObject(name2->Data());
377 TwoParticle_2d[c1][c2][term]->Sumw2();
378 TString *proname = new TString(name2->Data());
379 proname->Append("_pro");
380
381 if(Ktbin==10) {Ktlowbin=1; Kthighbin=TwoParticle_2d[c1][c2][term]->GetNbinsX();}
382 TwoParticle[c1][c2][term] = (TH1D*)TwoParticle_2d[c1][c2][term]->ProjectionY(proname->Data(),Ktlowbin,Kthighbin);// bins:5-6 (kt:0.2-0.3)
383
384 norm_2[term] = TwoParticle[c1][c2][term]->Integral(TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutLow),TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutHigh));
cbf4f1cb 385 //cout<<"2-pion norms "<<norm_2[term]<<endl;
264b1bb9 386 TwoParticle[c1][c2][term]->Scale(norm_2[0]/norm_2[term]);
387
388 TwoParticle[c1][c2][term]->SetMarkerStyle(20);
389 TwoParticle[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
390 TwoParticle[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
391 TwoParticle[c1][c2][term]->SetTitle("");
57e5b35f 392 //
264b1bb9 393
57e5b35f 394 TString *nameUM=new TString(name2->Data());
395 nameUM->Append("_UnitMult");
396 UnitMult_2d[c1][c2][term] = (TH2D*)MyList->FindObject(nameUM->Data());
397 UnitMult_2d[c1][c2][term]->Sumw2();
398 TString *pronameUM = new TString(nameUM->Data());
399 pronameUM->Append("_pro");
400 UnitMult[c1][c2][term] = (TH1D*)UnitMult_2d[c1][c2][term]->ProjectionY(pronameUM->Data(),13,13);// 11 means 1000 pions
401 norm_2_UM[term] = UnitMult[c1][c2][term]->Integral(UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.0),UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.2));
402 UnitMult[c1][c2][term]->Scale(norm_2_UM[0]/norm_2_UM[term]);
403 //
404 UnitMult[c1][c2][term]->SetMarkerStyle(20);
405 UnitMult[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
406 UnitMult[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
407 UnitMult[c1][c2][term]->SetTitle("");
264b1bb9 408 }// term
409
410
411 // 3-particle
412 for(int c3=0; c3<2; c3++){
413 for(int term=0; term<5; term++){
414
415 TString *name3 = new TString("ThreeParticle_Charge1_");
416 *name3 += c1;
417 name3->Append("_Charge2_");
418 *name3 += c2;
419 name3->Append("_Charge3_");
420 *name3 += c3;
421 name3->Append("_M_");
422 *name3 += Mbin;
423 name3->Append("_ED_");
424 *name3 += EDbin;
425 name3->Append("_Term_");
426 *name3 += term+1;
427
428 TString *nameNorm3=new TString(name3->Data());
429 nameNorm3->Append("_Norm");
430 //
431 TString *nameK3=new TString(name3->Data());
432 nameK3->Append("_Kfactor");
433 //
434 TString *nameTPN3=new TString(name3->Data());
435 nameTPN3->Append("_TwoPartNorm");
436 //
437 TString *nameNegTPN3=new TString(name3->Data());
438 nameNegTPN3->Append("_TwoPartNegNorm");
439 //
440 name3->Append("_1D");
441
442 ///////////////////////////////////////
443 // skip degenerate histograms
444 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
445 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
446 ////////////////////////////////////////
447
448 norm_3[term] = ((TH1D*)MyList->FindObject(nameNorm3->Data()))->Integral();
449 ThreeParticle[c1][c2][c3][term] = (TH1D*)MyList->FindObject(name3->Data());
450 ThreeParticle[c1][c2][c3][term]->Sumw2();
451 //if(c1==c2 && c1==c3) cout<<"3-pion norms "<<norm_3[term]<<endl;
452 ThreeParticle[c1][c2][c3][term]->Scale(norm_3[0]/norm_3[term]);
453 ThreeParticle[c1][c2][c3][term]->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
454 ThreeParticle[c1][c2][c3][term]->GetYaxis()->SetTitle("C_{3}");
455 ThreeParticle[c1][c2][c3][term]->SetMarkerStyle(20);
456 ThreeParticle[c1][c2][c3][term]->SetTitle("");
457 //
458 ThreeParticle[c1][c2][c3][term]->Rebin(ThreeParticleRebin);
459
460 //
461 if(term<4){
462 K3avg[c1][c2][c3][term] = (TProfile*)MyList->FindObject(nameK3->Data());
463 K3avg[c1][c2][c3][term]->Rebin(ThreeParticleRebin);
cbf4f1cb 464 if(MCcase || !FSICorrection) {
465 K3avg[c1][c2][c3][term]->Reset();
466 for(int ii=1; ii<=K3avg[c1][c2][c3][term]->GetNbinsX(); ii++) {
467 K3avg[c1][c2][c3][term]->Fill(K3avg[c1][c2][c3][term]->GetXaxis()->GetBinCenter(ii), 1);
468 }
469 }
264b1bb9 470 }
471 //
472 if(term==4 && c1==c2 && c1==c3){
473 TPFullWeight_ThreeParticle_2D[c1] = (TH2D*)MyList->FindObject(nameTPN3->Data());
474 TPFullWeight_ThreeParticle_2D[c1]->Scale(norm_3[0]/norm_3[term]);
475 TPFullWeight_ThreeParticle_2D[c1]->RebinY(ThreeParticleRebin);
476 //
477 TPNegFullWeight_ThreeParticle_2D[c1] = (TH2D*)MyList->FindObject(nameNegTPN3->Data());
478 TPNegFullWeight_ThreeParticle_2D[c1]->Scale(norm_3[0]/norm_3[term]);
479 TPNegFullWeight_ThreeParticle_2D[c1]->RebinY(ThreeParticleRebin);
480 //
481 for(int binx=1; binx<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsX(); binx++) {
482 for(int biny=1; biny<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsY(); biny++) {
483 TPFullWeight_ThreeParticle_2D[c1]->SetBinError(binx, biny, 0);
cbf4f1cb 484 if(binx!=4){
485 TPFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny));
486 TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny));
487 if(InterpCorrection){
488 double q3 = TPFullWeight_ThreeParticle_2D[c1]->GetYaxis()->GetBinCenter(biny);
489 double InterCorr = 1.024 - 0.2765*q3;
490 //if(q3<0.1) cout<<binx<<" "<<biny<<" "<<q3<<" "<<InterCorr<<endl;
491 TPFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny));
492 TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny));
493 }
494 }
264b1bb9 495 }
496 }
cbf4f1cb 497
264b1bb9 498 TString *proName=new TString(nameTPN3->Data()); TString *proNameNeg=new TString(nameNegTPN3->Data());
499 proName->Append("_pro"); proNameNeg->Append("_pro");
500 TPN_ThreeParticle[c1] = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin);
501 //
502 proName->Append("_FullWeight"); proNameNeg->Append("_FullWeight");
503 TPFullWeight_ThreeParticle[c1] = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), Gbin, Gbin);
504 TPNegFullWeight_ThreeParticle[c1] = (TH1D*)TPNegFullWeight_ThreeParticle_2D[c1]->ProjectionY(proNameNeg->Data(), Gbin, Gbin);
505 proName->Append("_FullWeightDen"); proNameNeg->Append("_FullWeightDen");
506 TH1D *tempDen = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), 4, 4);
507 TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[c1]->ProjectionY(proNameNeg->Data(), 4, 4);
508 // Add Pos with Neg weights
509 tempDen->Add(tempDenNeg);
510 TPFullWeight_ThreeParticle[c1]->Add(TPNegFullWeight_ThreeParticle[c1]);
511 //
cbf4f1cb 512 //TPFullWeight_ThreeParticle[c1]->Add(tempDen);// now added above in Interp section
264b1bb9 513 TPFullWeight_ThreeParticle[c1]->Divide(tempDen);
514 TPFullWeight_ThreeParticle[c1]->SetLineColor(1);
515 //
516 }
517
518 }// term 3-particle
519
520
521 // 4-particle
522 for(int c4=0; c4<2; c4++){
523 for(int term=0; term<13; term++){
524
525 TString *name4 = new TString("FourParticle_Charge1_");
526 *name4 += c1;
527 name4->Append("_Charge2_");
528 *name4 += c2;
529 name4->Append("_Charge3_");
530 *name4 += c3;
531 name4->Append("_Charge4_");
532 *name4 += c4;
533 name4->Append("_M_");
534 *name4 += Mbin;
535 name4->Append("_ED_");
536 *name4 += EDbin;
537 name4->Append("_Term_");
538 *name4 += term+1;
539
540 TString *nameNorm4=new TString(name4->Data());
541 nameNorm4->Append("_Norm");
542 //
543 TString *nameK4=new TString(name4->Data());
cbf4f1cb 544 nameK4->Append("_Kfactor");// or Weighted
264b1bb9 545 //
546 TString *nameTPN4=new TString(name4->Data());
547 nameTPN4->Append("_TwoPartNorm");
548 //
549 TString *nameNegTPN4=new TString(name4->Data());
550 nameNegTPN4->Append("_TwoPartNegNorm");
551 //
cbf4f1cb 552 TString *nameFitBuild=new TString(name4->Data());
553 nameFitBuild->Append("_FullBuildFromFits");
554 //
555 TString *namePartialFitBuild=new TString(name4->Data());
556 namePartialFitBuild->Append("_PartialBuildFromFits");
557 //
264b1bb9 558 name4->Append("_1D");
559 ///////////////////////////////////////
560 // skip degenerate histograms
561 if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
562 if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
563 if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
564 /////////////////////////////////////////
565 norm_4[term] = ((TH1D*)MyList->FindObject(nameNorm4->Data()))->Integral();
57e5b35f 566
264b1bb9 567 //if( (c1+c2+c3+c4)==4) cout<<"4-pion norms "<<norm_4[term]<<endl;
568 if(norm_4[term] > 0){
57e5b35f 569 //if(c1==c2 && c1==c3 && c1==c4) cout<<term<<" "<<norm_4[0]/norm_4[term]<<endl;
264b1bb9 570 FourParticle[c1][c2][c3][c4][term] = (TH1D*)MyList->FindObject(name4->Data());
571 FourParticle[c1][c2][c3][c4][term]->Sumw2();
572 FourParticle[c1][c2][c3][c4][term]->Scale(norm_4[0]/norm_4[term]);
573 FourParticle[c1][c2][c3][c4][term]->GetXaxis()->SetTitle("Q_{4} (GeV/c)");
574 FourParticle[c1][c2][c3][c4][term]->GetYaxis()->SetTitle("C_{4}");
575 FourParticle[c1][c2][c3][c4][term]->SetMarkerStyle(20);
576 FourParticle[c1][c2][c3][c4][term]->SetTitle("");
577 //
578 FourParticle[c1][c2][c3][c4][term]->Rebin(FourParticleRebin);
cbf4f1cb 579
264b1bb9 580 }else{
cbf4f1cb 581 cout<<"4-particle normalization = 0."<<endl;
264b1bb9 582 }
583 if(term<12){
584 K4avg[c1][c2][c3][c4][term] = (TProfile*)MyList->FindObject(nameK4->Data());
585 K4avg[c1][c2][c3][c4][term]->Rebin(FourParticleRebin);
cbf4f1cb 586 if(MCcase || !FSICorrection) {
587 K4avg[c1][c2][c3][c4][term]->Reset();
588 for(int ii=1; ii<=K4avg[c1][c2][c3][c4][term]->GetNbinsX(); ii++) {
589 K4avg[c1][c2][c3][c4][term]->Fill(K4avg[c1][c2][c3][c4][term]->GetXaxis()->GetBinCenter(ii), 1);
590 }
591 }
264b1bb9 592 }
593 if(term==12 && c1==c2 && c1==c3 && c1==c4){
57e5b35f 594
264b1bb9 595 TPFullWeight_FourParticle_2D[c1] = (TH2D*)MyList->FindObject(nameTPN4->Data());
596 TPFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]);
597 TPFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin);
598 //
599 TPNegFullWeight_FourParticle_2D[c1] = (TH2D*)MyList->FindObject(nameNegTPN4->Data());
600 TPNegFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]);
601 TPNegFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin);
602 //
57e5b35f 603 if(c1==0 && !MCcase){
cbf4f1cb 604 FullBuildFromFits_3D = (TH3D*)MyList->FindObject(nameFitBuild->Data());
605 FullBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]);
606 FullBuildFromFits_3D->RebinZ(FourParticleRebin);
607 //
608 PartialBuildFromFits_3D = (TH3D*)MyList->FindObject(namePartialFitBuild->Data());
609 PartialBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]);
610 PartialBuildFromFits_3D->RebinZ(FourParticleRebin);
611 }
612 //
264b1bb9 613 for(int binx=1; binx<=TPFullWeight_FourParticle_2D[c1]->GetNbinsX(); binx++) {
614 for(int biny=1; biny<=TPFullWeight_FourParticle_2D[c1]->GetNbinsY(); biny++) {
615 TPFullWeight_FourParticle_2D[c1]->SetBinError(binx, biny, 0);
cbf4f1cb 616 if(binx!=4){
617 TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny));
618 TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny));
619
620 if(InterpCorrection){// Interpolator correction
621 double q4 = TPFullWeight_FourParticle_2D[c1]->GetYaxis()->GetBinCenter(biny);
622 double InterCorr = 1.032 - 0.2452*q4;
623 TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny));
624 TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny));
625 }
626 }
264b1bb9 627 }
628 }
629 TString *proName=new TString(nameTPN4->Data()); TString *proNegName=new TString(nameNegTPN4->Data());
cbf4f1cb 630 TString *proNameFitBuild=new TString(nameFitBuild->Data()); TString *proNamePartialFitBuild=new TString(namePartialFitBuild->Data());
631 proName->Append("_pro"); proNegName->Append("_pro"); proNameFitBuild->Append("_pro"); proNamePartialFitBuild->Append("_pro");
264b1bb9 632 TPN_FourParticle[c1] = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin);
633 //
634 proName->Append("_FullWeight"); proNegName->Append("_FullWeight");
635 TPFullWeight_FourParticle[c1] = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), Gbin, Gbin);
636 TPNegFullWeight_FourParticle[c1] = (TH1D*)TPNegFullWeight_FourParticle_2D[c1]->ProjectionY(proNegName->Data(), Gbin, Gbin);
637 proName->Append("_FullWeightDen"); proNegName->Append("_FullWeightDen");
638 TH1D *tempDen = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), 4, 4);
639 TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[c1]->ProjectionY(proNegName->Data(), 4, 4);
640 //
641 // Add Pos and Neg Weights
642 tempDen->Add(tempDenNeg);
643 TPFullWeight_FourParticle[c1]->Add(TPNegFullWeight_FourParticle[c1]);
644 //
cbf4f1cb 645 //TPFullWeight_FourParticle[c1]->Add(tempDen);// now added above in Interp section
264b1bb9 646 TPFullWeight_FourParticle[c1]->Divide(tempDen);
647 TPFullWeight_FourParticle[c1]->SetLineColor(1);
648 /*TString *ErrName=new TString(nameTPN4->Data());
649 ErrName->Append("Err");
650 TH2D *temperr2D = (TH2D*)MyList->FindObject(ErrName->Data());
651 TH1D *temperr = (TH1D*)temperr2D->ProjectionY("tesst",4,4);
652 temperr->Rebin(FourParticleRebin);
653 cout.precision(8);
654 cout<<temperr->GetBinContent(3)<<endl;
655 cout<<(temperr->GetBinContent(5) / tempDen->GetBinContent(5))<<" "<<TPFullWeight_FourParticle[c1]->GetBinContent(5)<<endl;
656 */
57e5b35f 657 if(c1==0 && !MCcase){
cbf4f1cb 658 FullBuildFromFits = (TH1D*)FullBuildFromFits_3D->ProjectionZ(proNameFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin);
659 TH1D *tempDen2 = (TH1D*)FullBuildFromFits_3D->ProjectionZ("tempDen2", c3FitType, c3FitType, 4, 4);
660 tempDen2->Scale(1/100.);// It was filled 100 times with the same value
661 FullBuildFromFits->Add(tempDen2);
662 FullBuildFromFits->Divide(tempDen2);
57e5b35f 663 FullBuildFromFits->SetLineColor(4);
cbf4f1cb 664 //
665 PartialBuildFromFits = (TH1D*)PartialBuildFromFits_3D->ProjectionZ(proNamePartialFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin);
666 PartialBuildFromFits->Add(tempDen2);
667 PartialBuildFromFits->Divide(tempDen2);
668 PartialBuildFromFits->SetLineColor(2);
669 }
264b1bb9 670 }
671 }// term 4-particle
672
673 }// c4
674 }// c3
675 }// c2
676 }// c1
677
678 TH1D *fTPNRejects3pion = (TH1D*)MyList->FindObject("fTPNRejects3pion2");
679 TH1D *fTPNRejects4pion1 = (TH1D*)MyList->FindObject("fTPNRejects4pion1");
cbf4f1cb 680
681 //TH2D *c4QSFitNum2D = (TH2D*)MyList->FindObject("fc4QSFitNum");
682 //TH2D *c4QSFitDen2D = (TH2D*)MyList->FindObject("fc4QSFitDen");
683 //c4QSFitNum2D->RebinY(3);
684 //c4QSFitDen2D->RebinY(3);
685
264b1bb9 686 cout<<"Done getting Histograms"<<endl;
687
688 TF1 *Unity=new TF1("Unity","1",0,100);
689 Unity->SetLineStyle(2);
690
691
692 int ch1_2=0,ch2_2=0;
693 int ch1_3=0,ch2_3=0,ch3_3=0;
694 int ch1_4=0,ch2_4=0,ch3_4=0,ch4_4=0;
695
696 if(SameCharge && CHARGE==+1) {ch1_2=1; ch2_2=1; ch1_3=1; ch2_3=1; ch3_3=1; ch1_4=1; ch2_4=1; ch3_4=1; ch4_4=1;}
697 if(!SameCharge){
698 ch1_2=0; ch2_2=1;
699 //
700 if(CHARGE==-1) { ch1_3=0; ch2_3=0; ch3_3=1; }
701 else { ch1_3=0; ch2_3=1; ch3_3=1; }
702 //
703 if(CHARGE==-1 && MixedCharge4pionType==1) {ch1_4=0; ch2_4=0; ch3_4=0; ch4_4=1;}
704 if(CHARGE==+1 && MixedCharge4pionType==1) {ch1_4=0; ch2_4=1; ch3_4=1; ch4_4=1;}
705 if(MixedCharge4pionType==2) {ch1_4=0; ch2_4=0; ch3_4=1; ch4_4=1;}
706 }
707
708 ///////////////////////////////////////////////////////////
709 // 2-pion section
710 cout<<"2-pion section"<<endl;
711
712 TCanvas *can1 = new TCanvas("can1", "can1",11,53,700,500);
713 can1->SetHighLightColor(2);
714 can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
715 gStyle->SetOptFit(0111);
716 can1->SetFillColor(10);
717 can1->SetBorderMode(0);
718 can1->SetBorderSize(2);
719 can1->SetGridx();
720 can1->SetGridy();
721 can1->SetFrameFillColor(0);
722 can1->SetFrameBorderMode(0);
723 can1->SetFrameBorderMode(0);
724 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
725
726 TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
727 legend1->SetBorderSize(1);
728 legend1->SetTextSize(.04);// small .03; large .036
729 legend1->SetFillColor(0);
730
731
732 TH1D *TERM1_2=(TH1D*)TwoParticle[ch1_2][ch2_2][0]->Clone();
733 TH1D *TERM2_2=(TH1D*)TwoParticle[ch1_2][ch2_2][1]->Clone();
734 //
735 if(SameCharge){
736 TERM1_2->Multiply(MRC_SC_2[0]);
737 TERM2_2->Multiply(MRC_SC_2[1]);
738 }else {
739 TERM1_2->Multiply(MRC_MC_2[0]);
740 TERM2_2->Multiply(MRC_MC_2[1]);
741 }
cbf4f1cb 742
264b1bb9 743 TH1D *C2=(TH1D*)TERM1_2->Clone();
744 C2->Divide(TERM2_2);
745 C2->GetXaxis()->SetRangeUser(0,0.1);
cbf4f1cb 746 C2->SetMinimum(0.98); C2->SetMaximum(2.5);
747 /*C2->GetXaxis()->SetTitleSize(0.06);
748 C2->GetXaxis()->SetLabelSize(0.06);
749 C2->GetYaxis()->SetTitleSize(0.06);
750 C2->GetYaxis()->SetLabelSize(0.06);
751 C2->GetXaxis()->SetTitleOffset(1.05);
752 C2->GetYaxis()->SetTitleOffset(1.1);
753 C2->GetXaxis()->SetNdivisions(606);
754 C2->GetYaxis()->SetNdivisions(505);*/
264b1bb9 755 C2->Draw();
57e5b35f 756
757
264b1bb9 758 TH1D *C2QS=(TH1D*)TERM1_2->Clone();
759 C2QS->Add(TERM2_2, -(1-TwoFrac));
760 C2QS->Scale(1/TwoFrac);
761 for(int bin=1; bin<=C2QS->GetNbinsX(); bin++) {
cbf4f1cb 762 Float_t K2 = 1.0;
763 if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, C2QS->GetXaxis()->GetBinCenter(bin));
264b1bb9 764 C2QS->SetBinContent(bin, C2QS->GetBinContent(bin)/K2);
765 C2QS->SetBinError(bin, C2QS->GetBinError(bin)/K2);
766 }
767 C2QS->Divide(TERM2_2);
768 if(SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionSC);
769 if(!SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionMC);
cbf4f1cb 770 C2QS->SetMarkerColor(2); C2QS->SetLineColor(2);
771 if(!MCcase) C2QS->Draw("same");
772
57e5b35f 773
774 // To visualize the Qcut and Norm Regions
775 //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
776 //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);
777 //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);
778 //TERM1_2->Scale(1/TERM1_2->Integral());
779 //for(int bin=1; bin<=16; bin++) QcutRegion->SetBinContent(bin,TERM1_2->GetBinContent(bin));
780 //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
781 //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
782 //QcutRegion->SetFillColor(4);
783 //NormRegion1->SetFillColor(2);
784 //NormRegion2->SetFillColor(3);
785 //TERM1_2->GetYaxis()->SetTitle("Pair probability");
786 //TERM1_2->GetYaxis()->SetTitleOffset(1.3);
787 //TERM1_2->GetXaxis()->SetRangeUser(0,1.6);
788 //TERM1_2->Draw();
789 //QcutRegion->Draw("same");
790 //cout<<TERM1_2->Integral(1,16) / TERM1_2->Integral()<<endl;
791
cbf4f1cb 792 // Print 2-pion values
793 /*for(int bin=1; bin<=20; bin++){
794 cout<<C2QS->GetBinContent(bin)<<", ";
795 }
796 cout<<endl;
797 for(int bin=1; bin<=20; bin++){
798 cout<<C2QS->GetBinError(bin)<<", ";
799 }
800 cout<<endl;*/
801 double C2refPoints_0p65[20]={-0.218603, 1.00087, 1.00215, 1.00209, 1.00202, 1.00215, 1.00162, 1.00129, 1.00118, 1.00076, 1.00049, 1.00021, 0.999996, 0.999733, 0.999622, 0.999554, 0.999764, 0.999679, 0.999677, 0.999641};
802 double C2refPoints_e_0p65[20]={0.309152, 0.000473757, 0.000305599, 0.00022603, 0.000179596, 0.000149295, 0.000128009, 0.000112371, 0.00010047, 9.11363e-05, 8.37061e-05, 7.76871e-05, 7.27445e-05, 6.86217e-05, 6.5155e-05, 6.22148e-05, 5.9718e-05, 5.75613e-05, 5.57066e-05, 5.41028e-05};
803 double C2refPoints_0p75[20]={-0.135326, 0.968257, 0.984681, 0.991465, 0.99503, 0.997249, 0.998093, 0.998645, 0.999135, 0.999172, 0.999237, 0.999213, 0.999197, 0.999095, 0.999105, 0.999123, 0.999385, 0.999366, 0.999408, 0.999411};
804 double C2refPoints_e_0p75[20]={0.19138, 0.000422359, 0.000272193, 0.000201183, 0.000159788, 0.000132797, 0.000113845, 9.99271e-05, 8.93377e-05, 8.10334e-05, 7.44239e-05, 6.90702e-05, 6.46743e-05, 6.10077e-05, 5.79248e-05, 5.53103e-05, 5.30903e-05, 5.11725e-05, 4.95234e-05, 4.80974e-05};
805 TH1D *C2ref_0p65=(TH1D*)C2QS->Clone();
806 TH1D *C2ref_0p75=(TH1D*)C2QS->Clone();
807 for(int bin=1; bin<=20; bin++){
808 C2ref_0p65->SetBinContent(bin, C2refPoints_0p65[bin-1]);
809 C2ref_0p65->SetBinError(bin, C2refPoints_e_0p65[bin-1]);
810 C2ref_0p75->SetBinContent(bin, C2refPoints_0p75[bin-1]);
811 C2ref_0p75->SetBinError(bin, C2refPoints_e_0p75[bin-1]);
812 }
813 C2ref_0p65->SetMarkerColor(4); C2ref_0p65->SetLineColor(4);
814 C2ref_0p75->SetMarkerColor(6); C2ref_0p75->SetLineColor(6);
815 //C2ref_0p65->Draw("same");
816 //C2ref_0p75->Draw("same");
817
818
57e5b35f 819 TH1D *UnitMultC2=(TH1D*)UnitMult[0][0][0]->Clone();
820 TH1D *UnitMultC2den=(TH1D*)UnitMult[0][0][1]->Clone();
821 UnitMultC2->Divide(UnitMultC2den);
822 UnitMultC2->GetXaxis()->SetRangeUser(0,0.3);
823 UnitMultC2->SetMinimum(0.97);
824 UnitMultC2->SetMaximum(1.25);
825 UnitMultC2->Draw();
826
827 TH1D *UnitMultC2QS=(TH1D*)UnitMult[0][0][0]->Clone();
828 UnitMultC2QS->Add(UnitMult[0][0][1], -(1-TwoFrac));
829 UnitMultC2QS->Scale(1/TwoFrac);
830 for(int bin=1; bin<=UnitMultC2QS->GetNbinsX(); bin++) {
831 Float_t K2 = 1.0;
832 if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, UnitMultC2QS->GetXaxis()->GetBinCenter(bin));
833 UnitMultC2QS->SetBinContent(bin, UnitMultC2QS->GetBinContent(bin)/K2);
834 UnitMultC2QS->SetBinError(bin, UnitMultC2QS->GetBinError(bin)/K2);
835 }
836 UnitMultC2QS->Divide(UnitMultC2den);
837 UnitMultC2QS->SetMarkerColor(2); UnitMultC2QS->SetLineColor(2);
838 UnitMultC2QS->Draw("same");
264b1bb9 839
840 ///////////////////////////////////////////////////////////
841 // 3-pion
842 cout<<"3-pion section"<<endl;
843
844 TCanvas *can2 = new TCanvas("can2", "can2",600,53,700,500);
845 can2->SetHighLightColor(2);
846 can2->Range(-0.7838115,-1.033258,0.7838115,1.033258);
847 gStyle->SetOptFit(0111);
848 can2->SetFillColor(10);
849 can2->SetBorderMode(0);
850 can2->SetBorderSize(2);
851 can2->SetGridx();
852 can2->SetGridy();
853 can2->SetFrameFillColor(0);
854 can2->SetFrameBorderMode(0);
855 can2->SetFrameBorderMode(0);
856 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
857 TLegend *legend2 = (TLegend*)legend1->Clone();
858 legend2->SetX1(0.45); legend2->SetX2(0.95); legend2->SetY1(0.6); legend2->SetY2(0.95);
859
860 TH1D *TERM1_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][0]->Clone();
861 TH1D *TERM2_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][1]->Clone();
862 TH1D *TERM3_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][2]->Clone();
863 TH1D *TERM4_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][3]->Clone();
864 TH1D *TERM5_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][4]->Clone();
865 //
57e5b35f 866
264b1bb9 867 if(SameCharge){
868 TERM1_3->Multiply(MRC_SC_3[0]);
869 TERM2_3->Multiply(MRC_SC_3[1]);
870 TERM3_3->Multiply(MRC_SC_3[1]);
871 TERM4_3->Multiply(MRC_SC_3[1]);
872 TERM5_3->Multiply(MRC_SC_3[2]);
873 }else{
874 if(CHARGE==-1){// --+ but MRC is stored as -++
875 TERM1_3->Multiply(MRC_MC_3[0]);
876 TERM2_3->Multiply(MRC_MC_3[2]);
877 TERM3_3->Multiply(MRC_MC_3[1]);
878 TERM4_3->Multiply(MRC_MC_3[1]);
879 TERM5_3->Multiply(MRC_MC_3[3]);
880 }else{
881 TERM1_3->Multiply(MRC_MC_3[0]);
882 TERM2_3->Multiply(MRC_MC_3[1]);
883 TERM3_3->Multiply(MRC_MC_3[1]);
884 TERM4_3->Multiply(MRC_MC_3[2]);
885 TERM5_3->Multiply(MRC_MC_3[3]);
886 }
887 }
57e5b35f 888
264b1bb9 889
890 TH1D *N3QS = (TH1D*)TERM1_3->Clone();
891 N3QS->Add(TERM5_3, -f31);
892 N3QS->Add(TERM2_3, -f32);
893 N3QS->Add(TERM3_3, -f32);
894 N3QS->Add(TERM4_3, -f32);
895 N3QS->Scale(1/f33);
896 N3QS->Multiply(K3avg[ch1_3][ch2_3][ch3_3][0]);
897
898
899 TH1D *C3QS = (TH1D*)N3QS->Clone();
900 C3QS->Divide(TERM5_3);
901 if(SameCharge) C3QS->Multiply(C3muonCorrectionSC[0]);
902 else C3QS->Multiply(C3muonCorrectionMC[0]);
903
57e5b35f 904
264b1bb9 905 C3QS->SetMarkerStyle(20);
906 C3QS->SetMarkerColor(4);
907 C3QS->GetXaxis()->SetRangeUser(0,0.1);
908 C3QS->SetMinimum(0.9); C3QS->SetMaximum(5.0);
909
910 TH1D *n3QS = (TH1D*)N3QS->Clone();
911 TH1D *c3QS = (TH1D*)N3QS->Clone();
912 for(int bin=1; bin<=n3QS->GetNbinsX(); bin++){
913 if(n3QS->GetBinContent(bin)==0) continue;
914 double MuonCorr1=1.0, MuonCorr2=1.0, MuonCorr3=1.0, MuonCorr4=1.0;
915 if(SameCharge){
916 MuonCorr1 = C3muonCorrectionSC[0]->GetBinContent(bin);
917 MuonCorr2 = C3muonCorrectionSC[1]->GetBinContent(bin);
918 MuonCorr3 = MuonCorr2;
919 MuonCorr4 = MuonCorr2;
920 }else if(CHARGE==-1){
921 MuonCorr1 = C3muonCorrectionMC[0]->GetBinContent(bin);
922 MuonCorr2 = C3muonCorrectionMC[2]->GetBinContent(bin);// --
923 MuonCorr3 = C3muonCorrectionMC[1]->GetBinContent(bin);// -+
924 MuonCorr4 = MuonCorr3;// -+
925 }else{
926 MuonCorr1 = C3muonCorrectionMC[0]->GetBinContent(bin);
927 MuonCorr2 = C3muonCorrectionMC[1]->GetBinContent(bin);// -+
928 MuonCorr3 = MuonCorr2;
929 MuonCorr4 = C3muonCorrectionMC[2]->GetBinContent(bin);
930 }
931 double value = n3QS->GetBinContent(bin) * MuonCorr1;
932 value -= ( TERM2_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][1]->GetBinContent(bin) / TwoFrac * MuonCorr2;
933 value -= ( TERM3_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][2]->GetBinContent(bin) / TwoFrac * MuonCorr3;
934 value -= ( TERM4_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][3]->GetBinContent(bin) / TwoFrac * MuonCorr4;
cbf4f1cb 935 //value -= TERM5_3->GetBinContent(bin)*(MuonCorr1 - MuonCorr2 - MuonCorr3 - MuonCorr4);
936 value += 2*TERM5_3->GetBinContent(bin);
264b1bb9 937 //
938 n3QS->SetBinContent(bin, value);
939 c3QS->SetBinContent(bin, value + TERM5_3->GetBinContent(bin));
940 }
941 c3QS->Divide(TERM5_3);
942
943 C3QS->Draw();
944 c3QS->Draw("same");
cbf4f1cb 945 //int lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
946 //int highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
947 //double SF=C3QS->Integral(lowBin, highBin); SF /= TPFullWeight_ThreeParticle[ch1_3]->Integral(lowBin, highBin);
264b1bb9 948 //TPFullWeight_ThreeParticle[ch1_3]->Scale(SF);
cbf4f1cb 949 //
264b1bb9 950 if(SameCharge) TPFullWeight_ThreeParticle[ch1_3]->Draw("same");
951
cbf4f1cb 952 //cout<<TPFullWeight_ThreeParticle[ch1_3]->GetBinContent(2)<<endl;
264b1bb9 953 //TH1D *C3raw = (TH1D*)TERM1_3->Clone();
954 //C3raw->Divide(TERM5_3);
955 //C3raw->SetMarkerColor(4);
956 //C3raw->Draw("same");
957
958 legend2->AddEntry(C3QS,"C_{3}^{QS}","p");
959 legend2->AddEntry(c3QS,"c_{3}^{QS}{2-pion removal}","p");
960 if(SameCharge) legend2->AddEntry(TPFullWeight_ThreeParticle[ch1_3],"C_{3}^{QS} built","l");
961 legend2->Draw("same");
962
cbf4f1cb 963
264b1bb9 964 ///////////////////////////////////////
965 // C3QS built ratio
966 TCanvas *can2_2 = new TCanvas("can2_2", "can2_2",1200,53,700,500);
967 can2_2->SetHighLightColor(2);
968 can2_2->Range(-0.7838115,-1.033258,0.7838115,1.033258);
969 gStyle->SetOptFit(0111);
970 can2_2->SetFillColor(10);
971 can2_2->SetBorderMode(0);
972 can2_2->SetBorderSize(2);
973 can2_2->SetGridx();
974 can2_2->SetGridy();
975 can2_2->SetFrameFillColor(0);
976 can2_2->SetFrameBorderMode(0);
977 can2_2->SetFrameBorderMode(0);
978 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
979 TLegend *legend2_2 = (TLegend*)legend1->Clone();
980 legend2_2->SetX1(0.25); legend2_2->SetX2(0.6); legend2_2->SetY1(0.8); legend2_2->SetY2(0.95);
981
982 TH1D *C3QSratio = (TH1D*)C3QS->Clone();
983
984 //TH1D *Chi2NDF_PointSize_3 = new TH1D("Chi2NDF_PointSize_3","",40,-0.5,39.5);
985 //TH1D *Chi2NDF_FullSize_3 = new TH1D("Chi2NDF_FullSize_3","",40,-0.5,39.5);
986 TH1D *Chi2NDF_PointSize_3 = new TH1D("Chi2NDF_PointSize_3","",100,-0.5,99.5);
987 TH1D *Chi2NDF_FullSize_3 = new TH1D("Chi2NDF_FullSize_3","",100,-0.5,99.5);
988 Chi2NDF_PointSize_3->SetLineColor(4); Chi2NDF_FullSize_3->SetLineColor(2);
989 Chi2NDF_PointSize_3->SetMarkerColor(4); Chi2NDF_FullSize_3->SetMarkerColor(2);
cbf4f1cb 990 Chi2NDF_PointSize_3->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_3->GetYaxis()->SetTitle("#sqrt{#chi^{2}}");
264b1bb9 991
cbf4f1cb 992 if(SameCharge && CollisionType==0){
264b1bb9 993
994 TH1D *tempDen = (TH1D*)TPFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPFullWeight3_Den", 4, 4);
995 TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPNegFullWeight3_Den", 4, 4);
996 tempDen->Add(tempDenNeg);// Add Pos and Neg Den
997
998 for(int binG=5; binG<=104; binG++){// was 44
999 TString *proName=new TString("TPFullWeight3_");
1000 *proName += binG;
1001 TH1D *tempNum = (TH1D*)TPFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY(proName->Data(), binG, binG);
1002 proName->Append("_Neg");
1003 TH1D *tempNumNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY(proName->Data(), binG, binG);
1004 //
1005 // Add Pos and Neg Num
1006 tempNum->Add(tempNumNeg);
1007 //
cbf4f1cb 1008 //tempNum->Add(tempDen);// now added in histogram read section
264b1bb9 1009 tempNum->Divide(tempDen);
cbf4f1cb 1010 //lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
1011 //highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
1012 //SF=C3QS->Integral(lowBin, highBin);
1013 //SF /= tempNum->Integral(lowBin, highBin);
264b1bb9 1014 //tempNum->Scale(SF);
1015
1016 double Chi2=0;
1017 double NDF=0;
1018 for(int binQ3=3; binQ3<=3; binQ3++){
1019 if(tempNum->GetBinContent(binQ3) <=0) continue;
1020 double value = C3QS->GetBinContent(binQ3) - tempNum->GetBinContent(binQ3);
1021 double err = C3QS->GetBinError(binQ3);
1022 if(err <=0) continue;
1023
1024 Chi2 += pow(value / err,2);
1025 //
1026 NDF += 1;
1027 }
1028 //if(binG<25) {Chi2NDF_PointSize_3->SetBinContent(1 + 2*(binG-5), sqrt(Chi2)/NDF); Chi2NDF_PointSize_3->SetBinError(1 + 2*(binG-5), 0.001);}
1029 //else {Chi2NDF_FullSize_3->SetBinContent(1 + 2*(binG-25), sqrt(Chi2)/NDF); Chi2NDF_FullSize_3->SetBinError(1 + 2*(binG-25), 0.001);}
1030 if(binG<55) {Chi2NDF_PointSize_3->SetBinContent(1 + 2*(binG-5), sqrt(Chi2)/NDF); Chi2NDF_PointSize_3->SetBinError(1 + 2*(binG-5), 0.001);}
1031 else {Chi2NDF_FullSize_3->SetBinContent(1 + 2*(binG-55), sqrt(Chi2)/NDF); Chi2NDF_FullSize_3->SetBinError(1 + 2*(binG-55), 0.001);}
1032 //if(NDF>0) cout<<binG<<" "<<sqrt(Chi2)/NDF<<endl;
1033 }
1034 Chi2NDF_PointSize_3->SetMarkerStyle(20); Chi2NDF_FullSize_3->SetMarkerStyle(20);
57e5b35f 1035
264b1bb9 1036 C3QSratio->Divide(TPFullWeight_ThreeParticle[ch1_3]);
1037 //
1038 C3QSratio->SetMinimum(0.94); C3QSratio->SetMaximum(1.02);
1039 C3QSratio->GetYaxis()->SetTitleOffset(1.2);
1040 C3QSratio->GetYaxis()->SetTitle("C_{3}^{QS} / C_{3}(built)");
1041 C3QSratio->Draw();
1042 Unity->Draw("same");
cbf4f1cb 1043
264b1bb9 1044 Chi2NDF_PointSize_3->Draw();
1045 Chi2NDF_FullSize_3->Draw("same");
1046 legend2_2->AddEntry(Chi2NDF_PointSize_3,"R_{coh}=0 (Point Source)","p");
1047 legend2_2->AddEntry(Chi2NDF_FullSize_3,"R_{coh}=R_{ch}","p");
1048 legend2_2->Draw("same");
1049 }
1050
1051 // r3
1052 TH1D *r3;
cbf4f1cb 1053 if(SameCharge && CollisionType==0){
264b1bb9 1054 r3 = (TH1D*)n3QS->Clone();
1055 TPN_ThreeParticle[ch1_3]->Multiply(MRC_SC_3[2]);
1056 r3->Divide(TPN_ThreeParticle[ch1_3]);
1057 r3->GetXaxis()->SetRangeUser(0,0.08);
1058 r3->SetMinimum(0.5); r3->SetMaximum(2.5);
1059 r3->GetYaxis()->SetTitle("r_{3}");
1060 //
cbf4f1cb 1061 r3->Draw();
264b1bb9 1062 //TPN_ThreeParticle[ch1_3]->Draw();
1063 //fTPNRejects3pion->SetLineColor(2);
1064 //fTPNRejects3pion->Draw("same");
1065 }
1066
1067
1068 // Print 3-pion points
1069 for(int bin=1; bin<=10; bin++){
1070 //cout<<C3QS->GetBinContent(bin)<<", ";
1071 //cout<<c3QS->GetBinContent(bin)<<", ";
1072 //cout<<TPFullWeight_ThreeParticle[ch1_3]->GetBinContent(bin)<<", ";
1073 }
1074 //cout<<endl;
1075 for(int bin=1; bin<=10; bin++){
1076 //cout<<C3QS->GetBinError(bin)<<", ";
1077 //cout<<c3QS->GetBinError(bin)<<", ";
1078 }
1079 //cout<<endl;
1080
1081
1082 ////////////////////////////////////////////////////////////////
1083 // 4-pion
1084 cout<<"4-pion section"<<endl;
1085
cbf4f1cb 1086 TCanvas *can3 = new TCanvas("can3", "can3",11,600,700,500);// 60 was 600
264b1bb9 1087 can3->SetHighLightColor(2);
1088 can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1089 gStyle->SetOptFit(0111);
1090 can3->SetFillColor(10);
1091 can3->SetBorderMode(0);
1092 can3->SetBorderSize(2);
1093 can3->SetGridx();
1094 can3->SetGridy();
1095 can3->SetFrameFillColor(0);
1096 can3->SetFrameBorderMode(0);
1097 can3->SetFrameBorderMode(0);
1098 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
1099 TLegend *legend3 = (TLegend*)legend1->Clone();
57e5b35f 1100 legend3->SetX1(0.47); legend3->SetX2(0.98); legend3->SetY1(0.7); legend3->SetY2(0.98);
264b1bb9 1101
1102 TH1D *TERMS_4[13];
1103 for(int index=0; index<13; index++){
1104 TERMS_4[index] = (TH1D*)FourParticle[ch1_4][ch2_4][ch3_4][ch4_4][index]->Clone();
1105 //
1106 TH1D *MRC_temp;
1107 if(SameCharge){
1108 if(index==0) MRC_temp = (TH1D*)MRC_SC_4[0]->Clone();
1109 else if(index<=4) MRC_temp = (TH1D*)MRC_SC_4[1]->Clone();
1110 else if(index<=10) MRC_temp = (TH1D*)MRC_SC_4[2]->Clone();
1111 else if(index==11) MRC_temp = (TH1D*)MRC_SC_4[3]->Clone();
1112 else MRC_temp = (TH1D*)MRC_SC_4[4]->Clone();
1113 }else if(CHARGE==-1 && MixedCharge4pionType==1){
1114 if(index==0) MRC_temp = (TH1D*)MRC_MC1_4[0]->Clone();//0
1115 else if(index!=1 && index<=4) MRC_temp = (TH1D*)MRC_MC1_4[1]->Clone();//1
1116 else if(index==1) MRC_temp = (TH1D*)MRC_MC1_4[2]->Clone();//2
1117 else if(index==7 || index==9 || index==10) MRC_temp = (TH1D*)MRC_MC1_4[3]->Clone();//3
1118 else if(index!=12) MRC_temp = (TH1D*)MRC_MC1_4[4]->Clone();//4
1119 else MRC_temp = (TH1D*)MRC_MC1_4[5]->Clone();//5
1120 }else if(CHARGE==+1 && MixedCharge4pionType==1){
1121 if(index==0) MRC_temp = (TH1D*)MRC_MC1_4[0]->Clone();
1122 else if(index<=3) MRC_temp = (TH1D*)MRC_MC1_4[1]->Clone();
1123 else if(index==4) MRC_temp = (TH1D*)MRC_MC1_4[2]->Clone();
1124 else if(index==5 || index==6 || index==7) MRC_temp = (TH1D*)MRC_MC1_4[3]->Clone();
1125 else if(index!=12) MRC_temp = (TH1D*)MRC_MC1_4[4]->Clone();
1126 else MRC_temp = (TH1D*)MRC_MC1_4[5]->Clone();
1127 }else{// --++
1128 if(index==0) MRC_temp = (TH1D*)MRC_MC2_4[0]->Clone();
1129 else if(index<=4) MRC_temp = (TH1D*)MRC_MC2_4[1]->Clone();
cbf4f1cb 1130 else if(index==5 || index==10) MRC_temp = (TH1D*)MRC_MC2_4[2]->Clone();//2
1131 else if(index<=9) MRC_temp = (TH1D*)MRC_MC2_4[3]->Clone();// 3
1132 else if(index==11) MRC_temp = (TH1D*)MRC_MC2_4[4]->Clone();// 4
264b1bb9 1133 else MRC_temp = (TH1D*)MRC_MC2_4[5]->Clone();
1134 }
1135 //
1136 TERMS_4[index]->Multiply(MRC_temp);
1137 }
1138
1139 TH1D *C4raw=(TH1D*)TERMS_4[2]->Clone();
1140 C4raw->Divide(TERMS_4[12]);
1141 C4raw->GetXaxis()->SetRangeUser(0,0.2);
264b1bb9 1142
cbf4f1cb 1143
264b1bb9 1144 TH1D *N4QS=(TH1D*)TERMS_4[0]->Clone();
1145 //
1146 N4QS->Add(TERMS_4[1], -f43);
1147 N4QS->Add(TERMS_4[2], -f43);
1148 N4QS->Add(TERMS_4[3], -f43);
1149 N4QS->Add(TERMS_4[4], -f43);
1150 //
1151 N4QS->Add(TERMS_4[5], -f42);
1152 N4QS->Add(TERMS_4[6], -f42);
1153 N4QS->Add(TERMS_4[7], -f42);
1154 N4QS->Add(TERMS_4[8], -f42);
1155 N4QS->Add(TERMS_4[9], -f42);
1156 N4QS->Add(TERMS_4[10], -f42);
1157 //
1158 N4QS->Add(TERMS_4[12], -f41);
1159 //
1160 N4QS->Scale(1/f44);
cbf4f1cb 1161 //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Scale(1.005);// K scale variation
264b1bb9 1162 N4QS->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]);
cbf4f1cb 1163 //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]);// K factorization test (MC1)
1164 //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]);// K factorization test (MC2)
264b1bb9 1165
1166 TH1D *C4QS=(TH1D*)N4QS->Clone();
1167 C4QS->GetXaxis()->SetRangeUser(0,0.179);
1168 C4QS->SetMinimum(0.5);
1169 C4QS->SetMarkerColor(4);
1170
cbf4f1cb 1171
264b1bb9 1172 //
1173 /*TH1D *C4QS_basic=(TH1D*)TERMS_4[0]->Clone();
1174 for(int ii=1; ii<=C4QS_basic->GetNbinsX(); ii++){
1175 double value = C4QS_basic->GetBinContent(ii) - TERMS_4[12]->GetBinContent(ii);
1176 value /= f44;
1177 value += TERMS_4[12]->GetBinContent(ii);
1178 C4QS_basic->SetBinContent(ii, value);
1179 }
1180 C4QS_basic->Divide(TERMS_4[12]);
1181 //C4QS_basic->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]);
1182 C4QS_basic->SetMarkerColor(1);
1183 //C4QS_basic->Draw("same");
1184 */
1185 //
1186
1187 TH1D *C4_pieces[12];
1188 for(int i=0; i<12; i++){
1189 C4_pieces[i]=(TH1D*)TERMS_4[i]->Clone();
1190 C4_pieces[i]->Divide(TERMS_4[12]);
1191
1192 if(i==0) C4_pieces[i]->SetMarkerColor(1);
1193 else if(i<5) C4_pieces[i]->SetMarkerColor(2);
1194 else if(i<10) C4_pieces[i]->SetMarkerColor(3);
1195 else C4_pieces[i]->SetMarkerColor(4);
1196 }
1197 C4_pieces[0]->GetXaxis()->SetRangeUser(0,0.2);
1198 C4_pieces[0]->SetMinimum(0.5); C4_pieces[0]->SetMaximum(3);
1199 //C4_pieces[0]->Draw();
1200 //C4_pieces[1]->Draw("same");
1201 //C4_pieces[5]->Draw("same");
1202 //C4_pieces[11]->Draw("same");
1203 //
1204
1205 TH1D *C4_2pairs_FromSquares = (TH1D*)TERMS_4[5]->Clone();
1206 C4_2pairs_FromSquares->Divide(TERMS_4[12]);
1207 for(int bin=1; bin<=C4_2pairs_FromSquares->GetNbinsX(); bin++){
1208 double value = C4_2pairs_FromSquares->GetBinContent(bin);
1209 value -= 1;
1210 value *= value;
1211 value += 1;
1212 C4_2pairs_FromSquares->SetBinContent(bin, value);
1213 }
1214 C4_2pairs_FromSquares->SetMarkerColor(6);
1215 //C4_2pairs_FromSquares->Draw("same");
1216 //
1217 TH1D *C4QS_2pairs=(TH1D*)TERMS_4[11]->Clone();
1218 C4QS_2pairs->Add(TERMS_4[12], -pow(1-TwoFrac,2));// N1^4
1219 C4QS_2pairs->Add(TERMS_4[5], -2*TwoFrac*(1-TwoFrac));// N2 * N1^2
1220 C4QS_2pairs->Scale(1/pow(TwoFrac,2));
1221 C4QS_2pairs->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]);
1222 C4QS_2pairs->Divide(TERMS_4[12]);
1223 C4QS_2pairs->SetMarkerColor(6);
1224 //C4QS_2pairs->Draw("same");
1225 //
1226
1227
1228
1229 TH1D *n4QS = (TH1D*)N4QS->Clone();
1230 TH1D *c4QS = (TH1D*)N4QS->Clone();
1231 TH1D *c4QS_RemovalStage1 = (TH1D*)N4QS->Clone();
1232 TH1D *c4QS_RemovalStage2 = (TH1D*)N4QS->Clone();
1233 TH1D *c4QS_RemovalStage3 = (TH1D*)N4QS->Clone();
1234
1235
1236 for(int bin=1; bin<=n4QS->GetNbinsX(); bin++){
1237 if(n4QS->GetBinContent(bin)==0) continue;
1238 bool exitCode=kFALSE;
1239 for(int ii=0; ii<13; ii++) {if(TERMS_4[ii]->GetBinContent(bin) < 1) exitCode=kTRUE;}
1240 if(exitCode) continue;
1241 //cout<<TERMS_4[0]->GetBinContent(bin)<<endl;
1242 //
1243 double N4QSvalue = N4QS->GetBinContent(bin);
1244 double SubtractionTerm[11] = {0};
1245 double ErrorTerm[12] = {0};
1246 ErrorTerm[0] = n4QS->GetBinError(bin);
1247 //
1248 // 3-pion terms
1249 if(SameCharge || MixedCharge4pionType==1){
1250 SubtractionTerm[0] = (TERMS_4[1]->GetBinContent(bin) - f32*TERMS_4[5]->GetBinContent(bin) - f32*TERMS_4[6]->GetBinContent(bin) - f32*TERMS_4[8]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;// 5 6 8
1251 SubtractionTerm[0] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]->GetBinContent(bin);
1252 ErrorTerm[1] = sqrt(TERMS_4[1]->GetBinContent(bin) + f32*f32*TERMS_4[5]->GetBinContent(bin) + f32*f32*TERMS_4[6]->GetBinContent(bin) + f32*f32*TERMS_4[8]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1253 ErrorTerm[1] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]->GetBinContent(bin);
1254 //
1255 SubtractionTerm[1] = (TERMS_4[2]->GetBinContent(bin) - f32*TERMS_4[5]->GetBinContent(bin) - f32*TERMS_4[7]->GetBinContent(bin) - f32*TERMS_4[9]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1256 SubtractionTerm[1] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][2]->GetBinContent(bin);
1257 ErrorTerm[2] = sqrt(TERMS_4[2]->GetBinContent(bin) + f32*f32*TERMS_4[5]->GetBinContent(bin) + f32*f32*TERMS_4[7]->GetBinContent(bin) + f32*f32*TERMS_4[9]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1258 ErrorTerm[2] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][2]->GetBinContent(bin);
1259 //
1260 SubtractionTerm[2] = (TERMS_4[3]->GetBinContent(bin) - f32*TERMS_4[6]->GetBinContent(bin) - f32*TERMS_4[7]->GetBinContent(bin) - f32*TERMS_4[10]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1261 SubtractionTerm[2] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][3]->GetBinContent(bin);
1262 ErrorTerm[3] = sqrt(TERMS_4[3]->GetBinContent(bin) + f32*f32*TERMS_4[6]->GetBinContent(bin) + f32*f32*TERMS_4[7]->GetBinContent(bin) + f32*f32*TERMS_4[10]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1263 ErrorTerm[3] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][3]->GetBinContent(bin);
1264 //
1265 SubtractionTerm[3] = (TERMS_4[4]->GetBinContent(bin) - f32*TERMS_4[8]->GetBinContent(bin) - f32*TERMS_4[9]->GetBinContent(bin) - f32*TERMS_4[10]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1266 SubtractionTerm[3] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][4]->GetBinContent(bin);
1267 ErrorTerm[4] = sqrt(TERMS_4[4]->GetBinContent(bin) + f32*f32*TERMS_4[8]->GetBinContent(bin) + f32*f32*TERMS_4[9]->GetBinContent(bin) + f32*f32*TERMS_4[10]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1268 ErrorTerm[4] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][4]->GetBinContent(bin);
1269 //
1270 }
1271 // 2-pion terms
cbf4f1cb 1272 SubtractionTerm[4] = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin);
264b1bb9 1273 ErrorTerm[5] = sqrt( TERMS_4[5]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac;
1274 //
1275 SubtractionTerm[5] = ( TERMS_4[6]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][6]->GetBinContent(bin) / TwoFrac;
1276 ErrorTerm[6] = sqrt( TERMS_4[6]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][6]->GetBinContent(bin) / TwoFrac;
1277 //
1278 SubtractionTerm[6] = ( TERMS_4[7]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][7]->GetBinContent(bin) / TwoFrac;
1279 ErrorTerm[7] = sqrt( TERMS_4[7]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][7]->GetBinContent(bin) / TwoFrac;
1280 //
1281 SubtractionTerm[7] = ( TERMS_4[8]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][8]->GetBinContent(bin) / TwoFrac;
1282 ErrorTerm[8] = sqrt( TERMS_4[8]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][8]->GetBinContent(bin) / TwoFrac;
1283 //
1284 SubtractionTerm[8] = ( TERMS_4[9]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][9]->GetBinContent(bin) / TwoFrac;
1285 ErrorTerm[9] = sqrt( TERMS_4[9]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][9]->GetBinContent(bin) / TwoFrac;
1286 //
1287 SubtractionTerm[9] = ( TERMS_4[10]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][10]->GetBinContent(bin) / TwoFrac;
1288 ErrorTerm[10] = sqrt( TERMS_4[10]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][10]->GetBinContent(bin) / TwoFrac;
1289 //
1290
1291 //
1292 // 2 2-pion terms: cos(q12*x12)*cos(q34*x34)
1293 SubtractionTerm[10] = TERMS_4[11]->GetBinContent(bin);// N22
1294 SubtractionTerm[10] -= 2*(1-TwoFrac) * TERMS_4[5]->GetBinContent(bin);
1295 SubtractionTerm[10] += pow(1-TwoFrac,2) * TERMS_4[12]->GetBinContent(bin);
264b1bb9 1296 SubtractionTerm[10] /= pow(TwoFrac,2);
1297 SubtractionTerm[10] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]->GetBinContent(bin);
cbf4f1cb 1298 SubtractionTerm[10] *= C4muonCorrectionSC[3]->GetBinContent(bin);
1299 double intermed = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac;
1300 intermed *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin);
1301 intermed = intermed * C4muonCorrectionSC[2]->GetBinContent(bin);
1302 SubtractionTerm[10] -= 2*intermed;
264b1bb9 1303 SubtractionTerm[10] += 2*TERMS_4[12]->GetBinContent(bin);
1304 //
1305 ErrorTerm[11] = TERMS_4[11]->GetBinContent(bin);
1306 ErrorTerm[11] += pow(1-TwoFrac,4) * TERMS_4[12]->GetBinContent(bin);
1307 ErrorTerm[11] += pow(2*(1-TwoFrac)*TwoFrac,2) * TERMS_4[5]->GetBinContent(bin);
1308 ErrorTerm[11] = sqrt(ErrorTerm[11]);
1309 ErrorTerm[11] /= pow(TwoFrac,2);
1310 ErrorTerm[11] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]->GetBinContent(bin);
1311 ErrorTerm[11] = pow(ErrorTerm[11],2);
1312 ErrorTerm[11] += pow( 2*sqrt(TERMS_4[5]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac, 2);
1313 ErrorTerm[11] += 4*TERMS_4[12]->GetBinContent(bin);
1314 ErrorTerm[11] = sqrt(ErrorTerm[11]);
1315
1316 //
1317 double FinalValue[4]={0};
1318 double FinalValue_e[4]={0};
cbf4f1cb 1319 if(SameCharge) {
1320 //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1321 N4QSvalue *= C4muonCorrectionSC[0]->GetBinContent(bin);
1322 ErrorTerm[0] *= C4muonCorrectionSC[0]->GetBinContent(bin);
1323 }else if(MixedCharge4pionType==1) {
1324 //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1325 N4QSvalue *= C4muonCorrectionMC1[0]->GetBinContent(bin);
1326 ErrorTerm[0] *= C4muonCorrectionMC1[0]->GetBinContent(bin);
1327 }else {
1328 //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1329 N4QSvalue *= C4muonCorrectionMC2[0]->GetBinContent(bin);
1330 ErrorTerm[0] *= C4muonCorrectionMC2[0]->GetBinContent(bin);
1331 }
264b1bb9 1332 for(int ii=0; ii<4; ii++) {
1333 FinalValue[ii] = N4QSvalue;
1334 FinalValue_e[ii] = pow(ErrorTerm[0],2);
1335 }
1336
1337 if(SameCharge) {
cbf4f1cb 1338 //FinalValue[0] -= 4*(SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[1]->GetBinContent(bin);
1339 //FinalValue[0] += 6*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[2]->GetBinContent(bin);
1340 FinalValue[0] -= 4*(SubtractionTerm[0] * C4muonCorrectionSC[1]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1341 FinalValue[0] += 6*(SubtractionTerm[4] * C4muonCorrectionSC[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1342 FinalValue[0] -= 3*(SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin));
1343
264b1bb9 1344 FinalValue_e[0] += 4*pow(ErrorTerm[1]*C4muonCorrectionSC[1]->GetBinContent(bin),2);
1345 FinalValue_e[0] += 6*pow(ErrorTerm[5]*C4muonCorrectionSC[2]->GetBinContent(bin),2);
1346 FinalValue_e[0] += 3*pow(ErrorTerm[11]*C4muonCorrectionSC[3]->GetBinContent(bin),2);
1347 FinalValue_e[0] += (4*C4muonCorrectionSC[1]->GetBinContent(bin) + 6*C4muonCorrectionSC[2]->GetBinContent(bin) + 3*C4muonCorrectionSC[3]->GetBinContent(bin)) * TERMS_4[12]->GetBinContent(bin);
1348 //
1349 FinalValue[1] -= 6*SubtractionTerm[4] - 6*TERMS_4[12]->GetBinContent(bin);// 2-pion removal
1350 FinalValue_e[1] += 6*pow(ErrorTerm[5],2) + 6*TERMS_4[12]->GetBinContent(bin);
1351 FinalValue[2] = FinalValue[1] - 3*SubtractionTerm[10] + 3*TERMS_4[12]->GetBinContent(bin);// further 2-pair removal
1352 FinalValue_e[2] += FinalValue_e[1] + 3*pow(ErrorTerm[11],2) + 3*TERMS_4[12]->GetBinContent(bin);
1353 FinalValue[3] = SubtractionTerm[10];
1354 }else{
1355 if(MixedCharge4pionType==1 && CHARGE==-1){
cbf4f1cb 1356 //FinalValue[0] -= (SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin);// [2] is +++ MuonCorr
1357 FinalValue[0] -= SubtractionTerm[0] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin);// [2] is +++ MuonCorr
1358 FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1359 //
264b1bb9 1360 FinalValue_e[0] += pow(ErrorTerm[1]*C4muonCorrectionMC1[2]->GetBinContent(bin),2);
1361 FinalValue[1] -= 3*SubtractionTerm[4] - 3*TERMS_4[12]->GetBinContent(bin);
1362 FinalValue_e[1] += 3*pow(ErrorTerm[5],2) + 3*TERMS_4[12]->GetBinContent(bin);
1363 }else if(MixedCharge4pionType==1 && CHARGE==+1){
cbf4f1cb 1364 //FinalValue[0] -= (SubtractionTerm[3] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin);
1365 FinalValue[0] -= SubtractionTerm[3] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin);
1366 FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1367 //
264b1bb9 1368 FinalValue_e[0] += pow(ErrorTerm[4]*C4muonCorrectionMC1[2]->GetBinContent(bin),2);
1369 FinalValue[1] -= 3*SubtractionTerm[9] - 3*TERMS_4[12]->GetBinContent(bin);
1370 FinalValue_e[1] += 3*pow(ErrorTerm[10],2) + 3*TERMS_4[12]->GetBinContent(bin);
1371 }else{// --++ case
cbf4f1cb 1372 //FinalValue[0] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin);// old way
1373 FinalValue[0] -= 2*(SubtractionTerm[4] * C4muonCorrectionMC2[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1374 FinalValue[0] -= (SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin));
1375 FinalValue[0] -= 4*(SubtractionTerm[6] * C4muonCorrectionMC2[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
264b1bb9 1376 FinalValue_e[0] += 2*pow(ErrorTerm[5],2) + pow(ErrorTerm[11],2) + 2*TERMS_4[12]->GetBinContent(bin);
1377 //
1378 FinalValue[1] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin);
1379 FinalValue_e[1] += 2*pow(ErrorTerm[5],2) + 2*TERMS_4[12]->GetBinContent(bin);
1380 }
1381 }
1382
1383 for(int ii=0; ii<4; ii++) FinalValue_e[ii] = sqrt(FinalValue_e[ii]);
1384 //
1385
1386 n4QS->SetBinContent(bin, FinalValue[0] - TERMS_4[12]->GetBinContent(bin));
1387 n4QS->SetBinError(bin, sqrt( pow(FinalValue_e[0],2) + TERMS_4[12]->GetBinContent(bin)));
1388 c4QS->SetBinContent(bin, FinalValue[0]);
1389 c4QS->SetBinError(bin, FinalValue_e[0]);
1390 C4QS->SetBinContent(bin, N4QSvalue);
1391 C4QS->SetBinError(bin, ErrorTerm[0]);
1392 //
1393 c4QS_RemovalStage1->SetBinContent(bin, FinalValue[1]);
1394 c4QS_RemovalStage1->SetBinError(bin, FinalValue_e[1]);
1395 //
1396 c4QS_RemovalStage2->SetBinContent(bin, FinalValue[2]);
1397 c4QS_RemovalStage2->SetBinError(bin, FinalValue_e[2]);
1398 //
1399 c4QS_RemovalStage3->SetBinContent(bin, FinalValue[3]);
1400 c4QS_RemovalStage3->SetBinError(bin, FinalValue_e[3]);
1401
1402 }
1403
1404 C4QS->Divide(TERMS_4[12]);
1405 c4QS->Divide(TERMS_4[12]);
1406 c4QS->SetMarkerColor(1); c4QS->SetLineColor(1);
1407 //
1408 c4QS_RemovalStage1->Divide(TERMS_4[12]);
1409 c4QS_RemovalStage1->SetMarkerColor(2); c4QS_RemovalStage1->SetLineColor(2);
1410 c4QS_RemovalStage2->Divide(TERMS_4[12]);
1411 c4QS_RemovalStage2->SetMarkerColor(6); c4QS_RemovalStage2->SetLineColor(6);
1412 c4QS_RemovalStage3->Divide(TERMS_4[12]);
1413 c4QS_RemovalStage3->SetMarkerColor(7); c4QS_RemovalStage3->SetLineColor(7);
1414 //
1415 //
cbf4f1cb 1416
264b1bb9 1417
1418
57e5b35f 1419 if(SameCharge) {
1420 if(CollisionType==0) C4QS->SetMaximum(7);
1421 else C4QS->SetMaximum(12);
1422 }else if(MixedCharge4pionType==1) C4QS->SetMaximum(4);
264b1bb9 1423 else C4QS->SetMaximum(3);
1424
57e5b35f 1425 if(CollisionType!=0) C4QS->GetXaxis()->SetRangeUser(0,0.5);
cbf4f1cb 1426
264b1bb9 1427 C4QS->Draw();
1428 c4QS_RemovalStage1->Draw("same");
1429 if(SameCharge) c4QS_RemovalStage2->Draw("same");
1430 //c4QS_RemovalStage3->Draw("same");
1431 c4QS->Draw("same");
1432
1433
cbf4f1cb 1434 //cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(9)<<endl;
1435
57e5b35f 1436 if(SameCharge && !MCcase) {
cbf4f1cb 1437 //TPFullWeight_FourParticle[ch1_4]->Draw("same");
1438 FullBuildFromFits->Draw("same");
1439 PartialBuildFromFits->Draw("same");
1440 }
264b1bb9 1441 legend3->AddEntry(C4QS,"C_{4}^{QS}","p");
1442 legend3->AddEntry(c4QS_RemovalStage1,"c_{4}^{QS}{2-pion removal}","p");
1443 if(SameCharge) legend3->AddEntry(c4QS_RemovalStage2,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1444 if(SameCharge) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair+3-pion removal}","p");
1445 if(!SameCharge && MixedCharge4pionType==1) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+3-pion removal}","p");
1446 if(!SameCharge && MixedCharge4pionType==2) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1447 //legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1448
1449 if(SameCharge) legend3->AddEntry(TPFullWeight_FourParticle[ch1_4],"C_{4}^{QS} built","l");
1450 legend3->Draw("same");
1451
1452 //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Draw();
1453
cbf4f1cb 1454
1455
1456 //TH1D *c4QSFitNum = (TH1D*)c4QSFitNum2D->ProjectionY("c4QSFitNum",5,5);
1457 //TH1D *c4QSFitDen = (TH1D*)c4QSFitDen2D->ProjectionY("c4QSFitDen",5,5);
1458 //c4QSFitNum->Divide(c4QSFitDen);
1459 //for(int bin=1; bin<20; bin++){
1460 //double fc4 =
1461 //c4QSFitNum->SetBinContent(bin, (c4QSFitNum->GetBinContent(bin) - 1.0)*
1462 //}
1463 //c4QSFitNum->Draw("same");
1464
264b1bb9 1465 //C4QS_Syst->Draw("E2 same");
1466 //C4QS->Draw("same");
1467 //C4raw->Draw();
cbf4f1cb 1468 //TH1D *testTerm = (TH1D*)TERMS_4[11]->Clone();
1469 //testTerm->Divide(TERMS_4[12]);
1470 //testTerm->SetMinimum(1); testTerm->SetMaximum(1.4); testTerm->GetXaxis()->SetRangeUser(0,0.14);
1471 //testTerm->Draw();
1472
264b1bb9 1473 ////////////////////////////////////////////////////////////////
cbf4f1cb 1474 if(SameCharge && CollisionType==0){
264b1bb9 1475 TCanvas *can4 = new TCanvas("can4", "can4",600,600,700,500);
1476 can4->SetHighLightColor(2);
1477 can4->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1478 gStyle->SetOptFit(0111);
1479 can4->SetFillColor(10);
1480 can4->SetBorderMode(0);
1481 can4->SetBorderSize(2);
1482 can4->SetGridx();
1483 can4->SetGridy();
1484 can4->SetFrameFillColor(0);
1485 can4->SetFrameBorderMode(0);
1486 can4->SetFrameBorderMode(0);
1487 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
1488 TLegend *legend3_2 = (TLegend*)legend1->Clone();
1489 legend3_2->SetX1(0.45); legend3_2->SetX2(0.98); legend3_2->SetY1(0.6); legend3_2->SetY2(0.95);
1490 //
1491 //TH1D *Chi2NDF_PointSize_4 = new TH1D("Chi2NDF_PointSize_4","",40,-0.5,39.5);
1492 //TH1D *Chi2NDF_FullSize_4 = new TH1D("Chi2NDF_FullSize_4","",40,-0.5,39.5);
1493 TH1D *Chi2NDF_PointSize_4 = new TH1D("Chi2NDF_PointSize_4","",100,-0.5,99.5);
1494 TH1D *Chi2NDF_FullSize_4 = new TH1D("Chi2NDF_FullSize_4","",100,-0.5,99.5);
1495 Chi2NDF_PointSize_4->SetLineColor(4); Chi2NDF_FullSize_4->SetLineColor(2);
1496 Chi2NDF_PointSize_4->SetMarkerColor(4); Chi2NDF_FullSize_4->SetMarkerColor(2);
cbf4f1cb 1497 Chi2NDF_PointSize_4->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_4->GetYaxis()->SetTitle("#sqrt{#chi^{2}}");
264b1bb9 1498
1499
1500 TH1D *tempDen = (TH1D*)TPFullWeight_FourParticle_2D[ch1_4]->ProjectionY("TPFullWeight4_Den", 4, 4);
1501 TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[ch1_4]->ProjectionY("TPNegFullWeight4_Den", 4, 4);
1502 tempDen->Add(tempDenNeg);// Add Pos and Neg Weight
1503
1504 for(int binG=5; binG<=104; binG++){// 44
1505 TString *proName=new TString("TPFullWeight4_");
1506 *proName += binG;
1507 TH1D *tempNum = (TH1D*)TPFullWeight_FourParticle_2D[ch1_4]->ProjectionY(proName->Data(), binG, binG);
1508 proName->Append("_Neg");
1509 TH1D *tempNumNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[ch1_4]->ProjectionY(proName->Data(), binG, binG);
1510 //
1511 // Add Pos and Neg Weights
1512 tempNum->Add(tempNumNeg);
1513 //
cbf4f1cb 1514 //tempNum->Add(tempDen);// now added in InterpCorr section
264b1bb9 1515 tempNum->Divide(tempDen);
cbf4f1cb 1516
264b1bb9 1517
1518 double Chi2=0;
1519 double NDF=0;
1520 for(int binQ4=4; binQ4<=4; binQ4++){
1521 if(tempNum->GetBinContent(binQ4) <=0) continue;
1522 double value = C4QS->GetBinContent(binQ4) - tempNum->GetBinContent(binQ4);
1523 double err = pow(C4QS->GetBinError(binQ4),2);
cbf4f1cb 1524
264b1bb9 1525 err = sqrt(err);
1526 if(err<=0) continue;
1527
1528 Chi2 += pow(value / err,2);
1529 //
1530 NDF += 1;
1531 }
1532 //if(binG<25) {Chi2NDF_PointSize_4->SetBinContent(1 + 2*(binG-5), sqrt(fabs(Chi2))/NDF); Chi2NDF_PointSize_4->SetBinError(1 + 2*(binG-5), 0.001);}
1533 //else {Chi2NDF_FullSize_4->SetBinContent(1 + 2*(binG-25), sqrt(fabs(Chi2))/NDF); Chi2NDF_FullSize_4->SetBinError(1 + 2*(binG-25), 0.001);}
1534 if(binG<55) {Chi2NDF_PointSize_4->SetBinContent(1 + 2*(binG-5), sqrt(fabs(Chi2))/NDF); Chi2NDF_PointSize_4->SetBinError(1 + 2*(binG-5), 0.001);}
1535 else {Chi2NDF_FullSize_4->SetBinContent(1 + 2*(binG-55), sqrt(fabs(Chi2))/NDF); Chi2NDF_FullSize_4->SetBinError(1 + 2*(binG-55), 0.001);}
1536 }
1537 Chi2NDF_PointSize_4->SetMarkerStyle(20); Chi2NDF_FullSize_4->SetMarkerStyle(20);
1538
1539
1540 Chi2NDF_PointSize_4->Draw();
1541 Chi2NDF_FullSize_4->Draw("same");
1542 legend3_2->AddEntry(Chi2NDF_PointSize_4,"R_{coh}=0 (Point Source)","p");
1543 legend3_2->AddEntry(Chi2NDF_FullSize_4,"R_{coh}=R_{ch}","p");
1544 legend3_2->Draw("same");
1545
1546
1547
1548 }
1549
1550 // Print 4-pion points
1551 for(int bin=1; bin<=12; bin++){
cbf4f1cb 1552 //cout<<C4QS->GetBinContent(bin)<<", ";
264b1bb9 1553 //cout<<c4QS->GetBinContent(bin)<<", ";
1554 //cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(bin)<<", ";
1555 //cout<<C4raw->GetBinContent(bin)<<", ";
1556 //cout<<K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->GetBinContent(bin)<<", ";
1557 }
1558 cout<<endl;
1559 for(int bin=1; bin<=12; bin++){
cbf4f1cb 1560 //cout<<c4QS->GetBinContent(bin)<<", ";
264b1bb9 1561 //cout<<C4QS->GetBinError(bin)<<", ";
1562 //cout<<c4QS->GetBinError(bin)<<", ";
1563 //cout<<C4raw->GetBinError(bin)<<", ";
1564 }
1565 //cout<<endl;
1566 ////////////////////////////////////////////////////////////////
1567 // r4
1568
1569 TF1 *ChaoticLimit_r42 = new TF1("ChaoticLimit_r42","6",0,10);
1570 ChaoticLimit_r42->SetLineStyle(2);
1571 TH1D *r42;
cbf4f1cb 1572 if(SameCharge && CollisionType==0){
264b1bb9 1573 /*TCanvas *can5 = new TCanvas("can5", "can5",1200,600,700,500);
1574 can5->SetHighLightColor(2);
1575 can5->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1576 gStyle->SetOptFit(0111);
1577 can5->SetFillColor(10);
1578 can5->SetBorderMode(0);
1579 can5->SetBorderSize(2);
1580 can5->SetGridx();
1581 can5->SetGridy();
1582 can5->SetFrameFillColor(0);
1583 can5->SetFrameBorderMode(0);
1584 can5->SetFrameBorderMode(0);
1585 gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);*/
1586
1587 r42 = (TH1D*)n4QS->Clone();
1588 TPN_FourParticle[ch1_4]->Multiply(MRC_SC_4[4]);
1589 r42->Divide(TPN_FourParticle[ch1_4]);
1590 r42->GetXaxis()->SetRangeUser(0,0.08);
1591 r42->SetMinimum(0.5); r42->SetMaximum(14);
1592 r42->GetYaxis()->SetTitle("r_{4}{2}");
1593 //
1594 //r42->Draw();
1595 //ChaoticLimit_r42->Draw("same");
1596 //TPN_FourParticle[ch1_4]->Draw();
1597 //fTPNRejects4pion1->SetLineColor(2);
1598 //fTPNRejects4pion1->Draw("same");
1599 }
1600
1601 //double EA = exp(-pow(0.005 * 10/FmToGeV,2)/2.);
1602 //TF1 *C2mag = new TF1("C2mag","pow(1-x,2)*pow([0],2) + 2*x*(1-x)*[0]",0,1);
1603 //C2mag->FixParameter(0, EA);
1604 //TF1 *C4norm = new TF1("C4norm","( 6*pow(1-x,4)*pow([0],4) + 24*x*pow(1-x,3)*pow([0],3) + 4*(2*pow(1-x,3)*pow([0],3) + 6*x*pow(1-x,2)*pow([0],2)) + 3*pow(C2mag,2) + 6*C2mag + 1) / (6*pow(C2mag,2) + 8*pow(C2mag,1.5) + 3*pow(C2mag,2) + 6*C2mag + 1)",0,1);
1605 //C4norm->FixParameter(0, EA);
1606 //C4norm->Draw();
1607
1608 /*TF1 *C2mag = new TF1("C2mag","pow(1-[0],2)*exp(-pow(x*[1],2)/6.) + 2*[0]*(1-[0])*exp(-pow(x*[1],2)/12.)",0,0.2);
1609 C2mag->FixParameter(0, 0.25);
1610 C2mag->FixParameter(1, 10./FmToGeV);
1611 TF1 *C4norm = new TF1("C4norm","( 6*pow(1-[0],4)*exp(-pow(x*[1],2)/3.) + 24*[0]*pow(1-[0],3)*exp(-pow(x*[1],2)/4.) + 4*(2*pow(1-[0],3)*exp(-pow(x*[1],2)/4.) + 6*[0]*pow(1-[0],2)*exp(-pow(x*[1],2)/6.)) + 3*pow(C2mag,2) + 6*C2mag + 1) / (6*pow(C2mag,2) + 8*pow(C2mag,1.5) + 3*pow(C2mag,2) + 6*C2mag + 1)",0,0.2);
1612 C4norm->FixParameter(0, 0.25);
1613 C4norm->FixParameter(1, 10./FmToGeV);*/
1614 //C4norm->Draw();
1615
1616
1617 if(SaveToFile){
1618 TString *savefilename = new TString("OutFiles/OutFile");
1619 if(MCcase) savefilename->Append("MonteCarlo");
1620 //
1621 if(SameCharge) savefilename->Append("SC");
1622 else savefilename->Append("MC");
1623 //
1624 if(!SameCharge) *savefilename += MixedCharge4pionType_def;
1625 //
1626 if(CHARGE==1) savefilename->Append("_Pos_");
1627 else savefilename->Append("_Neg_");
1628 //
1629
1630 savefilename->Append("KT_");
1631 *savefilename += EDbin+1;
1632
1633 savefilename->Append("_M");
1634 *savefilename += Mbin;
1635 savefilename->Append(".root");
1636 TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
1637 //
1638 C2->Write("C2");
1639 C2QS->Write("C2QS");
1640 C3QS->Write("C3QS");
1641 c3QS->Write("c3QS");
1642 C4QS->Write("C4QS");
1643 c4QS->Write("c4QS");
1644 c4QS_RemovalStage1->Write("c4QS_RemovalStage1");
1645 if(SameCharge) {
1646 r3->Write("r3");
1647 r42->Write("r42");
1648 c4QS_RemovalStage2->Write("c4QS_RemovalStage2");
1649 TPFullWeight_ThreeParticle[ch1_3]->Write("C3QS_built");
1650 TPFullWeight_FourParticle[ch1_4]->Write("C4QS_built");
1651 //
1652 TPFullWeight_ThreeParticle_2D[ch1_3]->Write("C3QS_built2D");
1653 TPFullWeight_FourParticle_2D[ch1_4]->Write("C4QS_built2D");
1654 TPNegFullWeight_ThreeParticle_2D[ch1_3]->Write("C3QS_Negbuilt2D");
1655 TPNegFullWeight_FourParticle_2D[ch1_4]->Write("C4QS_Negbuilt2D");
1656 }
1657 //
1658 savefile->Close();
1659
1660 }
1661
1662
1663}
1664
1665//________________________________________________________________________
1666void SetFSICorrelations(){
1667 // read in 2-particle and 3-particle FSI correlations = K2 & K3
1668 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
1669 TFile *fsifile = new TFile("KFile.root","READ");
1670 if(!fsifile->IsOpen()) {
1671 cout<<"No FSI file found!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
1672 }
1673
1674 TH1D *temphistoSS[12];
1675 TH1D *temphistoOS[12];
1676 for(Int_t MB=0; MB<12; MB++) {
1677 TString *nameK2SS = new TString("K2ss_");
1678 *nameK2SS += MB;
1679 temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
1680 //
1681 TString *nameK2OS = new TString("K2os_");
1682 *nameK2OS += MB;
1683 temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
1684 //
1685 fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
1686 fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
1687 fFSIss[MB]->SetDirectory(0);
1688 fFSIos[MB]->SetDirectory(0);
1689 }
1690 //
1691
1692 fsifile->Close();
1693
1694}
1695
1696double Gamov(int chargeProduct, double qinv){
1697
1698 double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
1699
1700 return arg/(exp(arg)-1);
1701}
1702
1703void SetMomResCorrections(){
1704
1705 TString *momresfilename = new TString("MomResFile");
1706 if(FileSetting==7) momresfilename->Append("_10percentIncrease");
cbf4f1cb 1707 if(CollisionType!=0) momresfilename->Append("_ppAndpPb");
264b1bb9 1708 momresfilename->Append(".root");
1709
1710 TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
1711 TString *proName[28];
1712 for(int ii=0; ii<28; ii++){
1713 proName[ii] = new TString("MRC_pro_");
1714 *proName[ii] += ii;
1715 }
1716
1717 MRC_SC_2[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_SC_term1"))->ProjectionY(proName[0]->Data(), RbinMRC, RbinMRC));
1718 MRC_SC_2[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_SC_term2"))->ProjectionY(proName[1]->Data(), RbinMRC, RbinMRC));
1719 MRC_SC_2[0]->SetDirectory(0); MRC_SC_2[1]->SetDirectory(0);
1720 //
1721 MRC_MC_2[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_MC_term1"))->ProjectionY(proName[2]->Data(), RbinMRC, RbinMRC));
1722 MRC_MC_2[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_MC_term2"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC));
1723 MRC_MC_2[0]->SetDirectory(0); MRC_MC_2[1]->SetDirectory(0);
1724 //
1725 MRC_SC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term1"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
1726 MRC_SC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term2"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC));
1727 MRC_SC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term3"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC));
1728 MRC_SC_3[0]->SetDirectory(0); MRC_SC_3[1]->SetDirectory(0); MRC_SC_3[2]->SetDirectory(0);
1729 //
1730 MRC_MC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term1"))->ProjectionY(proName[7]->Data(), RbinMRC, RbinMRC));
1731 MRC_MC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term2"))->ProjectionY(proName[8]->Data(), RbinMRC, RbinMRC));
1732 MRC_MC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term3"))->ProjectionY(proName[9]->Data(), RbinMRC, RbinMRC));
1733 MRC_MC_3[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term4"))->ProjectionY(proName[10]->Data(), RbinMRC, RbinMRC));
1734 MRC_MC_3[0]->SetDirectory(0); MRC_MC_3[1]->SetDirectory(0); MRC_MC_3[2]->SetDirectory(0); MRC_MC_3[3]->SetDirectory(0);
1735 //
1736 MRC_SC_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term1"))->ProjectionY(proName[11]->Data(), RbinMRC, RbinMRC));
1737 MRC_SC_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term2"))->ProjectionY(proName[12]->Data(), RbinMRC, RbinMRC));
1738 MRC_SC_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term3"))->ProjectionY(proName[13]->Data(), RbinMRC, RbinMRC));
1739 MRC_SC_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term4"))->ProjectionY(proName[14]->Data(), RbinMRC, RbinMRC));
1740 MRC_SC_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term5"))->ProjectionY(proName[15]->Data(), RbinMRC, RbinMRC));
1741 MRC_SC_4[0]->SetDirectory(0); MRC_SC_4[1]->SetDirectory(0); MRC_SC_4[2]->SetDirectory(0); MRC_SC_4[3]->SetDirectory(0);
1742 MRC_SC_4[4]->SetDirectory(0);
1743 //
1744 MRC_MC1_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term1"))->ProjectionY(proName[16]->Data(), RbinMRC, RbinMRC));
1745 MRC_MC1_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term2"))->ProjectionY(proName[17]->Data(), RbinMRC, RbinMRC));
1746 MRC_MC1_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term3"))->ProjectionY(proName[18]->Data(), RbinMRC, RbinMRC));
1747 MRC_MC1_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term4"))->ProjectionY(proName[19]->Data(), RbinMRC, RbinMRC));
1748 MRC_MC1_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term5"))->ProjectionY(proName[20]->Data(), RbinMRC, RbinMRC));
1749 MRC_MC1_4[5] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term6"))->ProjectionY(proName[21]->Data(), RbinMRC, RbinMRC));
1750 MRC_MC1_4[0]->SetDirectory(0); MRC_MC1_4[1]->SetDirectory(0); MRC_MC1_4[2]->SetDirectory(0); MRC_MC1_4[3]->SetDirectory(0);
1751 MRC_MC1_4[4]->SetDirectory(0); MRC_MC1_4[5]->SetDirectory(0);
1752 //
1753 MRC_MC2_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term1"))->ProjectionY(proName[22]->Data(), RbinMRC, RbinMRC));
1754 MRC_MC2_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term2"))->ProjectionY(proName[23]->Data(), RbinMRC, RbinMRC));
1755 MRC_MC2_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term3"))->ProjectionY(proName[24]->Data(), RbinMRC, RbinMRC));
1756 MRC_MC2_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term4"))->ProjectionY(proName[25]->Data(), RbinMRC, RbinMRC));
1757 MRC_MC2_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term5"))->ProjectionY(proName[26]->Data(), RbinMRC, RbinMRC));
1758 MRC_MC2_4[5] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term6"))->ProjectionY(proName[27]->Data(), RbinMRC, RbinMRC));
1759 MRC_MC2_4[0]->SetDirectory(0); MRC_MC2_4[1]->SetDirectory(0); MRC_MC2_4[2]->SetDirectory(0); MRC_MC2_4[3]->SetDirectory(0);
1760 MRC_MC2_4[4]->SetDirectory(0); MRC_MC2_4[5]->SetDirectory(0);
1761
cbf4f1cb 1762 if(!MRC || MCcase_def){
264b1bb9 1763 for(int bin=1; bin<=MRC_SC_2[0]->GetNbinsX(); bin++){
1764 MRC_SC_2[0]->SetBinContent(bin, 1.0); MRC_SC_2[1]->SetBinContent(bin, 1.0);
1765 MRC_MC_2[0]->SetBinContent(bin, 1.0); MRC_MC_2[1]->SetBinContent(bin, 1.0);
1766 }
1767 for(int bin=1; bin<=MRC_SC_3[0]->GetNbinsX(); bin++){
1768 MRC_SC_3[0]->SetBinContent(bin, 1.0); MRC_SC_3[1]->SetBinContent(bin, 1.0); MRC_SC_3[2]->SetBinContent(bin, 1.0);
1769 MRC_MC_3[0]->SetBinContent(bin, 1.0); MRC_MC_3[1]->SetBinContent(bin, 1.0); MRC_MC_3[2]->SetBinContent(bin, 1.0);
1770 MRC_MC_3[3]->SetBinContent(bin, 1.0);
1771 }
1772 for(int bin=1; bin<=MRC_SC_4[0]->GetNbinsX(); bin++){
1773 MRC_SC_4[0]->SetBinContent(bin, 1.0); MRC_SC_4[1]->SetBinContent(bin, 1.0); MRC_SC_4[2]->SetBinContent(bin, 1.0);
1774 MRC_SC_4[3]->SetBinContent(bin, 1.0); MRC_SC_4[4]->SetBinContent(bin, 1.0);
1775 MRC_MC1_4[0]->SetBinContent(bin, 1.0); MRC_MC1_4[1]->SetBinContent(bin, 1.0); MRC_MC1_4[2]->SetBinContent(bin, 1.0);
1776 MRC_MC1_4[3]->SetBinContent(bin, 1.0); MRC_MC1_4[4]->SetBinContent(bin, 1.0); MRC_MC1_4[5]->SetBinContent(bin, 1.0);
1777 //
1778 MRC_MC2_4[0]->SetBinContent(bin, 1.0); MRC_MC2_4[1]->SetBinContent(bin, 1.0); MRC_MC2_4[2]->SetBinContent(bin, 1.0);
1779 MRC_MC2_4[3]->SetBinContent(bin, 1.0); MRC_MC2_4[4]->SetBinContent(bin, 1.0); MRC_MC2_4[5]->SetBinContent(bin, 1.0);
1780 }
1781 }
1782 MomResFile->Close();
1783
1784}
1785
1786
1787double cubicInterpolate (double p[4], double x) {
1788 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
1789}
1790double bicubicInterpolate (double p[4][4], double x, double y) {
1791 double arr[4];
1792 arr[0] = cubicInterpolate(p[0], y);
1793 arr[1] = cubicInterpolate(p[1], y);
1794 arr[2] = cubicInterpolate(p[2], y);
1795 arr[3] = cubicInterpolate(p[3], y);
1796 return cubicInterpolate(arr, x);
1797}
1798double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
1799 double arr[4];
1800 arr[0] = bicubicInterpolate(p[0], y, z);
1801 arr[1] = bicubicInterpolate(p[1], y, z);
1802 arr[2] = bicubicInterpolate(p[2], y, z);
1803 arr[3] = bicubicInterpolate(p[3], y, z);
1804 return cubicInterpolate(arr, x);
1805}
1806float fact(float n){
1807 return (n < 1.00001) ? 1 : fact(n - 1) * n;
1808}
1809//________________________________________________________________________________________
1810void SetMuonCorrections(){
1811 TString *name = new TString("MuonCorrection");
1812 if(FileSetting==8) name->Append("_92percent");
cbf4f1cb 1813 if(CollisionType!=0) name->Append("_ppAndpPb");
264b1bb9 1814 name->Append(".root");
1815 TFile *MuonFile=new TFile(name->Data(),"READ");
1816 TString *proName[22];
1817 for(int ii=0; ii<22; ii++){
1818 proName[ii] = new TString("MuonCorr_pro_");
1819 *proName[ii] += ii;
1820 }
1821 C2muonCorrectionSC = (TH1D*)(((TH2D*)MuonFile->Get("C2muonCorrection_SC"))->ProjectionY(proName[0]->Data(), RbinMRC, RbinMRC));
1822 C2muonCorrectionMC = (TH1D*)(((TH2D*)MuonFile->Get("C2muonCorrection_MC"))->ProjectionY(proName[1]->Data(), RbinMRC, RbinMRC));
1823 WeightmuonCorrection = (TH1D*)(((TH2D*)MuonFile->Get("WeightmuonCorrection"))->ProjectionY(proName[2]->Data(), RbinMRC, RbinMRC));
1824 C2muonCorrectionSC->SetDirectory(0); C2muonCorrectionMC->SetDirectory(0); WeightmuonCorrection->SetDirectory(0);
1825 //
1826 C3muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term1"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC));
1827 C3muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term2"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
1828 C3muonCorrectionMC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term1"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC));
1829 C3muonCorrectionMC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term2"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC));
1830 C3muonCorrectionMC[2] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term3"))->ProjectionY(proName[7]->Data(), RbinMRC, RbinMRC));
1831 C3muonCorrectionSC[0]->SetDirectory(0); C3muonCorrectionSC[1]->SetDirectory(0);
1832 C3muonCorrectionMC[0]->SetDirectory(0); C3muonCorrectionMC[1]->SetDirectory(0); C3muonCorrectionMC[2]->SetDirectory(0);
1833 //
1834 C4muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term1"))->ProjectionY(proName[8]->Data(), RbinMRC, RbinMRC));
1835 C4muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term2"))->ProjectionY(proName[9]->Data(), RbinMRC, RbinMRC));
1836 C4muonCorrectionSC[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term3"))->ProjectionY(proName[10]->Data(), RbinMRC, RbinMRC));
1837 C4muonCorrectionSC[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term4"))->ProjectionY(proName[11]->Data(), RbinMRC, RbinMRC));
1838 C4muonCorrectionSC[0]->SetDirectory(0); C4muonCorrectionSC[1]->SetDirectory(0);
1839 C4muonCorrectionSC[2]->SetDirectory(0); C4muonCorrectionSC[3]->SetDirectory(0);
1840 //
1841 C4muonCorrectionMC1[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term1"))->ProjectionY(proName[12]->Data(), RbinMRC, RbinMRC));
1842 C4muonCorrectionMC1[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term2"))->ProjectionY(proName[13]->Data(), RbinMRC, RbinMRC));
1843 C4muonCorrectionMC1[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term3"))->ProjectionY(proName[14]->Data(), RbinMRC, RbinMRC));
1844 C4muonCorrectionMC1[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term4"))->ProjectionY(proName[15]->Data(), RbinMRC, RbinMRC));
1845 C4muonCorrectionMC1[4] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term5"))->ProjectionY(proName[16]->Data(), RbinMRC, RbinMRC));
1846 C4muonCorrectionMC1[0]->SetDirectory(0); C4muonCorrectionMC1[1]->SetDirectory(0);
1847 C4muonCorrectionMC1[2]->SetDirectory(0); C4muonCorrectionMC1[3]->SetDirectory(0); C4muonCorrectionMC1[4]->SetDirectory(0);
1848 //
1849 C4muonCorrectionMC2[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term1"))->ProjectionY(proName[17]->Data(), RbinMRC, RbinMRC));
1850 C4muonCorrectionMC2[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term2"))->ProjectionY(proName[18]->Data(), RbinMRC, RbinMRC));
1851 C4muonCorrectionMC2[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term3"))->ProjectionY(proName[19]->Data(), RbinMRC, RbinMRC));
1852 C4muonCorrectionMC2[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term4"))->ProjectionY(proName[20]->Data(), RbinMRC, RbinMRC));
1853 C4muonCorrectionMC2[4] = (TH1D*)C4muonCorrectionSC[3]->Clone();
1854 C4muonCorrectionMC2[0]->SetDirectory(0); C4muonCorrectionMC2[1]->SetDirectory(0);
1855 C4muonCorrectionMC2[2]->SetDirectory(0); C4muonCorrectionMC2[3]->SetDirectory(0); C4muonCorrectionMC2[4]->SetDirectory(0);
1856 //
1857 //
cbf4f1cb 1858 if(!MuonCorrection || MCcase_def){
264b1bb9 1859 for(int bin=1; bin<=C2muonCorrectionSC->GetNbinsX(); bin++){
1860 C2muonCorrectionSC->SetBinContent(bin, 1.0);
1861 C2muonCorrectionMC->SetBinContent(bin, 1.0);
1862 WeightmuonCorrection->SetBinContent(bin, 1.0);
1863 }
1864 for(int bin=1; bin<=C3muonCorrectionSC[0]->GetNbinsX(); bin++){
1865 C3muonCorrectionSC[0]->SetBinContent(bin, 1.0); C3muonCorrectionSC[1]->SetBinContent(bin, 1.0);
1866 C3muonCorrectionMC[0]->SetBinContent(bin, 1.0); C3muonCorrectionMC[1]->SetBinContent(bin, 1.0); C3muonCorrectionMC[2]->SetBinContent(bin, 1.0);
1867 }
1868 for(int bin=1; bin<=C4muonCorrectionSC[0]->GetNbinsX(); bin++){
1869 C4muonCorrectionSC[0]->SetBinContent(bin, 1.0); C4muonCorrectionSC[1]->SetBinContent(bin, 1.0); C4muonCorrectionSC[2]->SetBinContent(bin, 1.0);
1870 C4muonCorrectionSC[3]->SetBinContent(bin, 1.0);
1871 C4muonCorrectionMC1[0]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[1]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[2]->SetBinContent(bin, 1.0);
1872 C4muonCorrectionMC1[3]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[4]->SetBinContent(bin, 1.0);
1873 //
1874 C4muonCorrectionMC2[0]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[1]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[2]->SetBinContent(bin, 1.0);
1875 C4muonCorrectionMC2[3]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[4]->SetBinContent(bin, 1.0);
1876 }
1877 }
1878
1879 MuonFile->Close();
1880}
1881//________________________________________________________________________
1882void SetFSIindex(Float_t R){
1883 if(!MCcase_def){
cbf4f1cb 1884 if(CollisionType==0){
264b1bb9 1885 if(Mbin_def==0) fFSIindex = 0;//0-5%
1886 else if(Mbin_def==1) fFSIindex = 1;//5-10%
1887 else if(Mbin_def<=3) fFSIindex = 2;//10-20%
1888 else if(Mbin_def<=5) fFSIindex = 3;//20-30%
1889 else if(Mbin_def<=7) fFSIindex = 4;//30-40%
1890 else if(Mbin_def<=9) fFSIindex = 5;//40-50%
1891 else if(Mbin_def<=12) fFSIindex = 6;//40-50%
1892 else if(Mbin_def<=15) fFSIindex = 7;//40-50%
1893 else if(Mbin_def<=18) fFSIindex = 8;//40-50%
1894 else fFSIindex = 8;//90-100%
1895 }else fFSIindex = 9;// pp and pPb
1896 }else{// FSI binning for MC
1897 if(R>=10.) fFSIindex = 0;
1898 else if(R>=9.) fFSIindex = 1;
1899 else if(R>=8.) fFSIindex = 2;
1900 else if(R>=7.) fFSIindex = 3;
1901 else if(R>=6.) fFSIindex = 4;
1902 else if(R>=5.) fFSIindex = 5;
1903 else if(R>=4.) fFSIindex = 6;
1904 else if(R>=3.) fFSIindex = 7;
1905 else if(R>=2.) fFSIindex = 8;
1906 else fFSIindex = 9;
1907 }
1908}
1909//________________________________________________________________________
1910Float_t FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
1911 // returns 2-particle Coulomb correlations = K2
1912 Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
1913 Int_t qbinH = qbinL+1;
1914 if(qbinL <= 0) return 1.0;
1915 if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) {// Scaled Gamov approximation
1916 int chargeproduct = 1;
1917 if(charge1!=charge2) {
1918 chargeproduct = -1;
1919 Float_t ScaleFac = (fFSIos[fFSIindex]->GetBinContent(fFSIos[fFSIindex]->GetNbinsX()-1) - 1);
1920 ScaleFac /= (Gamov(chargeproduct, fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(fFSIos[fFSIindex]->GetNbinsX()-1)) - 1);
1921 return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1);
1922 }else{
1923 Float_t ScaleFac = (fFSIss[fFSIindex]->GetBinContent(fFSIss[fFSIindex]->GetNbinsX()-1) - 1);
1924 ScaleFac /= (Gamov(chargeproduct, fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(fFSIss[fFSIindex]->GetNbinsX()-1)) - 1);
1925 return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1);
1926 }
1927 }
1928
1929 Float_t slope=0;
1930 if(charge1==charge2){
1931 slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
1932 slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
1933 return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
1934 }else {
1935 slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
1936 slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
1937 return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
1938 }
1939}