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