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 }