#include "TCanvas.h"
#include "TLatex.h"
#include "TPaveStats.h"
+#include "TASImage.h"
#define BohrR 1963.6885 // Bohr Radius for pions
#define FmToGeV 0.19733 // conversion to fm
using namespace std;
-int FileSetting=0;// 0(standard), 1(r3 Lambda), 2(TTC), 3(PID), 4(Pos B), 5(Neg B), 6(Gauss), 7(4 kt bins old TTC cuts)
+int FileSetting=0;// 0(standard), 1(r3 Lambda), 2(TTC), 3(PID), 4(Pos B), 5(Neg B), 6(GaussR8to5), 7(GaussR11to6), 8(4 kt bins old TTC cuts), 9 (FB5and7overlap), 10 (muon Variation: -0.5<NsigmaPion<2.0)
bool MCcase_def=kFALSE;// MC data?
int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case
-bool SameCharge_def=kTRUE;// 3-pion same-charge?
+bool SameCharge_def=1;// 3-pion same-charge?
int EDbin_def=0;// 0 or 1: Kt3 bin
//
+bool MuonCorrection=1;// correct for Muons?
bool IncludeG_def=kFALSE;// Include coherence?
bool IncludeEW_def=kTRUE;// Include EdgeWorth coefficients?
bool IncludeEWfromTherm_def=kFALSE;// Include EdgeWorth coefficients from Therminator?
int Mbin_def=0;// 0-9: centrality bin in widths of 5%
int Ktbin_def=10;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase
+double KShift=1.0;// K factor shift for testing (1.0 for default). TURN OFF STRONGSC FOR THIS TESTING
//
//
bool IncludeMJcorrection=kTRUE;// linear Mini-Jet correction for denominator of r3?
+bool IncludeStrongSC=kTRUE;// same-charge strong FSI
bool SaveToFile_def=kFALSE;// Save outputs to file?
-int SourceType=1;// 0=Gaussian, 1=Therminator, 2=Lorentzian (keep at 1 for default)
+int SourceType=0;// 0=Therminator, 1=Therminator with 50fm cut (keep at 0 for default), 2=Gaussian
bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin?
bool GofP=kFALSE;// Include momentum dependence of coherent fraction?
bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states?
TH3D *CoulCorrOmega0OS;
TH1D *CoulCorr2SS;
TH1D *CoulCorr2OS;
-//
-//static double Lednicky_qinv[74];
-//static double Lednicky_CoulStrong[74];
+TF1 *StrongSC;// same-charge pion strong FSI
+
-void ReadCoulCorrections(int, int, int, int);
+void ReadCoulCorrections(int, int, int);
void ReadCoulCorrections_Omega0();
-//void ReadLednickyFile(int);
+void ReadMuonCorrections(int);
void ReadMomResFile(int, double);
double CoulCorr2(int, double);
double CoulCorrOmega0(bool, double, double, double);
double cubicInterpolate(double[4], double);
double bicubicInterpolate(double[4][4], double, double);
double tricubicInterpolate(double[4][4][4], double, double, double);
+void DrawALICELogo(Bool_t, Float_t, Float_t, Float_t, Float_t);
//
void fcnC2_Global(int&, double*, double&, double[], int);
void fcnOSL(int&, double*, double&, double[], int);
+float fact(float);
const int BINRANGE_C2global=20;
double C2ssFitting[BINRANGE_C2global];
double AvgQinvSS[30];
double AvgQinvOS[30];
//
+TH1D *C2muonCorrection_SC;
+TH1D *C2muonCorrection_MC;
+TH1D *WeightmuonCorrection;
+TH1D *C3muonCorrection;
+
+
void Plot_PDCumulants(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool IncludeEWfromTherm=IncludeEWfromTherm_def, bool SameCharge=SameCharge_def, bool IncludeG=IncludeG_def, bool IncludeEW=IncludeEW_def, bool GRS=GRS_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
OneFrac = sqrt(TwoFrac);
ThreeFrac = pow(TwoFrac,3/2.);
+
+ // same-charge pion strong FSI (obtained from ratio of CoulStrong to Coul at R=7 Gaussian from Lednicky's code)
+ StrongSC=new TF1("StrongSC","[0]+[1]*exp(-pow([2]*x,2))",0,50);
+ StrongSC->FixParameter(0,9.99192e-01);
+ StrongSC->FixParameter(1,-8.01113e-03);
+ StrongSC->FixParameter(2,5.35912e-02);
+
-
- if(SourceType==0 && RValue > 8) {cout<<"Radius value too large!!!"<<endl; return;}
+ if(SourceType==2 && RValue > 8) {cout<<"Radius value too large!!!"<<endl; return;}
- cout<<"Mbin = "<<Mbin<<" Kt = "<<Ktbin<<" R input = "<<RValue<<" lambda input = "<<TwoFrac<<endl;
+ cout<<"Mbin = "<<Mbin<<" Kt = "<<Ktbin<<" SameCharge = "<<SameCharge<<endl;
// Core-Halo modeling of single-particle and triple-particle core fraction
- //float AvgN[10]={472.56, 390.824, 319.597, 265.933, 218.308, 178.865, 141.2, 115.88, 87.5556, 69.3235};// (avg total pion mult)/2. 2.0sigma
- /*
- float fc = (1+sqrt(1 + 4*AvgN[Mbin]*TwoFrac*(AvgN[Mbin]-1)))/(2*AvgN[Mbin]);// Two component model (core + non-interacting halo)
- cout<<"Before: "<<OneFrac<<" "<<ThreeFrac<<endl;
- OneFrac = fc;
- ThreeFrac = (fc*AvgN[Mbin]*(fc*AvgN[Mbin]-1)*(fc*AvgN[Mbin]-2))/(AvgN[Mbin]*(AvgN[Mbin]-1)*(AvgN[Mbin]-2));
- cout<<"After: "<<OneFrac<<" "<<ThreeFrac<<" "<<endl;
- */
+ //float AvgN[10]={472.56, 390.824, 319.597, 265.933, 218.308, 178.865, 141.2, 115.88, 87.5556, 69.3235};// 10h (avg total pion mult)/2. 2.0sigma
+ //float AvgN[10]={571.2, 472.7, 400.2, 325.2, 268.721, 220.3, 178.9, 143.4, 113.412, 88.0};// 11h (avg total pion mult)/2. 2.0sigma
+
//
// avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values. From unsmeared HIJING 0-5% with input signal
double AvgQinvSS_temp[30]={0.00421164, 0.00796583, 0.0128473, 0.0178262, 0.0228416, 0.0276507, 0.0325368, 0.0376114, 0.0424707, 0.0476274, 0.0526015, 0.0575645, 0.0625478, 0.0675416, 0.0725359, 0.0775219, 0.0825521, 0.0875211, 0.0925303, 0.0975319, 0.102544, 0.10753, 0.112499, 0.117545, 0.122522, 0.127522, 0.132499, 0.137514, 0.142495, 0.147521};
if(PbPbcase){// PbPb
if(MCcase) {
if(Mbin<=1){
- myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11.root","READ");
+ //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11.root","READ");
//myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11_Smeared.root","READ");
+ //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_KtgenSignal_Rinv11.root","READ");
+ myfile = new TFile("Results/PDC_pureHIJING_12a17a_fix_KtgenSignal_Rinv11.root","READ");
+ //myfile = new TFile("Results/PDC_HIJING_10h8.root","READ");
}else{
myfile = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p68R11_C2plots.root","READ");
}
//if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// v3 Standard
//
//if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins_FullTPC.root","READ");// FullTPC runlist
+ //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_HighNorm.root","READ");// High Norm variation 1.06<qinv<1.1
if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// Standard
if(FileSetting==1) myfile = new TFile("Results/PDC_11h_Lam0p7_and_Lam0p6.root","READ");// Lam=0.7 and Lam=0.6 (3ktbins, 0.035TTC)
- if(FileSetting==2) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// TTC variation
+ if(FileSetting==2) myfile = new TFile("Results/PDC_11h_3sigmaTTC.root","READ");// TTC variation
if(FileSetting==3) myfile = new TFile("Results/PDC_11h_PID1p5.root","READ");// PID variation (0.035TTC)
if(FileSetting==4) myfile = new TFile("Results/PDC_11h_PosB.root","READ");// Positive B field for r3num (0.035TTC)
if(FileSetting==5) myfile = new TFile("Results/PDC_11h_NegB.root","READ");// Negative B field for r3num (0.035TTC)
- if(FileSetting==6) myfile = new TFile("Results/PDC_11h_3sigmaTTCDenextrap_GaussDenextrap.root","READ");// Gaussian, 3sigmaTTC)
- if(FileSetting==7) myfile = new TFile("Results/PDC_11h_4ktbins.root","READ");// 4 kt bins (0.035TTC)
- if(FileSetting==8) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// 2 Kt3 bins
+ if(FileSetting==6) myfile = new TFile("Results/PDC_11h_GaussR8to5.root","READ");// Gaussian
+ if(FileSetting==7) myfile = new TFile("Results/PDC_11h_GaussR11to6.root","READ");// Gaussian
+ if(FileSetting==8) myfile = new TFile("Results/PDC_11h_4ktbins.root","READ");// 4 kt bins (0.035TTC)
+ if(FileSetting==9) myfile = new TFile("Results/PDC_11h_FB5and7overlap_l0p8.root","READ");// FB5and7overlap
+ if(FileSetting==10) myfile = new TFile("Results/PDC_11h_muonVar.root","READ");
}
}else{// pp
cout<<"pp not currently supported"<<endl;
return;
}
- ReadCoulCorrections(SourceType, RValue, bValue, KtbinFSI);
+ ReadCoulCorrections(RValue, bValue, KtbinFSI);
//ReadLednickyFile(RValue);
ReadMomResFile(RValueMomRes, TwoFracMomRes);
ReadCoulCorrections_Omega0();
+ ReadMuonCorrections(Mbin);
//
/////////////////////////////////////////////////////////
-
+
double NormQcutLow;
double NormQcutHigh;
TList *MyList;
if(!MCcase){
TDirectoryFile *tdir = (TDirectoryFile*)myfile->Get("PWGCF.outputChaoticityAnalysis.root");
- if(FileSetting!=1 && FileSetting !=6) MyList=(TList*)tdir->Get("ChaoticityOutput_1");
+ if(FileSetting!=1) MyList=(TList*)tdir->Get("ChaoticityOutput_1");
else MyList=(TList*)tdir->Get("ChaoticityOutput_2");
}else{
MyList=(TList*)myfile->Get("MyList");
}// term
}// sc
-
+
// 3-particle
for(int c3=0; c3<2; c3++){
for(int sc=0; sc<Sbins_3; sc++){
Two_pipi_mp->Divide(Two_ex[0][1][0][1]);
// Normalization range counting. Just to evaluate required normalization width in qinv.
- //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.1))<<endl;
//cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.06), Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.1))<<endl;
//cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.15), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.175))<<endl;
-
+ //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.15), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.175))<<endl;
const int SCOI_2=0;
for(int i=1; i<=thermDen_ss->GetNbinsX(); i++){
if(thermDen_ss->GetBinContent(i) > 0){
- double err = thermNumSq_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i);
- err -= pow(thermNum_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i),2);
- err /= C2Den_ss->GetBinContent(i);
- err = err + pow(sqrt(thermLargeR_ss->GetBinContent(i))/C2Den_ss->GetBinContent(i),2);
- err += pow(sqrt(C2Den_ss->GetBinContent(i))*thermLargeR_ss->GetBinContent(i)/pow(C2Den_ss->GetBinContent(i),2),2);
- err = sqrt(err);
+ double err = thermNumSq_ss->GetBinContent(i)/thermDen_ss->GetBinContent(i);
+ err -= pow(thermNum_ss->GetBinContent(i)/thermDen_ss->GetBinContent(i),2);
+ err /= thermDen_ss->GetBinContent(i);
+ err = sqrt(fabs(err));
C2Therm_ss->SetBinError(i, err);
}
if(thermDen_os->GetBinContent(i) > 0){
- double err = thermNumSq_os->GetBinContent(i)/C2Den_os->GetBinContent(i);
- err -= pow(thermNum_os->GetBinContent(i)/C2Den_os->GetBinContent(i),2);
- err /= C2Den_os->GetBinContent(i);
- err = err + pow(sqrt(thermLargeR_os->GetBinContent(i))/C2Den_os->GetBinContent(i),2);
- err += pow(sqrt(C2Den_os->GetBinContent(i))*thermLargeR_os->GetBinContent(i)/pow(C2Den_os->GetBinContent(i),2),2);
- err = sqrt(err);
+ double err = thermNumSq_os->GetBinContent(i)/thermDen_os->GetBinContent(i);
+ err -= pow(thermNum_os->GetBinContent(i)/thermDen_os->GetBinContent(i),2);
+ err /= thermDen_os->GetBinContent(i);
+ err = sqrt(fabs(err));
C2Therm_os->SetBinError(i, err);
}
}
for(int i=0; i<BINRANGE_C2global; i++){
C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.;
- if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= MomRes_C2_pp[i];
+ if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= (MomRes_C2_pp[i]-1)*MRCShift+1;
C2ssFitting_e[i] = pow(MomRes_C2_pp[i]*sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2);
C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.);
C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2));
C2ssFitting_e[i] += pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]);
C2osFitting[i] = temp_mp->GetBinContent(i+1);
- if(!GeneratedSignal && !MCcase) C2osFitting[i] *= MomRes_C2_mp[i];
+ if(!GeneratedSignal && !MCcase) C2osFitting[i] *= (MomRes_C2_mp[i]-1)*MRCShift+1;
C2osFitting_e[i] = pow(MomRes_C2_mp[i]*temp_mp->GetBinError(i+1),2);
C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1));
C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1));
//
K2SS_e[i] = sqrt(pow((K2SS[i]-1)*0.02,2) + pow((K2SS[i]-Gamov(+1, AvgQinvSS[i]))*0.02,2));
K2OS_e[i] = sqrt(pow((K2OS[i]-1)*0.02,2) + pow((K2OS[i]-Gamov(-1, AvgQinvOS[i]))*0.02,2));
- //K2SS_e[i] = 0.001;
- //K2OS_e[i] = 0.001;
+ //K2SS_e[i] = 0.0001;
+ //K2OS_e[i] = 0.0001;
//
if(iteration<2){
- par[0] = 1.0; par[1]=0.5; par[2]=0.5; par[3]=9.2; par[4] = .1; par[5] = .2; par[6] = .2; par[7] = 0.; par[8] = 0.; par[9] = 0.; par[10] = 0.01;
+ par[0] = 1.0; par[1]=0.5; par[2]=0.5; par[3]=9.2; par[4] = .1; par[5] = .2; par[6] = .0; par[7] = 0.; par[8] = 0.; par[9] = 1.; par[10] = 0.01;
stepSize[0] = 0.01; stepSize[1] = 0.01; stepSize[2] = 0.02; stepSize[3] = 0.2; stepSize[4] = 0.01; stepSize[5] = 0.001; stepSize[6] = 0.001; stepSize[7] = 0.001; stepSize[8]=0.001; stepSize[9]=0.01; stepSize[10]=0.01;
minVal[0] = 0.995; minVal[1] = 0.2; minVal[2] = 0.; minVal[3] = 0.1; minVal[4] = 0.001; minVal[5] = -10.; minVal[6] = -10.; minVal[7] = -10.; minVal[8]=-10; minVal[9] = 0.995; minVal[10] = 0;
- maxVal[0] = 1.1; maxVal[1] = 1.0; maxVal[2] = 0.99; maxVal[3] = 15.; maxVal[4] = 2.; maxVal[5] = 10.; maxVal[6] = 10.; maxVal[7] = 10.; maxVal[8]=10.; maxVal[9]=1.1; maxVal[10]=1.;
+ maxVal[0] = 1.03; maxVal[1] = 1.0; maxVal[2] = 0.99; maxVal[3] = 15.; maxVal[4] = 2.; maxVal[5] = 10.; maxVal[6] = 10.; maxVal[7] = 10.; maxVal[8]=10.; maxVal[9]=1.03; maxVal[10]=1.;
parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh";
parName[5] = "kappa_3"; parName[6] = "kappa_4"; parName[7] = "kappa_5"; parName[8]="kappa_6"; parName[9]="Norm_2"; parName[10]="P_{coh}";
MyMinuit.mnexcm( "RELease", tmp, 1, err );// Rcoh
}
- //MyMinuit.DefineParameter(3, parName[3].c_str(), 10.88, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
- //MyMinuit.DefineParameter(4, parName[4].c_str(), 5., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
- MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
- MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
- //
if(!IncludeEW){
if(!IncludeEWfromTherm){
// kappa3=0.24, kappa4=0.16 for testing
MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2
}
-
+ //MyMinuit.DefineParameter(3, parName[3].c_str(), 10.5, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
+ //MyMinuit.DefineParameter(4, parName[4].c_str(), 5., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
+ //MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
+ //MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
+ MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
+ MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
+ //
int ierflg=0;
double arglist[10];
//arglist[0]=2;// improve Minimization Strategy (1 is default)
}// iteration loop
-
+
TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone();
TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone();
TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone();
TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone();
TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone();
TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone();
-
+
for(int i=0; i<BINRANGE_C2global; i++){
C2_ss->SetBinContent(i+1, C2ssFitting[i]);
C2_os->SetBinContent(i+1, C2osFitting[i]);
C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]);
C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]);
C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]);
- double staterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
- double staterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
- C2_ss->SetBinError(i+1, staterror_ss);
- C2_os->SetBinError(i+1, staterror_os);
+ double Toterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
+ double Toterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
+ C2_ss->SetBinError(i+1, Toterror_ss);
+ C2_os->SetBinError(i+1, Toterror_os);
C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]);
C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]);
C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]);
K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0);
K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0);
}
-
+
C2_ss_Momsys->SetMarkerSize(0);
C2_ss_Momsys->SetFillStyle(1000);
C2_ss_Momsys->SetFillColor(kRed-10);
C2_os->SetMarkerSize(1.5);
C2_os->SetMarkerStyle(25);
C2_os->SetMarkerColor(4);
-
- TF1 *fitC2ss = new TF1("fitC2ss",C2ssFitFunction, 0.005,0.2, npar);
- TF1 *fitC2os = new TF1("fitC2os",C2osFitFunction, 0.005,0.2, npar);
+
+ TF1 *fitC2ss = new TF1("fitC2ss",C2ssFitFunction, 0.001,0.2, npar);
+ TF1 *fitC2os = new TF1("fitC2os",C2osFitFunction, 0.001,0.2, npar);
for(int i=0; i<npar; i++) {
fitC2ss->FixParameter(i,OutputPar[i]);
fitC2ss->SetParError(i,OutputPar_e[i]);
fitC2os->FixParameter(i,OutputPar[i]);
fitC2os->SetParError(i,OutputPar_e[i]);
}
+ TH1D *fitC2ss_h=new TH1D("fitC2ss_h","",C2_ss->GetNbinsX(),0,C2_ss->GetXaxis()->GetBinUpEdge(C2_ss->GetNbinsX()));
+ TH1D *fitC2os_h=new TH1D("fitC2os_h","",C2_os->GetNbinsX(),0,C2_os->GetXaxis()->GetBinUpEdge(C2_os->GetNbinsX()));
+ fitC2ss_h->SetLineWidth(2); fitC2os_h->SetLineWidth(2);
+ fitC2ss_h->SetLineColor(2); fitC2os_h->SetLineColor(4);
+
+ for(int bin=1; bin<=fitC2ss_h->GetNbinsX(); bin++){
+ double qinv=(bin-0.5)*0.005;
+ fitC2ss_h->SetBinContent(bin, fitC2ss->Eval(qinv));
+ fitC2os_h->SetBinContent(bin, fitC2os->Eval(qinv));
+ }
if(!MCcase){
C2_ss->DrawCopy();
C2_os_Ksys->DrawCopy("E2 same");
C2_ss->DrawCopy("same");
C2_os->DrawCopy("same");
- fitC2ss->DrawCopy("same");
fitC2os->SetLineColor(4);
- fitC2os->DrawCopy("same");
+ fitC2ss_h->DrawCopy("C same");
+ fitC2os_h->DrawCopy("C same");
MJ_bkg_ss->SetLineColor(1);
MJ_bkg_ss->Draw("same");
legend1->AddEntry(MJ_bkg_ss,"Non-femto bkg","l");
legend1->Draw("same");
}
+ //C2Therm_ss->DrawCopy();
+ //C2Therm_os->DrawCopy("same");
/*
////////// C2 systematics
TH1D *Coul_Omega0 = new TH1D("Coul_Omega0","",Q3BINS,0,0.2);
TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,0.2);
TH1D *c3_hist_STAR = new TH1D("c3_hist_STAR","",Q3BINS,0,0.2);
- TProfile *MomRes_2d = new TProfile("MomRes_2d","",Q3BINS,0,0.2, 0,20.,"");
- TProfile *MomRes_3d = new TProfile("MomRes_3d","",Q3BINS,0,0.2, 0,20.,"");
+ //TProfile *MomRes_2d = new TProfile("MomRes_2d","",Q3BINS,0,0.2, 0,20.,"");
+ TProfile *MomRes_3d_term1 = new TProfile("MomRes_3d_term1","",Q3BINS,0,0.2, 0,20.,"");
+ TProfile *MomRes_3d_term2 = new TProfile("MomRes_3d_term2","",Q3BINS,0,0.2, 0,20.,"");
+ TProfile *MomRes_3d_term5 = new TProfile("MomRes_3d_term5","",Q3BINS,0,0.2, 0,20.,"");
TH1D *r3_num_Q3 = new TH1D("r3_num_Q3","",Q3BINS,0,0.2);
TH1D *r3_den_Q3 = new TH1D("r3_den_Q3","",Q3BINS,0,0.2);
r3_num_Q3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;}
}
+ cout<<"here"<<endl;
// SC = species combination. Always 0 meaning pi-pi-pi.
int SCBin=0;
-
//
- ReadCoulCorrections(SourceType, RValue, bValue, 10);// switch to full kt range, 10.
+ ReadCoulCorrections(RValue, bValue, 10);// switch to full kt range, 10.
//ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation (STAR method testing)
TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",20,0,0.2);
TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",20,0,0.2);
if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
- //
- double G_12 = Gamov(CP12, Q12);// point-source Coulomb correlation
+ // point-source Coulomb correlation
+ double G_12 = Gamov(CP12, Q12);
double G_13 = Gamov(CP13, Q13);
double G_23 = Gamov(CP23, Q23);
- double K2_12 = CoulCorr2(CP12, Q12);// finite-source Coulomb+Strong correlation from Therminator.
+ // finite-source Coulomb+Strong correlation from Therminator.
+ double K2_12 = CoulCorr2(CP12, Q12);
double K2_13 = CoulCorr2(CP13, Q13);
double K2_23 = CoulCorr2(CP23, Q23);
double K3 = 1.;// 3-body Coulomb+Strong correlation
if(GRS) {// GRS approach
- K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
- if(CHARGE==+1 && !SameCharge) K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
+ if(SameCharge || CHARGE==-1) K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
+ else K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
}else {
- K3 = CoulCorrOmega0(SameCharge, Q12, Q13, Q23);
- if(CHARGE==+1 && !SameCharge) K3 = CoulCorrOmega0(SameCharge, Q23, Q12, Q13);
+ if(SameCharge || CHARGE==-1) K3 = CoulCorrOmega0(SameCharge, Q12, Q13, Q23);
+ else K3 = CoulCorrOmega0(SameCharge, Q23, Q12, Q13);
}
if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;}
if(K3==0) continue;
-
+ if(GeneratedSignal) K3 = K2_12*K2_13*K2_23;// no interpolation for Generated signal
+
double TERM1=Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k);// N3 (3-pion yield per q12,q13,q23 cell). 3-pions from same-event
double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
+ double muonCorr_C3=1.0, muonCorr_C2_12=1.0, muonCorr_C2_13=1.0, muonCorr_C2_23=1.0;
+ double muonCorr_W12=1.0, muonCorr_W13=1.0, muonCorr_W23=1.0;
+ if(MuonCorrection){
+ if(SameCharge){
+ muonCorr_C3 = C3muonCorrection->GetBinContent(Q3bin);
+ muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(Q12bin);
+ muonCorr_C2_13 = C2muonCorrection_SC->GetBinContent(Q13bin);
+ muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(Q23bin);
+ muonCorr_W12 = WeightmuonCorrection->GetBinContent(Q12bin);
+ muonCorr_W13 = WeightmuonCorrection->GetBinContent(Q13bin);
+ muonCorr_W23 = WeightmuonCorrection->GetBinContent(Q23bin);
+ }else{
+ muonCorr_C3 = C3muonCorrection->GetBinContent(Q3bin);
+ if(CHARGE==-1){// pi-pi-pi+
+ muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(Q12bin);
+ muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(Q13bin);
+ muonCorr_C2_23 = C2muonCorrection_MC->GetBinContent(Q23bin);
+ }else{// pi-pi+pi+
+ muonCorr_C2_12 = C2muonCorrection_MC->GetBinContent(Q12bin);
+ muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(Q13bin);
+ muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(Q23bin);
+ }
+ }
+ }
+
if(Q3>0.08 && Q3<0.09){// just for testing
if(Q12>0.08 || Q13>0.08 || Q23>0.08) OutTriplets++;
else InTriplets++;
TERM3 *= (MomRes3d[0][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
TERM4 *= (MomRes3d[0][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
TERM5 *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
- // Triple produce of 1-d momentum resolution corrections (less accurate).
- /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1];
- TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];
- TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term2_pp[k-1];
- TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term1_pp[k-1];
- TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];*/
- //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1]);
- //MomRes_3d->Fill(Q3, MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
+ MomRes_3d_term1->Fill(Q3, MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM1);
+ MomRes_3d_term2->Fill(Q3, MomRes3d[0][1]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM2);
+ MomRes_3d_term5->Fill(Q3, MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM5);
}else{
if(CHARGE==-1){// pi-pi-pi+
TERM1 *= (MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
TERM3 *= (MomRes3d[1][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
TERM4 *= (MomRes3d[1][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
TERM5 *= (MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
- /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1];
- TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];
- TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_mp[k-1];
- TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_mp[k-1];
- TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];*/
}else {// pi-pi+pi+
TERM1 *= (MomRes3d[1][0]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
TERM2 *= (MomRes3d[1][3]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
TERM3 *= (MomRes3d[1][2]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
TERM4 *= (MomRes3d[1][1]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
TERM5 *= (MomRes3d[1][4]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
- /*TERM1 *= MomRes_term1_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_pp[k-1];
- TERM2 *= MomRes_term1_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];
- TERM3 *= MomRes_term2_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_pp[k-1];
- TERM4 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_pp[k-1];
- TERM5 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];*/
}
- //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1]);// always treats 12 as ss pair
- //MomRes_3d->Fill(Q3, MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
+ MomRes_3d_term1->Fill(Q3, MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin), TERM1);
+ MomRes_3d_term2->Fill(Q3, MomRes3d[1][1]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM2);
+ MomRes_3d_term5->Fill(Q3, MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM5);
}
}
//
for(int LamType=0; LamType<2; LamType++){
- if(LamType==1){TwoFrac=TwoFracr3; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);}// r3 assignment
- else {
- TwoFrac = OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);// Assign to C2 global fit values found previously
- }
+ if(LamType==1) TwoFrac=TwoFracr3;// r3
+ else TwoFrac = 0.7;//OutputPar[1];// c3 (newly fixed to 0.7)
+
+ if(FileSetting==9) TwoFrac=0.8;// for FB5and7overlap choose a fixed higher lambda
+
+ OneFrac=pow(TwoFrac,0.5); // 0.495 to bring baseline up
+ ThreeFrac=pow(TwoFrac,1.5);
+
+ // finite-multiplicity method
+ //OneFrac = (1+sqrt(1 + 4*AvgN[Mbin]*TwoFrac*(AvgN[Mbin]-1)))/(2*AvgN[Mbin]);
+ //ThreeFrac = (OneFrac*AvgN[Mbin]*(OneFrac*AvgN[Mbin]-1)*(OneFrac*AvgN[Mbin]-2))/(AvgN[Mbin]*(AvgN[Mbin]-1)*(AvgN[Mbin]-2));
+
// Purify. Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
N3_QS[LamType] = TERM1;
N3_QS[LamType] -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
N3_QS[LamType] -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
N3_QS[LamType] /= ThreeFrac;
N3_QS[LamType] /= K3;
+ N3_QS[LamType] *= muonCorr_C3;
if(LamType==0) num_QS->Fill(Q3, N3_QS[LamType]);
// Isolate 3-pion cumulant
value_num[LamType] = N3_QS[LamType];
- value_num[LamType] -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac;
- value_num[LamType] -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac;
- value_num[LamType] -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac;
+ value_num[LamType] -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac * muonCorr_C2_12;
+ value_num[LamType] -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac * muonCorr_C2_13;
+ value_num[LamType] -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac * muonCorr_C2_23;
value_num[LamType] += 2*TERM5;
+ //value_num[LamType] += 0.004*TERM5;
// r3 denominator
if(LamType==1 && SameCharge) {
double denMJ = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - (MJ_bkg_ss->Eval(Q12)-1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
denMJ *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - (MJ_bkg_ss->Eval(Q13)-1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
denMJ *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - (MJ_bkg_ss->Eval(Q23)-1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
+ // Strong same-charge correction
+ double Coul_12 = K2_12/StrongSC->Eval(Q12*1000/2.);// Coulomb used online
+ double Coul_13 = K2_13/StrongSC->Eval(Q13*1000/2.);// Coulomb used online
+ double Coul_23 = K2_23/StrongSC->Eval(Q23*1000/2.);// Coulomb used online
+ double denCoulOnly = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*Coul_12 - (1-TwoFrac))/(TwoFrac*Coul_12);
+ denCoulOnly *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*Coul_13 - (1-TwoFrac))/(TwoFrac*Coul_13);
+ denCoulOnly *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*Coul_23 - (1-TwoFrac))/(TwoFrac*Coul_23);
+ // K shift testing
+ double K2back_12 = (K2_12-1)/KShift+1;
+ double K2back_13 = (K2_13-1)/KShift+1;
+ double K2back_23 = (K2_23-1)/KShift+1;
+ double denKshiftBack = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*K2back_12 - (1-TwoFrac))/(TwoFrac*K2back_12);
+ denKshiftBack *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*K2back_13 - (1-TwoFrac))/(TwoFrac*K2back_13);
+ denKshiftBack *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*K2back_23 - (1-TwoFrac))/(TwoFrac*K2back_23);
double den_ratio = sqrt(fabs(denMRC2))/sqrt(fabs(denMRC1));
value_den[LamType] *= den_ratio;// apply Momentum Resolution correction systematic variation if any.
double den_ratioMJ = sqrt(fabs(denMJ))/sqrt(fabs(denMRC1));
if(IncludeMJcorrection) value_den[LamType] *= den_ratioMJ;// Non-femto bkg correction
-
+ if(IncludeStrongSC){
+ double den_ratioStrong = sqrt(fabs(denMRC1))/sqrt(fabs(denCoulOnly));
+ value_den[LamType] *= den_ratioStrong;
+ }
+ double den_ratioKShift = sqrt(fabs(denMRC1))/sqrt(fabs(denKshiftBack));
+ value_den[LamType] *= den_ratioKShift;
+ //
+ value_den[LamType] *= sqrt(muonCorr_W12*muonCorr_W13*muonCorr_W23);// muon correction
+ //
//value_den[LamType] -= TPNRejects->GetBinContent(i,j,k);// Testing for 0-5% only to estimate the effect when C2^QS < 1.0
- value_den[LamType] *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;// additional momentum resolution correction
- }
+ value_den[LamType] *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;// momentum resolution correction to combinatorics
+ }// !MCcase
}
-
+
+
// errors
- N3_QS_e[LamType] = TERM1;
- N3_QS_e[LamType] += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2);
- N3_QS_e[LamType] += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2));
+ N3_QS_e[LamType] = fabs(TERM1);
+ N3_QS_e[LamType] += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(fabs(TERM5)),2);
+ N3_QS_e[LamType] += (pow((1-OneFrac),2)*fabs(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(fabs(TERM5)),2));
N3_QS_e[LamType] /= pow(ThreeFrac,2);
N3_QS_e[LamType] /= pow(K3,2);
if(LamType==0) num_QS_e[Q3bin-1] += N3_QS_e[LamType];
- // errors
value_num_e[LamType] = N3_QS_e[LamType];
- value_num_e[LamType] += (pow(1/K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(TERM5),2));
- value_num_e[LamType] += (pow(1/K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(TERM5),2));
- value_num_e[LamType] += (pow(1/K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(TERM5),2));
- value_num_e[LamType] += pow(2*sqrt(TERM5),2);
+ value_num_e[LamType] += (pow(1/K2_12/TwoFrac*sqrt(fabs(TERM2)),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(fabs(TERM5)),2));
+ value_num_e[LamType] += (pow(1/K2_13/TwoFrac*sqrt(fabs(TERM3)),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(fabs(TERM5)),2));
+ value_num_e[LamType] += (pow(1/K2_23/TwoFrac*sqrt(fabs(TERM4)),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(fabs(TERM5)),2));
+ value_num_e[LamType] += pow(2*sqrt(fabs(TERM5)),2);
if(LamType==0) c3_e[Q3bin-1] += value_num_e[LamType] + TERM5;// add baseline stat error
else r3_e[Q3bin-1] += value_num_e[LamType];
+
}
// Fill histograms
double arg12 = Q12*radius_exp;
double arg13 = Q13*radius_exp;
double arg23 = Q23*radius_exp;
- /*Float_t EW12 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
- EW12 += 0.16/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
- Float_t EW13 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
- EW13 += 0.16/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
- Float_t EW23 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
- EW23 += 0.16/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);*/
+ Float_t EW12 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
+ EW12 += 0.45/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
+ Float_t EW13 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
+ EW13 += 0.45/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
+ Float_t EW23 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
+ EW23 += 0.45/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
- Float_t EW12 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
- EW12 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
- Float_t EW13 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
- EW13 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
- Float_t EW23 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
- EW23 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
+ //Float_t EW12 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
+ //EW12 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
+ //Float_t EW13 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
+ //EW13 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
+ //Float_t EW23 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
+ //EW23 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
//
double cumulant_exp=0, term1_exp=0;
if(SameCharge) {
GenSignalExpected_num->Fill(Q3, term1_exp);
GenSignalExpected_den->Fill(Q3, TERM5);
///////////////////////////////////////////////////////////
-
+
}
}
}
Cterm1->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
Cterm1->GetYaxis()->SetTitle("C_{3}");
//Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}");
- Cterm1->SetMinimum(0.97);
- Cterm1->SetMaximum(1.7);// 6.1 for same-charge
+ Cterm1->SetMinimum(0.99);
+ Cterm1->SetMaximum(1.08);// 6.1 for same-charge
Cterm1->GetXaxis()->SetRangeUser(0,.095);
Cterm1->GetYaxis()->SetTitleOffset(1.4);
Cterm1->Draw();
c3_hist->SetLineColor(2);
c3_hist->SetMaximum(3);
c3_hist->SetMinimum(.9);
- if(!MCcase) c3_hist->Draw("same");
+ if(!MCcase && !GeneratedSignal) c3_hist->Draw("same");
//legend2->AddEntry(c3_hist,"#font[12]{c}_{3} (cumulant correlation)","p");
-
+ if(SameCharge){
+ //for(int ii=0; ii<10; ii++) cout<<c3_hist->GetBinContent(ii+1)<<", ";
+ TF1 *c3Fit=new TF1("c3Fit","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,0.2);
+ //TF1 *c3Fit=new TF1("c3Fit","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + (-0.12/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + (0.43/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1);
+ c3Fit->SetParameter(0,1);
+ c3Fit->SetParameter(1,2);
+ c3Fit->SetParameter(2,6);
+ //c3Fit->FixParameter(1,0.8545);
+ //c3Fit->FixParameter(2,8.982);
+ //c3_hist->Fit(c3Fit,"IME","",0,0.11);
+ //cout<<"Chi2/NDF = "<<c3Fit->GetChisquare()/c3Fit->GetNDF()<<endl;
+ }
+
GenSignalExpected_num->SetMarkerStyle(20);
//GenSignalExpected_num->Draw("same");
//legend2->AddEntry(GenSignalExpected_num,"#kappa_{3}=0.24, #kappa_{4}=0.16, #lambda=0.68, R=6 fm","p");
//MomRes_2d->SetTitle("");
//MomRes_2d->Draw();
//legend2->AddEntry(MomRes_2d,"2D: Triple pair product","p");
- //MomRes_3d->Draw("same");
+ MomRes_3d_term2->SetLineColor(2);
+ MomRes_3d_term5->SetLineColor(4);
+ //MomRes_3d_term1->Draw();
+ //MomRes_3d_term2->Draw("same");
+ //MomRes_3d_term5->Draw("same");
//legend2->AddEntry(MomRes_3d,"3D","p");
-
+
//legend2->Draw("same");
-
+
+ //cout<<c3_hist->Integral(8,10)/3.<<endl;
/*
////////// C3 systematics
- // C3 --- base, M0, (0.035 TTC for all)
- //double C3_base[10]={0, 1.63072, 1.6096, 1.43731, 1.28205, 1.17777, 1.11973, 1.07932, 1.05459, 1.04029};
- //double C3_base_e[10]={0, 0.022528, 0.00387504, 0.00115667, 0.000423799, 0.000238973, 0.000143993, 9.71502e-05, 7.02007e-05, 5.31267e-05};
- // C3 --- base, M0, Pos B (0.035 TTC for all)
- //double C3_base[10]={0, 1.62564, 1.60461, 1.438, 1.28234, 1.17827, 1.12009, 1.07953, 1.05474, 1.04037};
- //double C3_base_e[10]={0, 0.0239233, 0.00409002, 0.0012215, 0.000446701, 0.000251769, 0.000151651, 0.000102284, 7.38993e-05, 5.59212e-05};
- // C3 --+ base, M0, (0.035 TTC for all)
- //double C3_base[10]={0, 1.66016, 1.41961, 1.25229, 1.16459, 1.11254, 1.08012, 1.05768, 1.04265, 1.0332};
- //double C3_base_e[10]={0, 0.00539779, 0.00111398, 0.000387926, 0.000192906, 0.00011428, 7.24765e-05, 5.09126e-05, 3.76421e-05, 2.87533e-05};
-
- // C3 --- base, M0, (New Standard: 0.04 TTC )
- //double C3_base[10]={0, 1.63903, 1.60254, 1.43381, 1.2794, 1.17603, 1.11806, 1.07806, 1.05345, 1.03936};
- //double C3_base_e[10]={0, 0.0322796, 0.00438361, 0.00122249, 0.000424557, 0.000233965, 0.000139058, 9.28136e-05, 6.66159e-05, 5.01816e-05};
- // C3 --- base, M1, (New Standard: 0.04 TTC )
- //double C3_base[10]={0, 1.6127, 1.65561, 1.49508, 1.33093, 1.21458, 1.14708, 1.099, 1.06864, 1.05064};
- //double C3_base_e[10]={0, 0.0425447, 0.0061176, 0.00172948, 0.000600699, 0.000329342, 0.000194832, 0.000129427, 9.25599e-05, 6.95395e-05};
- // C3 --- base, M2, (New Standard: 0.04 TTC )
- //double C3_base[10]={0, 1.6509, 1.71863, 1.54652, 1.38092, 1.25226, 1.17549, 1.12068, 1.08408, 1.06236};
- //double C3_base_e[10]={0, 0.0981473, 0.0143699, 0.00404612, 0.00141189, 0.000770764, 0.000453944, 0.000300452, 0.000214068, 0.000160448};
- // C3 --- base, M9, (New Standard: 0.04 TTC )
- //double C3_base[10]={0, 2.41982, 2.18303, 1.93205, 1.80399, 1.60955, 1.49305, 1.38565, 1.29873, 1.23626};
- //double C3_base_e[10]={0, 1.60227, 0.177274, 0.0501455, 0.018456, 0.00998147, 0.00583719, 0.00379465, 0.00264116, 0.0019391};
- //
- // C3 --+ base, M0, (New Standard: 0.04 TTC )
- //double C3_base[10]={0, 1.66087, 1.41943, 1.25081, 1.16313, 1.11143, 1.07917, 1.05701, 1.04215, 1.0328};
- //double C3_base_e[10]={0, 0.00584743, 0.00111278, 0.000374009, 0.00018315, 0.000107523, 6.78669e-05, 4.75116e-05, 3.50489e-05, 2.67328e-05};
+ // C3 --+ base, M0, Kt3 0
+ double C3_base[10]={0, 1.64715, 1.40709, 1.24344, 1.15809, 1.1071, 1.07544, 1.0534, 1.03881, 1.02974};
+ double C3_base_e[10]={0, 0.00348782, 0.000841611, 0.00031699, 0.000163779, 9.95652e-05, 6.43811e-05, 4.60347e-05, 3.45978e-05, 2.68396e-05};
//
// HIJING C3 --- base, M0
//double C3_base[10]={0, 0.970005, 1.00655, 1.00352, 1.00291, 1.00071, 1.0002, 0.999524, 0.999404, 0.999397};
//double C3_base[10]={0, 1.34345, 1.04226, 1.0278, 0.99582, 1.00554, 1.00296, 1.00057, 1.00271, 1.00152};
//double C3_base_e[10]={0, 0.363559, 0.0503531, 0.0170117, 0.00679185, 0.00419035, 0.00264603, 0.00184587, 0.00136663, 0.00104772};
// HIJING C3 --- base, M3
- double C3_base[10]={0, 0.890897, 1.05222, 1.00461, 1.01364, 0.998981, 1.00225, 1.00305, 1.00235, 1.00043};
- double C3_base_e[10]={0, 0.297124, 0.0604397, 0.0195066, 0.00812906, 0.00490835, 0.00310751, 0.00217408, 0.00160575, 0.00123065};
+ //double C3_base[10]={0, 0.890897, 1.05222, 1.00461, 1.01364, 0.998981, 1.00225, 1.00305, 1.00235, 1.00043};
+ //double C3_base_e[10]={0, 0.297124, 0.0604397, 0.0195066, 0.00812906, 0.00490835, 0.00310751, 0.00217408, 0.00160575, 0.00123065};
TH1D *C3baseHisto=(TH1D*)Cterm1->Clone();
for(int ii=0; ii<10; ii++){
C3baseHisto->SetBinContent(ii+1, C3_base[ii]);
C3_Sys->SetMaximum(0.01);
//C3_Sys->Draw();
TF1 *C3lineSys=new TF1("C3lineSys","pol3",0,0.099);
- //C3_Sys->Fit(C3lineSys,"MEQ","",0,0.099);
+ C3_Sys->Fit(C3lineSys,"MEQ","",0,0.099);
if(MCcase){
// Plotting +++ added with --- for HIJING
*/
/*
////////// c3 systematics
- // c3 --- base, M0, (New Standard: 0.04 TTC )
- double c3_base[10]={0, 1.86014, 1.47533, 1.23733, 1.09944, 1.04145, 1.01693, 1.00715, 1.00253, 1.00111};
- double c3_base_e[10]={0, 0.104645, 0.0120917, 0.00333303, 0.00118126, 0.0006692, 0.000405246, 0.000274163, 0.000198507, 0.000150258};
+ // c3 --- base, M0, (0.04 TTC )
+ //double c3_base[10]={0, 1.86014, 1.47533, 1.23733, 1.09944, 1.04145, 1.01693, 1.00715, 1.00253, 1.00111};
+ //double c3_base_e[10]={0, 0.104645, 0.0120917, 0.00333303, 0.00118126, 0.0006692, 0.000405246, 0.000274163, 0.000198507, 0.000150258};
+ // c3 --- base, M0, Kt3 0 (New Standard)
+ double c3_base[10]={0, 1.00906, 1.00013, 1.00203, 1.00001, 1.00017, 1.00001, 0.999721, 0.999778, 0.999844};
+ double c3_base_e[10]={0, 0.00726286, 0.00196376, 0.00081047, 0.00044086, 0.000276571, 0.000182574, 0.000132467, 0.000100557, 7.85009e-05};
cout<<"c3 values:"<<endl;
for(int ii=0; ii<10; ii++){
c3_Sys->Fit(c3lineSys,"MEQ","",0,0.099);
*/
-
- /*
- TPad *pad1 = new TPad("pad1","pad1",0.0,0.0,1.,1.);
+
+ /*TPad *pad1 = new TPad("pad1","pad1",0.0,0.0,1.,1.);
gPad->SetGridx(0);
gPad->SetGridy(0);
gPad->SetTickx();
TF1 *QuarticFit = new TF1("QuarticFit","[0]*(1-[1]*pow(x,4))",0,.1);
QuarticFit->SetParameter(0,2); QuarticFit->SetParameter(1,0);
QuarticFit->SetLineColor(4);
- QuarticFit->SetParName(0,"I(Q^{4})"); QuarticFit->SetParName(1,"a(Q^{4})");
+ QuarticFit->SetParName(0,"I^{Quartic}"); QuarticFit->SetParName(1,"a^{Quartic}");
TF1 *QuadraticFit = new TF1("QuadraticFit","[0]*(1-[1]*pow(x,2))",0,.1);
QuadraticFit->SetParameter(0,2); QuadraticFit->SetParameter(1,0);
- QuadraticFit->SetParName(0,"I(Q^{2})"); QuadraticFit->SetParName(1,"a(Q^{2})");
+ QuadraticFit->SetParName(0,"I^{Quadratic}"); QuadraticFit->SetParName(1,"a^{Quadratic}");
+ TLegend *legend3 = new TLegend(.2,.85,.5,.95,NULL,"brNDC");
+ legend3->SetBorderSize(1);
+ legend3->SetTextSize(.04);// small .03; large .06
+ legend3->SetFillColor(0);
if(SameCharge){
- r3_num_Q3->SetMinimum(0); r3_num_Q3->SetMaximum(2.5);
- r3_num_Q3->GetXaxis()->SetRangeUser(0.01,0.09);
+ r3_num_Q3->SetMinimum(0); r3_num_Q3->SetMaximum(2.5);// 0 to 2.5
+ r3_num_Q3->GetXaxis()->SetRangeUser(0.0,0.1);
r3_num_Q3->Draw();
+
+
+ // HIJING standalone
+ // Kt3 1
+ /*double Gen_r3_m_M0[10]={0, 1.97963, 2.10248, 2.04465, 1.96697, 2.02295, 1.92281, 2.10031, 2.22338, 2.49729};
+ double Gen_r3_m_M0_e[10]={0, 0.132726, 0.0668652, 0.0451347, 0.042838, 0.0546967, 0.0754362, 0.133624, 0.286307, 0.628381};
+ double Gen_r3_p_M0[10]={0, 1.91771, 1.9653, 2.00742, 2.02393, 1.90624, 1.93554, 1.66, 1.79584, 0.301761};
+ double Gen_r3_p_M0_e[10]={0, 0.133654, 0.0667213, 0.0451512, 0.0428925, 0.0547591, 0.0754764, 0.13365, 0.28638, 0.628441};*/
+ // Kt3 2
+ /*double Gen_r3_m_M0[10]={0, 1.12993, 2.09715, 1.91886, 2.08493, 2.10931, 2.00286, 1.99898, 1.78549, 1.91861};
+ double Gen_r3_m_M0_e[10]={0, 0.237903, 0.115454, 0.0725178, 0.0630867, 0.0730326, 0.0886163, 0.12885, 0.213654, 0.379273};
+ double Gen_r3_p_M0[10]={0, 2.05766, 1.97408, 2.05182, 2.02431, 2.11783, 1.93294, 1.97525, 2.21833, 2.31318};
+ double Gen_r3_p_M0_e[10]={0, 0.24238, 0.11515, 0.0725077, 0.0629938, 0.0729842, 0.0886386, 0.12884, 0.213737, 0.379196};*/
+
+
+ // HIJING + ALICE
+ // Kt3 1
+ /*double Gen_r3_m_M0[10]={0, 2.03468, 2.05783, 1.97757, 2.03809, 1.95703, 2.02915, 1.80055, 1.97664, 1.49573};
+ double Gen_r3_m_M0_e[10]={0, 0.284164, 0.0923288, 0.0543837, 0.045952, 0.0565222, 0.0748221, 0.128994, 0.271954, 0.599077};
+ double Gen_r3_p_M0[10]={0, 1.74611, 1.92369, 2.11024, 2.02823, 2.06235, 2.00127, 1.91551, 1.90576, 2.521};
+ double Gen_r3_p_M0_e[10]={0, 0.279962, 0.0928866, 0.0548168, 0.0462253, 0.0569129, 0.0752683, 0.129803, 0.2736, 0.602535};
+ double Gen_r3_m_M1[10]={0, 1.71341, 1.8735, 1.9352, 1.96466, 1.93852, 1.69509, 1.51402, 1.6014, -0.802941};
+ double Gen_r3_m_M1_e[10]={0, 0.334562, 0.108368, 0.0640739, 0.0540579, 0.0663174, 0.0876532, 0.149302, 0.307982, 0.640548};
+ double Gen_r3_p_M1[10]={0, 1.50861, 2.09508, 2.05338, 1.96178, 1.97999, 2.07162, 1.8607, 1.91805, 1.10468};
+ double Gen_r3_p_M1_e[10]={0, 0.3399, 0.109449, 0.0646241, 0.0544419, 0.0667711, 0.0882695, 0.150195, 0.309844, 0.643923};*/
+ // Kt3 2
+ /*double Gen_r3_m_M0[10]={0, 1.04272, 2.05961, 1.93025, 2.08203, 2.07565, 2.29192, 1.93009, 2.68715, 1.71175};
+ double Gen_r3_m_M0_e[10]={0, 3.09343, 0.30016, 0.117634, 0.0757208, 0.0794904, 0.090823, 0.128235, 0.21308, 0.40023};
+ double Gen_r3_p_M0[10]={0, 1.83276, 1.88969, 2.03778, 2.0907, 2.11919, 2.04992, 2.02593, 1.95209, 1.68264};
+ double Gen_r3_p_M0_e[10]={0, 3.59, 0.298056, 0.118354, 0.0762849, 0.079909, 0.0912179, 0.128866, 0.213999, 0.402173};
+ double Gen_r3_m_M1[10]={0, 5.40628, 1.52822, 1.93258, 2.13338, 2.05811, 2.02963, 2.1204, 2.04906, 1.9021};
+ double Gen_r3_m_M1_e[10]={0, 4.3025, 0.350163, 0.13871, 0.0897272, 0.0938973, 0.107045, 0.151557, 0.25029, 0.446325};
+ double Gen_r3_p_M1[10]={0, 3.15883, 1.86772, 2.24914, 2.12118, 2.09175, 2.01447, 2.36802, 2.57239, 2.68729};
+ double Gen_r3_p_M1_e[10]={0, 4.6622, 0.348237, 0.140294, 0.0903003, 0.09446, 0.107692, 0.152575, 0.251939, 0.448919};*/
+
+ //double r3Sys_M1[10]={0, 0.097, 0.056, 0.063, 0.097, 0.17, 0.32, 0.66, 1.4, 3.4};
+ // muon variation
+ /*double r3_muonVar_M1[10]={0, 1.73244, 1.80921, 1.77852, 1.7192, 1.62059, 1.50122, 1.37656, 1.01344, 0.755781};
+ double r3_muonVar_M1_e[10]={0, 0.160786, 0.051075, 0.031982, 0.0293467, 0.0390241, 0.0514494, 0.0760959, 0.112944, 0.154592};
+ TH1D *r3_muonVar=(TH1D*)r3_num_Q3->Clone();
+ for(int i=0; i<10; i++){
+ r3_muonVar->SetBinContent(i+1,r3_muonVar_M1[i]);
+ r3_muonVar->SetBinError(i+1,sqrt(pow(r3_muonVar_M1_e[i],2)+pow(r3Sys_M1[i],2)));
+ }
+ r3_muonVar->SetMarkerStyle(21);
+ r3_muonVar->SetMarkerColor(2); r3_muonVar->SetLineColor(2);
+ r3_muonVar->SetFillStyle(1000); r3_muonVar->SetFillColor(kRed-10);
+ r3_muonVar->Draw("E2 same");
+ r3_num_Q3->Draw("same");
+ legend3->AddEntry(r3_num_Q3,"-2.0<N#sigma_{Pion}<2.0 (default)","p");
+ legend3->AddEntry(r3_muonVar,"-0.5<N#sigma_{Pion}<2.0","p");
+ legend3->Draw("same");*/
+
+ // muon correction
+ /*double r3_muonCorr_M1[10]={0, 1.8462, 1.84371, 1.79934, 1.77217, 1.76725, 1.79545, 1.7986, 2.11717, 2.86177};
+ double r3_muonCorr_M1_e[10]={0, 0.0838277, 0.0266269, 0.0168108, 0.0156166, 0.0211333, 0.0286836, 0.0450595, 0.075618, 0.141419};
+ TH1D *r3_muonCorr=(TH1D*)r3_num_Q3->Clone();
+ for(int i=0; i<10; i++){
+ r3_muonCorr->SetBinContent(i+1,r3_muonCorr_M1[i]);
+ r3_muonCorr->SetBinError(i+1,sqrt(pow(r3_muonCorr_M1_e[i],2) + pow(r3Sys_M1[i],2)));
+ }
+ r3_muonCorr->SetMarkerStyle(21);
+ r3_muonCorr->SetMarkerColor(2); r3_muonCorr->SetLineColor(2);
+ r3_muonCorr->SetFillStyle(1000); r3_muonCorr->SetFillColor(kRed-10);
+ r3_muonCorr->Draw("E2 same");
+ r3_num_Q3->Draw("same");
+ legend3->AddEntry(r3_num_Q3,"No Muon Correction","p");
+ legend3->AddEntry(r3_muonCorr,"Correction Applied","p");
+ legend3->Draw("same");*/
+
+ // muon correction for muon variation data
+ /*double r3_muonCorr_M1[10]={0, 1.71115, 1.8128, 1.79761, 1.76112, 1.70699, 1.63908, 1.63417, 1.49861, 1.75565};
+ double r3_muonCorr_M1_e[10]={0, 0.155228, 0.0491558, 0.0308125, 0.0283357, 0.0378035, 0.0501422, 0.0756057, 0.119169, 0.202093};
+ TH1D *r3_muonCorr=(TH1D*)r3_num_Q3->Clone();
+ for(int i=0; i<10; i++){
+ r3_muonCorr->SetBinContent(i+1,r3_muonCorr_M1[i]);
+ r3_muonCorr->SetBinError(i+1,sqrt(pow(r3_muonCorr_M1_e[i],2) + pow(r3Sys_M1[i],2)));
+ }
+ r3_muonCorr->SetMarkerStyle(21);
+ r3_muonCorr->SetMarkerColor(2); r3_muonCorr->SetLineColor(2);
+ r3_muonCorr->SetFillStyle(1000); r3_muonCorr->SetFillColor(kRed-10);
+ r3_muonCorr->Draw("E2 same");
+ r3_num_Q3->Draw("same");
+ legend3->AddEntry(r3_num_Q3,"Muon Corrected. Default PID","p");
+ legend3->AddEntry(r3_muonCorr,"Muon Corrected. Varied PID","p");
+ legend3->Draw("same");*/
+
+ /*for(int i=1; i<=10; i++){
+ cout<<r3_num_Q3->GetBinContent(i)<<", ";
+ }
+ cout<<endl;
+ for(int i=1; i<=10; i++){
+ cout<<r3_num_Q3->GetBinError(i)<<", ";
+ }*/
+
/*
- r3_num_Q3->Fit(QuarticFit,"IME","",0,0.1);
+ TH1D *Merged_SanityCheck=(TH1D*)r3_num_Q3->Clone();
+ for(int i=1; i<=10; i++){
+ double value = (Gen_r3_m_M0[i-1]+Gen_r3_p_M0[i-1])/2.;
+ double value_e = sqrt(pow(Gen_r3_m_M0_e[i-1],2)+pow(Gen_r3_p_M0_e[i-1],2))/2.;
+ //double value = (Gen_r3_m_M0[i-1]+Gen_r3_p_M0[i-1]+Gen_r3_m_M1[i-1]+Gen_r3_p_M1[i-1])/4.;
+ //double value_e = sqrt(pow(Gen_r3_m_M0_e[i-1],2)+pow(Gen_r3_p_M0_e[i-1],2)+pow(Gen_r3_m_M1_e[i-1],2)+pow(Gen_r3_p_M1_e[i-1],2))/4.;
+ Merged_SanityCheck->SetBinContent(i,value);
+ Merged_SanityCheck->SetBinError(i,value_e);
+ }
+ gPad->SetTopMargin(0.02); gPad->SetLeftMargin(0.1);
+ gPad->SetGridx(0); gPad->SetGridy(0);
+ Merged_SanityCheck->GetXaxis()->SetTitleOffset(1.2);
+ Merged_SanityCheck->GetYaxis()->SetTitleOffset(1.3);
+ Merged_SanityCheck->SetMinimum(0.3); Merged_SanityCheck->SetMaximum(2.68);
+ Merged_SanityCheck->Draw();
+
+
+ //r3_num_Q3->Fit(QuarticFit,"IME","",0,0.1);
+ Merged_SanityCheck->Fit(QuarticFit,"IME","",0,0.1);
gPad->Update();
- TPaveStats *Quartic_stats=(TPaveStats*)r3_num_Q3->FindObject("stats");
- Quartic_stats->SetX1NDC(0.2);
- Quartic_stats->SetX2NDC(0.5);
- Quartic_stats->SetY1NDC(.8);
- Quartic_stats->SetY2NDC(.9);
+ //TPaveStats *Quartic_stats=(TPaveStats*)r3_num_Q3->FindObject("stats");
+ TPaveStats *Quartic_stats=(TPaveStats*)Merged_SanityCheck->FindObject("stats");
+ Quartic_stats->SetX1NDC(0.15);
+ Quartic_stats->SetX2NDC(0.52);
+ Quartic_stats->SetY1NDC(.2);
+ Quartic_stats->SetY2NDC(.3);
- TH1D *r3_clone=(TH1D*)r3_num_Q3->Clone();
+ //TH1D *r3_clone=(TH1D*)r3_num_Q3->Clone();
+ TH1D *r3_clone=(TH1D*)Merged_SanityCheck->Clone();
r3_clone->Fit(QuadraticFit,"IME","",0,0.1);
gPad->Update();
TPaveStats *Quadratic_stats=(TPaveStats*)r3_clone->FindObject("stats");
- Quadratic_stats->SetX1NDC(0.55);
- Quadratic_stats->SetX2NDC(0.85);
- Quadratic_stats->SetY1NDC(.8);
- Quadratic_stats->SetY2NDC(.9);
+ Quadratic_stats->SetX1NDC(0.54);
+ Quadratic_stats->SetX2NDC(0.91);
+ Quadratic_stats->SetY1NDC(.2);
+ Quadratic_stats->SetY2NDC(.3);
QuarticFit->Draw("same");
QuadraticFit->Draw("same");
Quartic_stats->Draw("same");
Quadratic_stats->Draw("same");
+
+ TF1 *ChaoticLimit = new TF1("ChaoticLimit","2.0",0,1);
+ ChaoticLimit->SetLineStyle(3);
+ ChaoticLimit->Draw("same");
+
+ TLatex *Specif_SanityCheck=new TLatex(0.2,0.35,"HIJING With Simulated QS+FSI Weights");
+ Specif_SanityCheck->SetNDC();
+ Specif_SanityCheck->SetTextFont(42);
+ Specif_SanityCheck->SetTextSize(0.04);
+ Specif_SanityCheck->Draw();
+ //TLatex *Specif_Kt3=new TLatex(0.2,0.45,"0.16<K_{t,3}<0.3 GeV/c");
+ TLatex *Specif_Kt3=new TLatex(0.2,0.45,"0.3<K_{t,3}<1.0 GeV/c");
+ Specif_Kt3->SetNDC();
+ Specif_Kt3->SetTextFont(42);
+ Specif_Kt3->SetTextSize(0.04);
+ Specif_Kt3->Draw();
+
+ legend3->AddEntry(QuarticFit,"Quartic fit","l");
+ legend3->AddEntry(QuadraticFit,"Quadratic fit","l");
+ legend3->Draw("same");
+ //DrawALICELogo(kFALSE, .72, .4, .87, .55);
*/
/*
//////////
C2_os_Ksys->Write("C2_os_Ksys");
fitC2ss->Write("fitC2ss");
fitC2os->Write("fitC2os");
+ fitC2ss_h->Write("fitC2ss_h");
+ fitC2os_h->Write("fitC2os_h");
if(IncludeEWfromTherm) {
fitC2ssTherm->Write("fitC2ssTherm");
fitC2osTherm->Write("fitC2osTherm");
}
-void ReadCoulCorrections(int ST, int RVal, int bVal, int kt){
+void ReadCoulCorrections(int RVal, int bVal, int kt){
///////////////////////
- TString *fname;
- if(FileSetting!=6) fname = new TString("KFile.root");
- if(FileSetting==6) fname = new TString("KFile_Gauss.root");
-
- TFile *File=new TFile(fname->Data(),"READ");
- if(ST==0){// Gaussian
- if(RVal < 3 || RVal > 10) cout<<"Coulomb Correlation Gaussian radius outside of range!!!!!!!!!!!!!!!!"<<endl;
- TH2D *tempG_ss = (TH2D*)File->Get("K2ssG");
- CoulCorr2SS = (TH1D*)tempG_ss->ProjectionY("CoulCorr2SS",RVal-2, RVal-2);
- TH2D *tempG_os = (TH2D*)File->Get("K2osG");
- CoulCorr2OS = (TH1D*)tempG_os->ProjectionY("CoulCorr2OS",RVal-2, RVal-2);
+ TString *fname2;
+ if(FileSetting==6) fname2 = new TString("KFile_GaussR8to5.root");
+ else if(FileSetting==7) fname2 = new TString("KFile_GaussR11to6.root");
+ else{
+ if(SourceType==0) fname2 = new TString("KFile.root");
+ if(SourceType==1) fname2 = new TString("KFile_50fmCut.root");
+ if(SourceType==2) fname2 = new TString("KFile_GaussR11to6.root");
}
- if(ST==1){//Therminator
- if(bVal!=2 && bVal!=3 && bVal!=5 && bVal!=7 && bVal!=8 && bVal!=9) cout<<"Therminator bVal not acceptable in 2-particle Coulomb read"<<endl;
-
- if(kt==10){// kt integrated
- TH2D *tempT_ss = (TH2D*)File->Get("K2ssT");
- TH2D *tempT_os = (TH2D*)File->Get("K2osT");
- CoulCorr2SS = (TH1D*)tempT_ss->ProjectionY("CoulCorr2SS",bBin, bBin);
- CoulCorr2OS = (TH1D*)tempT_os->ProjectionY("CoulCorr2OS",bBin, bBin);
- }else{
- if(kt < 1 || kt > 6) cout<<"kt bin out of range in 2-particle Coulomb read"<<endl;
- TH3D *tempT3_ss = (TH3D*)File->Get("K2ssT_kt");
- TH3D *tempT3_os = (TH3D*)File->Get("K2osT_kt");
- CoulCorr2SS = (TH1D*)tempT3_ss->ProjectionZ("CoulCorr2SS",bBin, bBin, kt,kt);
- CoulCorr2OS = (TH1D*)tempT3_os->ProjectionZ("CoulCorr2OS",bBin, bBin, kt,kt);
+
+ TFile *File=new TFile(fname2->Data(),"READ");
+ if(bVal!=2 && bVal!=3 && bVal!=5 && bVal!=7 && bVal!=8 && bVal!=9) cout<<"Therminator bVal not acceptable in 2-particle Coulomb read"<<endl;
+
+ if(kt==10){// kt integrated
+ TH2D *tempT_ss = (TH2D*)File->Get("K2ssT");
+ TH2D *tempT_os = (TH2D*)File->Get("K2osT");
+ CoulCorr2SS = (TH1D*)tempT_ss->ProjectionY("CoulCorr2SS",bBin, bBin);
+ CoulCorr2OS = (TH1D*)tempT_os->ProjectionY("CoulCorr2OS",bBin, bBin);
+ }else{
+ if(kt < 1 || kt > 6) cout<<"kt bin out of range in 2-particle Coulomb read"<<endl;
+ TH3D *tempT3_ss = (TH3D*)File->Get("K2ssT_kt");
+ TH3D *tempT3_os = (TH3D*)File->Get("K2osT_kt");
+ CoulCorr2SS = (TH1D*)tempT3_ss->ProjectionZ("CoulCorr2SS",bBin, bBin, kt,kt);
+ CoulCorr2OS = (TH1D*)tempT3_os->ProjectionZ("CoulCorr2OS",bBin, bBin, kt,kt);
+ }
+
+ if(IncludeStrongSC){// include strong FSI for pi+pi+ as per Lednicky code with R=7 (ratio of Coul+Strong to Coul)
+ for(int ii=1; ii<=CoulCorr2SS->GetNbinsX(); ii++){
+ double newValue = CoulCorr2SS->GetBinContent(ii) * StrongSC->Eval(CoulCorr2SS->GetBinCenter(ii)*1000/2.);// k* not qinv
+ CoulCorr2SS->SetBinContent(ii, newValue);
}
}
+ // K factor shift for testing
+ for(int ii=1; ii<=CoulCorr2SS->GetNbinsX(); ii++){
+ double newValueSS = (CoulCorr2SS->GetBinContent(ii)-1)*KShift+1;// k* not qinv
+ CoulCorr2SS->SetBinContent(ii, newValueSS);
+ double newValueOS = (CoulCorr2OS->GetBinContent(ii)-1)*KShift+1;// k* not qinv
+ CoulCorr2OS->SetBinContent(ii, newValueOS);
+ }
+
CoulCorr2SS->SetDirectory(0);
CoulCorr2OS->SetDirectory(0);
File->Close();
indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
indexH = indexL+1;
- if(indexH >= Nlines) return 1.0;
+ if(indexH >= CoulCorr2SS->GetNbinsX()) return 1.0;
if(chargeproduct==1){
slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
return slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
}
- /*
- int Qbin=CoulCorr2SS->GetXaxis()->FindBin(Q2);
- if(chargeproduct==1) return CoulCorr2SS->GetBinContent(Qbin);
- else return CoulCorr2OS->GetBinContent(Qbin);
- */
}
//----------------------------------------------------------------------
for(int i=1; i<BINRANGE_C2global; i++){
qinvSS = AvgQinvSS[i];
-
+ //if(qinvSS>0.08) continue;
+
if(!GofP) Dp=fabs(par[2])/(1-fabs(par[2]));// p independent
else Dp = fabs(par[2])/(1-fabs(par[2]));
//else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
CSS *= par[1]*K2SS[i];
+ if(MuonCorrection) CSS /= C2muonCorrection_SC->GetBinContent(i+1);
CSS += 1-par[1];
CSS *= par[0];
//
COS = 1;
if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
COS *= par[1]*K2OS[i];
+ if(MuonCorrection) COS /= C2muonCorrection_MC->GetBinContent(i+1);
COS += 1-par[1];
COS *= par[9];// different Norm
//
else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
//else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
CSS *= par[1]*K2SS[qbin];
+ if(MuonCorrection) CSS /= C2muonCorrection_SC->GetBinContent(qbin+1);
CSS += 1-par[1];
CSS *= par[0];
double COS = 1;
if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
COS *= par[1]*K2OS[qbin];
+ if(MuonCorrection) COS /= C2muonCorrection_MC->GetBinContent(qbin+1);
COS += 1-par[1];
COS *= par[9];
return COS;
void ReadCoulCorrections_Omega0(){
// read in 3d 3-particle coulomb correlations = K3
- TFile *coulfile;
- if(FileSetting!=6) coulfile = new TFile("KFile.root","READ");
- if(FileSetting==6) coulfile = new TFile("KFile_Gauss.root","READ");
-
+ TString *fname;
+ if(FileSetting==6) fname = new TString("KFile_GaussR8to5.root");
+ else if(FileSetting==7) fname = new TString("KFile_GaussR11to6.root");
+ else{
+ if(SourceType==0) fname = new TString("KFile.root");
+ if(SourceType==1) fname = new TString("KFile_50fmCut.root");
+ if(SourceType==2) fname = new TString("KFile_GaussR11to6.root");
+ }
+ TFile *coulfile=new TFile(fname->Data(),"READ");
+
TString *name=new TString("K3ss_");
*name += bBin-1;
CoulCorrOmega0SS = (TH3D*)coulfile->Get(name->Data());
}
}
-double CoulCorrOmega0(bool SC, double Q_12, double Q_13, double Q_23){
+double CoulCorrOmega0(bool SC, double Q_12, double Q_13, double Q_23){// 12 is same-charge pair
int Q12bin = CoulCorrOmega0SS->GetXaxis()->FindBin(Q_12);
int Q13bin = CoulCorrOmega0SS->GetZaxis()->FindBin(Q_13);
int Q23bin = CoulCorrOmega0SS->GetYaxis()->FindBin(Q_23);
int index23L = int(fabs(Q_23 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
int index23H = index23L+1;
- if(SC){
- if(CoulCorrOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
-
- if(Tricubic){
- //////////////////////////
- // Tricubic Interpolation
- double arr[4][4][4]={{{0}}};
- for(int x=0; x<4; x++){
- for(int y=0; y<4; y++){
- for(int z=0; z<4; z++){
- arr[x][y][z] = CoulCorrOmega0SS->GetBinContent((index12L)+x, (index23L)+y, (index13L)+z);
- }
+
+ if(SC) {if(CoulCorrOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic same-charge Omega0 bin"<<endl;}}
+ else {if(CoulCorrOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic mixed-charge Omega0 bin"<<endl;}}
+
+ if(Tricubic){
+ //////////////////////////
+ // Tricubic Interpolation
+ double arr[4][4][4]={{{0}}};
+ for(int x=0; x<4; x++){
+ for(int y=0; y<4; y++){
+ for(int z=0; z<4; z++){
+ if(SC) arr[x][y][z] = CoulCorrOmega0SS->GetBinContent((index12L)+x, (index23L)+y, (index13L)+z);
+ else arr[x][y][z] = CoulCorrOmega0OS->GetBinContent((Q12bin)+x, (Q23bin)+y, (Q13bin)+z);
}
}
- return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
- }else{
- ///////////////////////////
- // Trilinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
- //
- double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
- xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
- double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
- yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
- double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
- zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
-
+ }
+ return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
+ }else{
+ ///////////////////////////
+ // Trilinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
+ //
+ double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
+ xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
+ double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
+ yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
+ double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
+ zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
+
+
+ double c00=0, c10=0, c01=0, c11=0;
+ if(SC){
// interpolate along x
- double c00=0, c10=0, c01=0, c11=0;
if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
- //
- double c0=0, c1=0;
- if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
- if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
-
- if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
- else {
- cout<<"Not all Omega0 bins well defined. Use GRS instead"<<endl;
- return CoulCorrGRS(kTRUE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
- }
- }
-
- }else {// mixed charge. Q12bin always designated as the same-charge pair
- if(CoulCorrOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
-
- if(Tricubic){
- //////////////////////////
- // Tricubic Interpolation
- double arr[4][4][4]={{{0}}};
- for(int x=0; x<4; x++){
- for(int y=0; y<4; y++){
- for(int z=0; z<4; z++){
- arr[x][y][z] = CoulCorrOmega0OS->GetBinContent((Q12bin)+x, (Q23bin)+y, (Q13bin)+z);
- }
- }
- }
- return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
- }else{
- ///////////////////////////
- // TriLinear Interpolation. See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
- //
- double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
- xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
- double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
- yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
- double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
- zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
- // interpolate along x
- double c00=0, c10=0, c01=0, c11=0;
+ }else {
if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
- //
- double c0=0, c1=0;
- if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
- if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
-
- if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
- else {
- cout<<"Not all Omega0 bins well defined. Use GRS instead"<<endl;
- return CoulCorrGRS(kFALSE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
- }
+ }
+ //
+ double c0=0, c1=0;
+ if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
+ if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
+
+ if(c0>0 && c1>0) {
+ double valueOm = (c0*(1-zd) + c1*zd);
+ if(SC) valueOm *= StrongSC->Eval(Q_12*1000/2.) * StrongSC->Eval(Q_23*1000/2.) * StrongSC->Eval(Q_13*1000/2.);
+ else valueOm *= StrongSC->Eval(Q_12*1000/2.);
+ return valueOm;
+ }else {
+ cout<<"Not all Omega0 bins well defined. Use GRS instead"<<endl;
+ return CoulCorrGRS(kTRUE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
}
}
-
}
double Gamov(int chargeProduct, double qinv){
MomResDev_term2_pp[i] = 1.;
}
-
TString *momresfilename = new TString("MomResFile");
momresfilename->Append("_Offline_");
if(FileSetting==3) momresfilename->Append("PID1p5");
+ else if(FileSetting==9) momresfilename->Append("FB5and7overlap");
else momresfilename->Append("11h");
momresfilename->Append(".root");// no corresponding file for 10h
arr[3] = bicubicInterpolate(p[3], y, z);
return cubicInterpolate(arr, x);
}
+float fact(float n){
+ return (n < 1.00001) ? 1 : fact(n - 1) * n;
+}
+void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2)
+{
+ // revision on July 24th or 25th 2012
+ // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
+ x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
+ // Printf("%f %f %f %f", x1, x2, y1, y2);
+
+ TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
+ myPadLogo->SetLeftMargin(0);
+ myPadLogo->SetTopMargin(0);
+ myPadLogo->SetRightMargin(0);
+ myPadLogo->SetBottomMargin(0);
+ myPadLogo->Draw();
+ myPadLogo->cd();
+ TASImage *myAliceLogo = new TASImage((prel) ? "~/Pictures/2011-Nov-24-ALICE_PRELIMARY_logo_BLACK_small_usage_design.eps" : "~/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.eps");
+ myAliceLogo->Draw();
+}
+//________________________________________________________________________________________
+void ReadMuonCorrections(int mbin){
+ TString *name = new TString("MuonCorrection_");
+ if(mbin<=1) *name += 0;
+ else if(mbin<=3) *name += 1;
+ else if(mbin<=5) *name += 2;
+ else if(mbin<=7) *name += 3;
+ else if(mbin<=9) *name += 3;
+ else *name += 4;
+ name->Append(".root");
+ TFile *muonfile=new TFile(name->Data(),"READ");
+ C2muonCorrection_SC = (TH1D*)muonfile->Get("C2muonCorrection_SC");
+ C2muonCorrection_MC = (TH1D*)muonfile->Get("C2muonCorrection_MC");
+ WeightmuonCorrection = (TH1D*)muonfile->Get("WeightmuonCorrection");
+ if(SameCharge_def) C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_SC");
+ else C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_MC");
+ //
+ C2muonCorrection_SC->SetDirectory(0);
+ C2muonCorrection_MC->SetDirectory(0);
+ WeightmuonCorrection->SetDirectory(0);
+ C3muonCorrection->SetDirectory(0);
+ //
+ muonfile->Close();
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff