fd15e966b7d67e349fc94c489ff705d0a99e44fd
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / Fit_c3.C
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 "TProfile3D.h"
19 #include "TMath.h"
20 #include "TText.h"
21 #include "TRandom3.h"
22 #include "TArray.h"
23 #include "TLegend.h"
24 #include "TStyle.h"
25 #include "TMinuit.h"
26 #include "TCanvas.h"
27 #include "TLatex.h"
28 #include "TPaveStats.h"
29 #include "TASImage.h"
30 #include "TGraph.h"
31 #include "TSpline.h"
32 #include "TVirtualFitter.h"
33
34 #define BohrR 1963.6885 // Bohr Radius for pions
35 #define FmToGeV 0.19733 // conversion to fm
36 #define PI 3.1415926
37 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
38
39 using namespace std;
40
41 //
42 int CollisionType_def=0;// PbPb, pPb, pp
43 int FitType=0;// (0)Edgeworth, (1)Laguerre, (2)Levy
44 //
45 int Mbin=0;// 0-9: centrality bin in widths of 5%
46 int CHARGE=-1;// -1 or +1: + or - pions for same-charge case, --+ or -++,  ---+ or -+++
47 //
48 int EDbin=0;// 0 or 1: Kt3 bin
49 double G_def = 0.;// coherent %
50 double Rcoh_def = 1.;// Radius of Gaussian coherent component in fm
51 //
52 bool MRC=1;// Momentum Resolution Corrections?
53 bool MuonCorrection=1;// correct for Muons?
54 //
55 int f_choice=0;// 0(Core/Halo), 1(40fm), 2(70fm), 3(100fm)
56 //
57 //
58 //
59 //
60 bool SaveToFile_def=0;
61 int fFSIindex=0;
62 float TwoFrac;// Lambda parameter
63 float OneFrac;// Lambda^{1/2}
64 float ThreeFrac;// Lambda^{3/2}
65 double Qcut_pp = 0.6;// 0.6
66 double Qcut_PbPb = 0.1;
67 double NormQcutLow_pp = 1.0;
68 double NormQcutHigh_pp = 1.5;
69 double NormQcutLow_PbPb = .15;
70 double NormQcutHigh_PbPb = .2;// was .175
71
72
73 const int BINRANGE_3=60;// q12,q13,q23
74 int BINLIMIT_3;
75 double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
76 double A_3_e[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
77 double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
78 double BinCenters[400];
79 double Chi2_c3;
80 double NFitPoints_c3;
81 double Q3LimitLow;
82 double Q3LimitHigh;
83 void fcn_c3(int&, double*, double&, double[], int);
84
85 int RbinMRC;
86
87 TH1D *fFSIss[12];
88 TH1D *fFSIos[12];
89
90 void SetFSICorrelations();
91 void SetFSIindex(Float_t);
92 Float_t FSICorrelation(Int_t, Int_t, Float_t);
93 void SetMuonCorrections();
94 void SetMomResCorrections();
95 double Gamov(int, double);
96
97 //
98 float fact(float);
99
100
101
102 TH1D *MRC_SC_3[3];
103 TH1D *C3muonCorrectionSC[2];
104
105 double AvgQinvSS[30];
106 double AvgQinvOS[30];
107 double BinCentersQ4[20];
108
109
110
111 void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, double G=G_def, double Rcoh=Rcoh_def){
112   
113   SaveToFile_def=SaveToFile;
114   CollisionType_def=CollisionType;
115   G /= 100.;
116   G_def=G;
117   Rcoh_def=Rcoh;
118
119   TFile *_file0;
120   if(CollisionType==0){// PbPb
121     //_file0 = new TFile("Results/PDC_11h_3Dhistos.root","READ");
122     //_file0 = new TFile("Results/PDC_11h_c3FitBuild.root","READ");
123     _file0 = new TFile("Results/PDC_11h_LowNorm_HighNorm.root","READ");
124   }else if(CollisionType==1){// pPb
125     //_file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm.root","READ");
126     _file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm_HighNorm.root","READ");
127   }else{// pp
128     _file0 = new TFile("Results/PDC_10bcde_kMB_3Dhisto_LowNorm_HighNorm.root","READ");
129   }
130
131   TList *MyList;
132   TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root");
133   if(CollisionType==0) MyList=(TList*)tdir->Get("FourPionOutput_1");
134   else MyList=(TList*)tdir->Get("FourPionOutput_2");
135   //MyList=(TList*)_file0->Get("MyList");
136
137   _file0->Close();
138
139
140   if(CollisionType==0) {Q3LimitLow = 0.01; Q3LimitHigh = 0.08;}// 0.01 and 0.08 
141   else {Q3LimitLow = 0.01; Q3LimitHigh = 0.25;}// 0.01 and 0.25
142   
143   
144   //
145   TwoFrac=0.7;
146   OneFrac = sqrt(TwoFrac);
147   ThreeFrac = pow(TwoFrac, 1.5);
148   
149   //// Core/Halo, 40fm, 70fm, 100fm
150   float ThermShift_f33[4]={0., 0.06933, 0.01637, 0.006326};
151   float ThermShift_f32[4]={0., -0.0185, -0.005301, -0.002286};
152   float ThermShift_f31[4]={0., -0.01382, -0.0004682, 0.0005337};
153   float f33prime = ThreeFrac;
154   float f32prime = TwoFrac*(1-OneFrac);
155   float f31prime = pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2);
156   f33prime += ThermShift_f33[f_choice];
157   f32prime += ThermShift_f32[f_choice];
158   f31prime += ThermShift_f31[f_choice];
159   float f33 = f33prime;
160   float f32 = f32prime/TwoFrac;
161   float f31 = f31prime - 3*((1-TwoFrac)*(1-OneFrac) + ThermShift_f32[f_choice]*(1-TwoFrac)/TwoFrac);
162   //cout<<f33 + 3*f32 + f31<<endl;
163
164  
165   cout<<"Mbin = "<<Mbin<<"   KT3 = "<<EDbin<<endl;
166   
167   if(CollisionType==0){
168     if(Mbin==0) {RbinMRC=10;}
169     else if(Mbin==1) {RbinMRC=9;}
170     else if(Mbin<=3) {RbinMRC=8;}
171     else if(Mbin<=5) {RbinMRC=7;}
172     else {RbinMRC=6;}
173   }else{
174     RbinMRC=2;
175   }
176
177
178   if(CollisionType==0) BINLIMIT_3=20;
179   else BINLIMIT_3=30;
180  
181  // bin centers from QS+FSI
182   double BinCenterPbPbCentral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509};
183   double BinCenterpPbAndpp[40]={0.00359275, 0.016357, 0.0257109, 0.035445, 0.045297, 0.0552251, 0.0651888, 0.0751397, 0.0851088, 0.0951108, 0.105084, 0.115079, 0.12507, 0.135059, 0.145053, 0.155049, 0.16505, 0.175038, 0.185039, 0.195034, 0.205023, 0.215027, 0.225024, 0.235023, 0.245011, 0.255017, 0.265017, 0.275021, 0.285021, 0.295017, 0.305018, 0.315018, 0.325013, 0.335011, 0.345016, 0.355019, 0.365012, 0.375016, 0.385017, 0.395016};
184   if(CollisionType==0){
185     for(int i=0; i<40; i++) BinCenters[i] = BinCenterPbPbCentral[i];
186   }else{
187     for(int i=0; i<40; i++) BinCenters[i] = BinCenterpPbAndpp[i];
188   }
189   
190   // extend BinCenters for high q
191   for(int index=40; index<400; index++){
192     if(CollisionType==0) BinCenters[index] = (index+0.5)*(0.005);
193     else BinCenters[index] = (index+0.5)*(0.010);
194   }
195   // Set 0's to 3-particle fit arrays
196   for(int i=1; i<=BINLIMIT_3; i++){// bin number
197     for(int j=1; j<=BINLIMIT_3; j++){// bin number
198       for(int k=1; k<=BINLIMIT_3; k++){// bin number
199         A_3[i-1][j-1][k-1]=0;
200         A_3_e[i-1][j-1][k-1]=0;
201         B_3[i-1][j-1][k-1]=0;
202       }
203     }
204   }
205
206
207   //
208   TH3D *ThreeParticle[2][2][2][5];// ch1,ch2,ch3,term
209   TProfile3D *K3avg[2][2][2][4];
210   double norm_3[5]={0};
211   //
212
213   gStyle->SetOptStat(0);
214   gStyle->SetOptDate(0);
215   gStyle->SetOptFit(0111);
216
217
218   
219   SetFSIindex(10.);
220   SetFSICorrelations();
221   SetMomResCorrections();
222   SetMuonCorrections();
223   //
224   /////////////////////////////////////////////////////////
225   
226
227  
228   TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
229   //
230
231   cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
232
233   
234   
235   ///////////////////////////////////
236   // Get Histograms
237   
238   for(int term=0; term<5; term++){
239     
240     TString *name3 = new TString("ThreeParticle_Charge1_");
241     *name3 += 0;
242     name3->Append("_Charge2_");
243     *name3 += 0;
244     name3->Append("_Charge3_");
245     *name3 += 0;
246     name3->Append("_M_");
247     *name3 += Mbin;
248     name3->Append("_ED_");
249     *name3 += EDbin;
250     name3->Append("_Term_");
251     *name3 += term+1;
252     
253     TString *nameNorm3=new TString(name3->Data());
254     nameNorm3->Append("_Norm");
255     //
256     TString *nameK3=new TString(name3->Data());
257     nameK3->Append("_Kfactor3D");
258     //
259     name3->Append("_3D");
260     
261     
262     
263     norm_3[term] = ((TH1D*)MyList->FindObject(nameNorm3->Data()))->Integral();
264     ThreeParticle[0][0][0][term] = (TH3D*)MyList->FindObject(name3->Data());
265     ThreeParticle[0][0][0][term]->Sumw2();
266     //if(0==0 && 0==0) cout<<"3-pion norms  "<<norm_3[term]<<endl;
267     ThreeParticle[0][0][0][term]->Scale(norm_3[0]/norm_3[term]);
268     ThreeParticle[0][0][0][term]->SetMarkerStyle(20);
269     ThreeParticle[0][0][0][term]->SetTitle("");
270     //
271     
272     //
273     if(term<4){
274       K3avg[0][0][0][term] = (TProfile3D*)MyList->FindObject(nameK3->Data());
275     }
276     
277     //
278   }// term 
279   
280   
281   
282
283   cout<<"Done getting Histograms"<<endl;
284   
285   TF1 *Unity=new TF1("Unity","1",0,100);
286   Unity->SetLineStyle(2);
287
288
289   int ch1=0,ch2=0,ch3=0;
290   
291   
292   
293   ///////////////////////////////////////////////////////////
294   // 3-pion 
295   cout<<"3-pion section"<<endl;
296  
297   TCanvas *can2 = new TCanvas("can2", "can2",600,53,700,500);
298   can2->SetHighLightColor(2);
299   can2->Range(-0.7838115,-1.033258,0.7838115,1.033258);
300   gStyle->SetOptFit(0111);
301   can2->SetFillColor(10);
302   can2->SetBorderMode(0);
303   can2->SetBorderSize(2);
304   can2->SetGridx();
305   can2->SetGridy();
306   can2->SetFrameFillColor(0);
307   can2->SetFrameBorderMode(0);
308   can2->SetFrameBorderMode(0);
309   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
310  
311   int Q3BINS=12;
312   float Q3HistoLimit=0.12;
313   if(CollisionType!=0){ Q3BINS=60; Q3HistoLimit=0.6;}
314
315   TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,Q3HistoLimit);
316   TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,Q3HistoLimit);
317   TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",Q3BINS,0,Q3HistoLimit);
318   TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",Q3BINS,0,Q3HistoLimit);
319   double c3_e[100]={0};
320   
321   
322   double value_num; 
323   double value_num_e;
324   double N3_QS;
325   double N3_QS_e;
326   
327   for(int i=2; i<=ThreeParticle[0][0][0][0]->GetNbinsX(); i++){// bin number
328     double Q12 = BinCenters[i-1];// true center
329     //int Q12bin = int(Q12/0.01)+1;
330     //
331     for(int j=2; j<=ThreeParticle[0][0][0][0]->GetNbinsY(); j++){// bin number
332       double Q13 = BinCenters[j-1];// true center
333       //int Q13bin = int(Q13/0.01)+1;
334       //
335       for(int k=2; k<=ThreeParticle[0][0][0][0]->GetNbinsZ(); k++){// bin number
336         double Q23 = BinCenters[k-1];// true center
337         //int Q23bin = int(Q23/0.01)+1;
338         //
339                 
340         if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
341         if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
342         
343         double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
344         int Q3bin = c3_hist->GetXaxis()->FindBin(Q3);
345         
346         //
347         double K3 = K3avg[0][0][0][0]->GetBinContent(i,j,k);
348         double K2_12 = K3avg[0][0][0][1]->GetBinContent(i,j,k);
349         double K2_13 = K3avg[0][0][0][2]->GetBinContent(i,j,k);
350         double K2_23 = K3avg[0][0][0][3]->GetBinContent(i,j,k);
351
352         
353         if(K3==0) continue;
354
355         double TERM1=ThreeParticle[ch1][ch2][ch3][0]->GetBinContent(i,j,k);// N3 (3-pion yield per q12,q13,q23 cell). 3-pions from same-event
356         double TERM2=ThreeParticle[ch1][ch2][ch3][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
357         double TERM3=ThreeParticle[ch1][ch2][ch3][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
358         double TERM4=ThreeParticle[ch1][ch2][ch3][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event
359         double TERM5=ThreeParticle[ch1][ch2][ch3][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
360         
361         
362         if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
363         if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
364         
365         //
366         if(MRC){
367           TERM1 *= MRC_SC_3[0]->GetBinContent(MRC_SC_3[0]->GetXaxis()->FindBin(Q3));
368           TERM2 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3));
369           TERM3 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3));
370           TERM4 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3));
371           TERM5 *= MRC_SC_3[2]->GetBinContent(MRC_SC_3[2]->GetXaxis()->FindBin(Q3));
372         }
373         double MuonCorr1=1, MuonCorr2=1, MuonCorr3=1, MuonCorr4=1;
374         if(MuonCorrection){
375           MuonCorr1 = C3muonCorrectionSC[0]->GetBinContent(C3muonCorrectionSC[0]->GetXaxis()->FindBin(Q3));
376           MuonCorr2 = C3muonCorrectionSC[1]->GetBinContent(C3muonCorrectionSC[0]->GetXaxis()->FindBin(Q3));
377           MuonCorr3 = MuonCorr2;
378           MuonCorr4 = MuonCorr2;
379         }
380
381                 
382         
383         // Purify.  Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
384         N3_QS = TERM1;
385         N3_QS -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
386         N3_QS -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
387         N3_QS /= ThreeFrac;
388         N3_QS *= K3;
389         N3_QS *=  MuonCorr1;
390
391         
392         // Isolate 3-pion cumulant
393         value_num = N3_QS;
394         value_num -= (TERM2 - (1-TwoFrac)*TERM5)*K2_12/TwoFrac * MuonCorr2;
395         value_num -= (TERM3 - (1-TwoFrac)*TERM5)*K2_13/TwoFrac * MuonCorr3;
396         value_num -= (TERM4 - (1-TwoFrac)*TERM5)*K2_23/TwoFrac * MuonCorr4;
397         value_num += 2*TERM5;
398         
399         
400         
401         
402         // errors
403         N3_QS_e = TERM1;
404         N3_QS_e += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2);
405         N3_QS_e += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2));
406         N3_QS_e /= pow(ThreeFrac,2);
407         N3_QS_e *= pow(K3,2);
408         //
409         value_num_e = N3_QS_e;
410         value_num_e += (pow(K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)*K2_12/TwoFrac*sqrt(TERM5),2));
411         value_num_e += (pow(K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)*K2_13/TwoFrac*sqrt(TERM5),2));
412         value_num_e += (pow(K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)*K2_23/TwoFrac*sqrt(TERM5),2));
413         value_num_e += pow(2*sqrt(TERM5),2);
414         
415         c3_e[Q3bin-1] += value_num_e + TERM5;// add baseline stat error
416
417
418         // Fill histograms
419         c3_hist->Fill(Q3, value_num + TERM5);// for cumulant correlation function
420         Combinatorics_1d->Fill(Q3, TERM5);
421         
422         //
423         A_3[i-1][j-1][k-1] = value_num + TERM5;
424         B_3[i-1][j-1][k-1] = TERM5;
425         A_3_e[i-1][j-1][k-1] = sqrt(value_num_e + TERM5);
426         //if(i==j && i==k && i==4) cout<<A_3[i-1][j-1][k-1]<<"  "<<B_3[i-1][j-1][k-1]<<"  "<<A_3_e[i-1][j-1][k-1]<<endl;
427         ///////////////////////////////////////////////////////////
428         
429       }
430     }
431   }
432   
433   
434  
435   ////////////////////////////
436
437   // Intermediate steps
438   for(int i=0; i<Q3BINS; i++) {c3_hist->SetBinError(i+1, sqrt(c3_e[i]));}
439
440   c3_hist->Divide(Combinatorics_1d);
441
442   ///////////////////////////////////////////////////////////////////////////////////////////////////
443   
444   
445  
446   c3_hist->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})");
447   c3_hist->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}");
448   c3_hist->GetYaxis()->SetTitleSize(0.045); c3_hist->GetXaxis()->SetTitleSize(0.045);
449   c3_hist->GetYaxis()->SetTitleOffset(1.1);
450   c3_hist->GetXaxis()->SetRangeUser(0,Q3LimitHigh);
451   c3_hist->GetYaxis()->SetRangeUser(0.9,4);
452   c3_hist->SetMarkerStyle(25);
453   c3_hist->SetMarkerColor(2);
454   c3_hist->SetLineColor(2);
455   c3_hist->SetMaximum(3);
456   c3_hist->SetMinimum(.7);
457   c3_hist->Draw();
458   //legend2->AddEntry(c3_hist,"#font[12]{#bf{c}}_{3} (cumulant correlation)","p");
459   
460   
461
462   const int npar_c3=7;
463   TMinuit MyMinuit_c3(npar_c3);
464   MyMinuit_c3.SetFCN(fcn_c3);
465   int ierflg_c3=0;
466   double arglist_c3 = 0;
467   MyMinuit_c3.mnexcm("SET NOWarnings",&arglist_c3,1, ierflg_c3);  
468   arglist_c3 = -1;
469   MyMinuit_c3.mnexcm("SET PRint",&arglist_c3,1, ierflg_c3);
470   //arglist_c3=2;// improve Minimization Strategy (1 is default)
471   //MyMinuit_c3.mnexcm("SET STR",&arglist_c3,1,ierflg_c3);
472   //arglist_c3 = 0;
473   //MyMinuit_c3.mnexcm("SCAN", &arglist_c3,1,ierflg_c3);
474   arglist_c3 = 5000;
475   MyMinuit_c3.mnexcm("MIGRAD", &arglist_c3 ,1,ierflg_c3);
476   
477
478   TF1 *c3Fit1D_Expan;
479   if(FitType==0) {
480     c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2.,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2.,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)) + [5]/(120.*pow(2.,2.5))*(32.*pow(x/sqrt(3.)*[2]/0.19733,5) - 160.*pow(x/sqrt(3.)*[2]/0.19733,3) + 120*x/sqrt(3.)*[2]/0.19733) ,3))",0,1);
481   }else if(FitType==1){
482     c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-[2]*x/0.19733 * sqrt(3.)/2.) * pow(1 + [3]*([2]*x/sqrt(3.)/0.19733 - 1) + [4]/2*(pow([2]*x/sqrt(3.)/0.19733,2) - 4*[2]*x/sqrt(3.)/0.19733 + 2) + [5]/6.*(-pow(x/sqrt(3.)*[1]/0.19733,3) + 9*pow(x/sqrt(3.)*[1]/0.19733,2) - 18*x/sqrt(3.)*[1]/0.19733 + 6),3))",0,1);
483   }else{
484     c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733, [3])))",0,1);
485   }
486   
487
488   double OutputPar_c3[npar_c3]={0};
489   double OutputPar_c3_e[npar_c3]={0};
490   
491   double par_c3[npar_c3];               // the start values
492   double stepSize_c3[npar_c3];          // step sizes 
493   double minVal_c3[npar_c3];            // minimum bound on parameter 
494   double maxVal_c3[npar_c3];            // maximum bound on parameter
495   string parName_c3[npar_c3];
496   //          1.0              1.0              10.              0.              0.              0.             1.5  
497   par_c3[0] = 1.0; par_c3[1] = 1.0; par_c3[2] = 10.; par_c3[3] = 0.; par_c3[4] = 0.; par_c3[5] = 0; par_c3[6] = 1.5;
498   stepSize_c3[0] = 0.01; stepSize_c3[1] = 0.1; stepSize_c3[2] = 0.1; stepSize_c3[3] = 0.01; stepSize_c3[4] = 0.01; stepSize_c3[5] = 0.01; stepSize_c3[6] = 0.1;
499   minVal_c3[0] = 0.8; minVal_c3[1] = 0.2; minVal_c3[2] = 4.; minVal_c3[3] = -2; minVal_c3[4] = -1; minVal_c3[5] = -1; minVal_c3[6] = 0.5;
500   maxVal_c3[0] = 1.1; maxVal_c3[1] = 1000.; maxVal_c3[2] = 100.; maxVal_c3[3] = +2; maxVal_c3[4] = +1; maxVal_c3[5] = +1; maxVal_c3[6] = 2.5;
501   parName_c3[0] = "N"; parName_c3[1] = "#lambda_{3}"; parName_c3[2] = "R_{inv}"; parName_c3[3] = "coeff_{3}"; parName_c3[4] = "coeff_{4}"; parName_c3[5] = "coeff_{5}"; parName_c3[6] = "#alpha";
502
503   if(CollisionType==0){ 
504     if(FitType!=0) {
505       par_c3[2]=15.; minVal_c3[2] = 8.;
506     }
507   }else{
508     if(FitType==0) {par_c3[2] = 2.0; minVal_c3[2] = 1.0;}
509     else {
510       par_c3[1] = 4.0; minVal_c3[1] = 1.0;
511       //
512       par_c3[2] = 1.3; minVal_c3[2] = .9; maxVal_c3[2] = 10.;
513     }
514   }
515   
516   if(FitType==0) {par_c3[6] = 2.;}
517   if(FitType==1) {par_c3[6] = 1.;}
518   if(FitType==2) {par_c3[3] = 0; par_c3[4] = 0; par_c3[5] = 0;}
519
520   if(FitType==2) {maxVal_c3[1] = 10.;}
521
522   for (int i=0; i<npar_c3; i++){
523     MyMinuit_c3.DefineParameter(i, parName_c3[i].c_str(), par_c3[i], stepSize_c3[i], minVal_c3[i], maxVal_c3[i]);
524   }
525   if(FitType==0 || FitType==1) { MyMinuit_c3.FixParameter(6);}
526   if(FitType==2){
527     MyMinuit_c3.FixParameter(3);
528     MyMinuit_c3.FixParameter(4);
529     MyMinuit_c3.FixParameter(5);
530   }
531   MyMinuit_c3.FixParameter(0);
532   //MyMinuit_c3.FixParameter(1);
533   //MyMinuit_c3.FixParameter(2);
534   //MyMinuit_c3.FixParameter(3);
535   //MyMinuit_c3.FixParameter(4);
536   MyMinuit_c3.FixParameter(5);
537   
538   /////////////////////////////////////////////////////////////
539   // Do the minimization!
540   cout<<"Start Three-d fit"<<endl;
541   MyMinuit_c3.Migrad();// Minuit's best minimization algorithm
542   cout<<"End Three-d fit"<<endl;
543   /////////////////////////////////////////////////////////////
544   MyMinuit_c3.mnexcm("SHOw PARameters", &arglist_c3, 1, ierflg_c3);
545   cout<<"c3 Fit: Chi2 = "<<Chi2_c3<<"   NDF = "<<(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
546   cout<<" Chi2/NDF = "<<Chi2_c3 / (NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
547
548   for(int i=0; i<npar_c3; i++){
549     MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]);
550   }
551   
552   cout<<"Tij Norm = "<<pow(OutputPar_c3[1], 1/3.)<<endl;
553   
554   if(FitType!=2){
555     c3Fit1D_Expan->FixParameter(0,OutputPar_c3[0]);
556     c3Fit1D_Expan->FixParameter(1,OutputPar_c3[1]);
557     c3Fit1D_Expan->FixParameter(2,OutputPar_c3[2]);
558     c3Fit1D_Expan->FixParameter(3,OutputPar_c3[3]);
559     c3Fit1D_Expan->FixParameter(4,OutputPar_c3[4]);
560     c3Fit1D_Expan->FixParameter(5,OutputPar_c3[5]);
561   }else{// Levy
562     c3Fit1D_Expan->FixParameter(0,OutputPar_c3[0]);
563     c3Fit1D_Expan->FixParameter(1,OutputPar_c3[1]);
564     c3Fit1D_Expan->FixParameter(2,OutputPar_c3[2]);
565     c3Fit1D_Expan->FixParameter(3,OutputPar_c3[6]);
566   }
567   c3Fit1D_Expan->SetLineStyle(1);
568   //c3Fit1D_Expan->Draw("same");
569   
570  
571   // project 3D EW fit onto 1D histogram
572   for(int i=2; i<=BINLIMIT_3; i++){// bin number
573     double Q12 = BinCenters[i-1];// true center
574     for(int j=2; j<=BINLIMIT_3; j++){// bin number
575       double Q13 = BinCenters[j-1];// true center
576       for(int k=2; k<=BINLIMIT_3; k++){// bin number
577         double Q23 = BinCenters[k-1];// true center
578         //      
579         double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
580         
581         if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
582         if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
583         
584         double TERM5=ThreeParticle[ch1][ch2][ch3][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
585         if(TERM5==0) continue;
586         
587         
588         if(MRC) TERM5 *= MRC_SC_3[2]->GetBinContent(MRC_SC_3[2]->GetXaxis()->FindBin(Q3));
589         //
590         double radius = OutputPar_c3[2]/FmToGeV;
591         double arg12 = Q12*radius;
592         double arg13 = Q13*radius;
593         double arg23 = Q23*radius;
594         double Expan12=1, Expan13=1, Expan23=1;
595         if(FitType==0){
596           Expan12 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg12,3) - 12*pow(arg12,1));
597           Expan12 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg12,4) -48*pow(arg12,2) + 12);
598           Expan12 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg12,5) - 160.*pow(arg12,3) + 120*arg12);
599           //
600           Expan13 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg13,3) - 12*pow(arg13,1));
601           Expan13 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg13,4) -48*pow(arg13,2) + 12);
602           Expan13 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg13,5) - 160.*pow(arg13,3) + 120*arg13);
603           //
604           Expan23 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg23,3) - 12*pow(arg23,1));
605           Expan23 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg23,4) -48*pow(arg23,2) + 12);
606           Expan23 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg23,5) - 160.*pow(arg23,3) + 120*arg23);
607         }else if(FitType==1){
608           Expan12 += OutputPar_c3[3]*(arg12 - 1);
609           Expan12 += OutputPar_c3[4]/2.*(pow(arg12,2) - 4*arg12 + 2);
610           Expan12 += OutputPar_c3[5]/6.*(-pow(arg12,3) + 9*pow(arg12,2) - 18*arg12 + 6);
611           //
612           Expan13 += OutputPar_c3[3]*(arg13 - 1);
613           Expan13 += OutputPar_c3[4]/2.*(pow(arg13,2) - 4*arg13 + 2);
614           Expan13 += OutputPar_c3[5]/6.*(-pow(arg13,3) + 9*pow(arg13,2) - 18*arg13 + 6);
615           //
616           Expan23 += OutputPar_c3[3]*(arg23 - 1);
617           Expan23 += OutputPar_c3[4]/2.*(pow(arg23,2) - 4*arg23 + 2);
618           Expan23 += OutputPar_c3[5]/6.*(-pow(arg23,3) + 9*pow(arg23,2) - 18*arg23 + 6);
619         }else{}
620         
621         //
622         double t12_coh = exp(-pow(Rcoh/FmToGeV * Q12,2)/2.);
623         double t23_coh = exp(-pow(Rcoh/FmToGeV * Q23,2)/2.);
624         double t13_coh = exp(-pow(Rcoh/FmToGeV * Q13,2)/2.);
625         double C = 2*OutputPar_c3[1] * pow(1-G,3)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg13,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan12*Expan13*Expan23;
626         C += 2*pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg13,OutputPar_c3[6]))/2.)*Expan12*Expan13*t23_coh;
627         C += 2*pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan12*Expan23*t13_coh;
628         C += 2*pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg13,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan13*Expan23*t12_coh;
629         C += 1.0;
630         C *= TERM5;
631         C *= OutputPar_c3[0];
632         //if(Q3<0.018) continue;
633         GenSignalExpected_num->Fill(Q3, C);
634         GenSignalExpected_den->Fill(Q3, TERM5);
635         //if(Q3<0.02) cout<<Q3<<"  "<<TERM5<<endl;
636         ///////////////////////////////////////////////////////////
637         
638         
639         
640       }
641     }
642   }
643   
644
645   GenSignalExpected_num->Sumw2();
646   GenSignalExpected_num->Divide(GenSignalExpected_den);
647  
648   TSpline3 *c3Fit1D_ExpanSpline = new TSpline3(GenSignalExpected_num);
649   c3Fit1D_ExpanSpline->SetLineWidth(2);
650   double xpoints[1000], ypoints[1000];
651   bool splineOnly=kFALSE;
652   for(int iii=0; iii<1000; iii++){
653     xpoints[iii] = 0 + (iii+0.5)*0.001;
654     //ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]);// to skip spline
655     splineOnly=kTRUE;// to skip 1D approximation
656     if(CollisionType==0) splineOnly=kTRUE;
657     if(CollisionType!=0 && xpoints[iii] > 0.06) splineOnly=kTRUE;
658     if(!splineOnly){
659       ypoints[iii] = c3Fit1D_Expan->Eval(xpoints[iii]);
660       if(c3Fit1D_Expan->Eval(xpoints[iii])<2. && fabs(c3Fit1D_Expan->Eval(xpoints[iii])-c3Fit1D_ExpanSpline->Eval(xpoints[iii])) < 0.01) splineOnly=kTRUE;
661     }
662     else {ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]); splineOnly=kTRUE;}
663   }
664   TGraph *gr_c3Spline = new TGraph(1000,xpoints,ypoints);
665   gr_c3Spline->SetLineWidth(2);
666   if(CollisionType==0) gr_c3Spline->SetLineColor(1);
667   if(CollisionType==1) gr_c3Spline->SetLineColor(2);
668   if(CollisionType==2) gr_c3Spline->SetLineColor(4);
669   
670   gr_c3Spline->Draw("c same");
671   
672   /*
673   double ChiSqSum_1D=0, SumNDF_1D=0;
674   for(int bin=1; bin<=300; bin++){
675     double binCenter = c3_hist->GetXaxis()->GetBinCenter(bin);
676     if(binCenter > Q3Limit) continue;
677     if(c3_hist->GetBinError(bin)==0) continue;
678     if(binCenter < 0.01) continue;
679     int grIndex=1;
680     for(int gr=0; gr<999; gr++){
681       if(binCenter > xpoints[gr] && (binCenter < xpoints[gr+1])) {grIndex=gr; break;}
682     }
683
684     ChiSqSum_1D += pow((c3_hist->GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2);
685     //cout<<c3_hist->GetBinContent(bin)<<"  "<<ypoints[grIndex]<<"  "<<c3_hist->GetBinError(bin)<<endl;
686     cout<<pow((c3_hist->GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2)<<endl;
687     SumNDF_1D++;
688   }
689   cout<<"1D Chi2/NDF = "<<ChiSqSum_1D / (SumNDF_1D-5.)<<endl;
690   */
691   
692   
693  
694
695   if(SaveToFile){
696     TString *savefilename = new TString("FitFiles/FitFile_CT");
697     *savefilename += CollisionType;
698     savefilename->Append("_FT");
699     *savefilename += FitType;
700     savefilename->Append("_R");
701     *savefilename += Rcoh;
702     savefilename->Append("_G");
703     *savefilename += int((G+0.001)/0.02);
704     savefilename->Append(".root");
705     TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
706     MyMinuit_c3.Write("MyMinuit_c3");
707     //
708     savefile->Close();
709   }
710   
711   
712 }
713
714 //________________________________________________________________________
715 void SetFSICorrelations(){
716   // read in 2-particle and 3-particle FSI correlations = K2 & K3
717   // 2-particle input histo from file is binned in qinv.  3-particle in qinv of each pair
718   TFile *fsifile = new TFile("KFile.root","READ");
719   if(!fsifile->IsOpen()) {
720     cout<<"No FSI file found!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
721   }
722   
723   TH1D *temphistoSS[12];
724   TH1D *temphistoOS[12];
725   for(Int_t MB=0; MB<12; MB++) {
726     TString *nameK2SS = new TString("K2ss_");
727     *nameK2SS += MB;
728     temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
729     //
730     TString *nameK2OS = new TString("K2os_");
731     *nameK2OS += MB;
732     temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
733     //
734     fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
735     fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
736     fFSIss[MB]->SetDirectory(0);
737     fFSIos[MB]->SetDirectory(0);
738   }
739   //
740   
741   fsifile->Close();
742   
743 }
744
745 double Gamov(int chargeProduct, double qinv){
746   
747   double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
748   
749   return arg/(exp(arg)-1);
750 }
751
752 void SetMomResCorrections(){
753  
754   TString *momresfilename = new TString("MomResFile");
755   if(CollisionType_def!=0) momresfilename->Append("_ppAndpPb");
756   momresfilename->Append(".root");
757   
758   TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
759   
760   TString *proName[28];
761   for(int ii=0; ii<28; ii++){
762     proName[ii] = new TString("MRC_pro_");
763     *proName[ii] += ii;
764   }
765
766   
767   //
768   MRC_SC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term1"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
769   MRC_SC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term2"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC));
770   MRC_SC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term3"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC));
771   MRC_SC_3[0]->SetDirectory(0); MRC_SC_3[1]->SetDirectory(0); MRC_SC_3[2]->SetDirectory(0);
772   //
773   
774   if(!MRC){
775     for(int bin=1; bin<=MRC_SC_3[0]->GetNbinsX(); bin++){
776       MRC_SC_3[0]->SetBinContent(bin, 1.0); MRC_SC_3[1]->SetBinContent(bin, 1.0); MRC_SC_3[2]->SetBinContent(bin, 1.0);
777     }
778     
779   }
780   MomResFile->Close();
781   
782 }
783
784
785 float fact(float n){
786   return (n < 1.00001) ? 1 : fact(n - 1) * n;
787 }
788 //________________________________________________________________________________________
789 void SetMuonCorrections(){
790   TString *name = new TString("MuonCorrection");
791   if(CollisionType_def!=0) name->Append("_ppAndpPb");
792   
793   name->Append(".root");
794   TFile *MuonFile=new TFile(name->Data(),"READ");
795   TString *proName[22];
796   for(int ii=0; ii<22; ii++){
797     proName[ii] = new TString("MuonCorr_pro_");
798     *proName[ii] += ii;
799   }
800  
801   //
802   C3muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term1"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC));
803   C3muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term2"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
804   
805   C3muonCorrectionSC[0]->SetDirectory(0); C3muonCorrectionSC[1]->SetDirectory(0);
806  
807   //
808   if(!MuonCorrection){
809     for(int bin=1; bin<=C3muonCorrectionSC[0]->GetNbinsX(); bin++){
810       C3muonCorrectionSC[0]->SetBinContent(bin, 1.0); C3muonCorrectionSC[1]->SetBinContent(bin, 1.0);
811     }
812   }
813   
814   MuonFile->Close();
815 }
816 //________________________________________________________________________
817 void SetFSIindex(Float_t R){
818   if(CollisionType_def==0){
819     if(Mbin==0) fFSIindex = 0;//0-5%
820     else if(Mbin==1) fFSIindex = 1;//5-10%
821     else if(Mbin<=3) fFSIindex = 2;//10-20%
822     else if(Mbin<=5) fFSIindex = 3;//20-30%
823     else if(Mbin<=7) fFSIindex = 4;//30-40%
824     else if(Mbin<=9) fFSIindex = 5;//40-50%
825     else if(Mbin<=12) fFSIindex = 6;//40-50%
826     else if(Mbin<=15) fFSIindex = 7;//40-50%
827     else if(Mbin<=18) fFSIindex = 8;//40-50%
828     else fFSIindex = 8;//90-100%
829   }else fFSIindex = 9;// pp and pPb
830   
831 }
832 //________________________________________________________________________
833 Float_t FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
834   // returns 2-particle Coulomb correlations = K2
835   Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
836   Int_t qbinH = qbinL+1;
837   if(qbinL <= 0) return 1.0;
838   if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) {// Scaled Gamov approximation 
839     int chargeproduct = 1;
840     if(charge1!=charge2) {
841       chargeproduct = -1;
842       Float_t ScaleFac = (fFSIos[fFSIindex]->GetBinContent(fFSIos[fFSIindex]->GetNbinsX()-1) - 1);
843       ScaleFac /= (Gamov(chargeproduct, fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(fFSIos[fFSIindex]->GetNbinsX()-1)) - 1);
844       return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1); 
845     }else{
846       Float_t ScaleFac = (fFSIss[fFSIindex]->GetBinContent(fFSIss[fFSIindex]->GetNbinsX()-1) - 1);
847       ScaleFac /= (Gamov(chargeproduct, fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(fFSIss[fFSIindex]->GetNbinsX()-1)) - 1);
848       return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1);
849     }
850   }
851   
852   Float_t slope=0;
853   if(charge1==charge2){
854     slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
855     slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
856     return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
857   }else {
858     slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
859     slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
860     return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
861   }
862 }
863 //__________________________________________________________________________
864 void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){
865
866   double q12=0, q13=0, q23=0;
867   double Expan12=0, Expan13=0, Expan23=0;
868   double C=0;
869   double Rch=par[2]/FmToGeV;
870   double SumChi2=0;
871   //double lnL=0, term1=0, term2=0;
872   NFitPoints_c3=0;
873   //double SumChi2_test=0;
874
875   for(int i=0; i<=BINLIMIT_3; i++){// q12
876     for(int j=0; j<=BINLIMIT_3; j++){// q13
877       for(int k=0; k<=BINLIMIT_3; k++){// q23
878         
879         if(B_3[i][j][k] == 0) continue;
880         if(A_3[i][j][k] == 0) continue;
881         if(A_3_e[i][j][k] == 0) continue;
882
883         q12 = BinCenters[i];
884         q13 = BinCenters[j];
885         q23 = BinCenters[k];
886         double q3 = sqrt(pow(q12,2)+pow(q13,2)+pow(q23,2));
887         if(q3 > Q3LimitHigh) continue;
888         if(q3 < Q3LimitLow) continue;
889         //
890         double arg12 = q12*Rch;
891         double arg13 = q13*Rch;
892         double arg23 = q23*Rch;
893         if(FitType==0){// Edgeworth expansion
894           Expan12 = 1;
895           Expan12 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg12,3) - 12*pow(arg12,1));
896           Expan12 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg12,4) -48*pow(arg12,2) + 12);
897           Expan12 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg12,5) - 160.*pow(arg12,3) + 120*arg12);
898           //
899           Expan13 = 1;
900           Expan13 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg13,3) - 12*pow(arg13,1));
901           Expan13 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg13,4) -48*pow(arg13,2) + 12);
902           Expan13 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg13,5) - 160.*pow(arg13,3) + 120*arg13);
903           //
904           Expan23 = 1;
905           Expan23 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg23,3) - 12*pow(arg23,1));
906           Expan23 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg23,4) -48*pow(arg23,2) + 12);
907           Expan23 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg23,5) - 160.*pow(arg23,3) + 120*arg23);
908         }else if(FitType==1){// Laguerre expansion
909           Expan12 = 1;
910           Expan12 += par[3]*(arg12 - 1);
911           Expan12 += par[4]/2.*(pow(arg12,2) - 4*arg12 + 2);
912           Expan12 += par[5]/6.*(-pow(arg12,3) + 9*pow(arg12,2) - 18*arg12 + 6);
913           //
914           Expan13 = 1;
915           Expan13 += par[3]*(arg13 - 1);
916           Expan13 += par[4]/2.*(pow(arg13,2) - 4*arg13 + 2);
917           Expan13 += par[5]/6.*(-pow(arg13,3) + 9*pow(arg13,2) - 18*arg13 + 6);
918           //
919           Expan23 = 1;
920           Expan23 += par[3]*(arg23 - 1);
921           Expan23 += par[4]/2.*(pow(arg23,2) - 4*arg23 + 2);
922           Expan23 += par[5]/6.*(-pow(arg23,3) + 9*pow(arg23,2) - 18*arg23 + 6);
923         }else{Expan12=1.0; Expan13=1.0; Expan23=1.0;}
924         //
925         double t12_coh = exp(-pow(Rcoh_def/FmToGeV * q12,2)/2.);
926         double t23_coh = exp(-pow(Rcoh_def/FmToGeV * q23,2)/2.);
927         double t13_coh = exp(-pow(Rcoh_def/FmToGeV * q13,2)/2.);
928         C = 2*par[1] * pow(1-G_def,3)*exp(-(pow(arg12,par[6])+pow(arg13,par[6])+pow(arg23,par[6]))/2.)*Expan12*Expan13*Expan23;
929         C += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg13,par[6]))/2.)*Expan12*Expan13*t23_coh;
930         C += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg23,par[6]))/2.)*Expan12*Expan23*t13_coh;
931         C += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg13,par[6])+pow(arg23,par[6]))/2.)*Expan13*Expan23*t12_coh;
932         C += 1.0;
933         C *= par[0];// Norm
934         
935
936         double error = pow(A_3_e[i][j][k]/B_3[i][j][k],2);
937         error += pow(sqrt(B_3[i][j][k])*A_3[i][j][k]/pow(B_3[i][j][k],2),2);
938         error = sqrt(error);
939         SumChi2 += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2);
940         
941         //if(q3<0.05) SumChi2_test += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2);
942         //
943         /*if(A_3[i][j][k] >= 1) term1 = C*(A_3[i][j][k]+B_3[i][j][k])/(A_3[i][j][k]*(C+1));
944         else term1 = 0;
945         term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1));
946         
947         if(term1 > 0.0 && term2 > 0.0){
948           lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2);
949         }else if(term1==0 && term2 > 0.0){
950           lnL += B_3[i][j][k]*log(term2);
951         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
952         */
953
954         NFitPoints_c3++;
955         
956       }
957     }
958   }
959   //f = -2.0*lnL;// log-liklihood minimization
960   f = SumChi2;// Chi2 minimization
961   Chi2_c3 = f;
962   //Chi2_c3 = SumChi2_test;
963   
964 }