| cbf4f1cb |
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 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame