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 }