--- /dev/null
+//gROOT->Reset();
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+#include "RooFit.h"
+#include "RooVoigtian.h"
+#include "RooAbsReal.h"
+#include "RooRealVar.h"
+#include "RooComplex.h"
+#include "RooMath.h"
+
+#include "TVector2.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TMath.h"
+#include "TText.h"
+#include "TRandom3.h"
+#include "TArray.h"
+#include "TLegend.h"
+#include "TStyle.h"
+#include "TMinuit.h"
+
+#include "TCanvas.h"
+#include "TLatex.h"
+#include "TImage.h"
+#include "TGaxis.h"
+#include "TGraphErrors.h"
+#include "TPaveStats.h"
+
+using namespace std;
+
+
+double myLevyPtFunc(Double_t *x, Double_t *par);
+double myMtExpFunc(Double_t *x, Double_t *par);
+double myPtExpFunc(Double_t *x, Double_t *par);
+double myBoltzFunc(Double_t *x, Double_t *par);
+double myPowerLawPtFunc(Double_t *x, Double_t *par);
+double myPowerLawFunc(Double_t *x, Double_t *par);
+double myBoltsBWFunc(Double_t *x, Double_t *par);
+double IntegrandBG(const double *x, const double *p);
+
+void Plot_XiStarResults(){
+
+ TGaxis::SetMaxDigits(3);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(01111);
+
+
+ TFile *filesCutVar[13];
+ TH1D *Spectra_CutVar[13];
+ TF1 *LevyFits_CutVar[13];
+ TH1D *Yields_CutVar[13];
+ TH1D *Eff_CutVar[13];
+ for(int i=0; i<13; i++){
+ TString *name=new TString("CutVariation_");
+ *name += i;
+ name->Append(".root");
+ filesCutVar[i]=new TFile(name->Data());
+ Spectra_CutVar[i] = (TH1D*)filesCutVar[i]->Get("h_Spectrum");
+ Spectra_CutVar[i]->SetDirectory(0);
+ LevyFits_CutVar[i] = (TF1*)filesCutVar[i]->Get("myLevy");
+ Yields_CutVar[i] = (TH1D*)filesCutVar[i]->Get("RawYields");
+ Eff_CutVar[i] = (TH1D*)filesCutVar[i]->Get("Efficiency");
+ }
+ //
+ TFile *filesFitVar[4];
+ TH1D *Spectra_FitVar[4];
+ TF1 *LevyFits_FitVar[4];
+ TH1D *Yields_FitVar[4];
+ for(int i=0; i<4; i++){
+ TString *name=new TString("FitVariation_");
+ *name += i+1;
+ name->Append(".root");
+ filesFitVar[i]=new TFile(name->Data());
+ Spectra_FitVar[i] = (TH1D*)filesFitVar[i]->Get("h_Spectrum");
+ Spectra_FitVar[i]->SetDirectory(0);
+ LevyFits_FitVar[i] = (TF1*)filesFitVar[i]->Get("myLevy");
+ Yields_FitVar[i] = (TH1D*)filesFitVar[i]->Get("RawYields");
+ }
+
+ double FVavg=0;
+ for(int bin=2; bin<=9; bin++){
+ // FV 1 has largest spectrum point variation
+ FVavg += fabs(Spectra_FitVar[0]->GetBinContent(bin)-Spectra_FitVar[1]->GetBinContent(bin))/Spectra_FitVar[0]->GetBinContent(bin);
+ }
+ cout<<"Average point-by-point spectrum variation due to fitting: "<<FVavg/8<<endl;
+
+ double TrackVavg=0;
+ for(int bin=2; bin<=9; bin++){
+ // CutVar 5 has largest spectrum point variation in the "tracking systematic" category
+ TrackVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[5]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin);
+ }
+ cout<<"Average point-by-point spectrum variation due to Tracking: "<<TrackVavg/8<<endl;
+
+ double TopologVavg=0;
+ for(int bin=2; bin<=9; bin++){
+ // CutVar 8 has largest spectrum point variation in the "topological systematic" category
+ TopologVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[8]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin);
+ }
+ cout<<"Average point-by-point spectrum variation due to Topological Cuts: "<<TopologVavg/8<<endl;
+
+ ////////////////////////////////////
+ TH1D *FinalSpectrum = (TH1D*)Spectra_CutVar[0]->Clone();
+ FinalSpectrum->SetMarkerSize(1.);
+ double StatError[8]={0};
+ for(int ptbin=1; ptbin<=8; ptbin++) StatError[ptbin-1] = FinalSpectrum->GetBinError(ptbin+1);// add 1 since first bin is a dummy
+
+
+ TH1D *PureSysFromCutVar[12];
+ double MeanSysCutVar[12]={0};
+ double sigmaSysCutVar[12]={0};
+ for(int i=0; i<12; i++) {
+ TString *name=new TString("PureSysFromCutVar");
+ *name += i+1;
+ PureSysFromCutVar[i] = new TH1D(name->Data(),"",8,0.5,8.5);
+ PureSysFromCutVar[i]->GetXaxis()->SetTitle("pt bin");
+ PureSysFromCutVar[i]->GetYaxis()->SetTitle("|y_{s}-y_{v}|/y_{s}");
+ }
+
+ //////////////////////////////////////
+ // Sys error assignment for pt bins to be used for mean value fits
+
+ int CutVarLow=1, CutVarHigh=13;
+ for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13
+ for(int ptbin=2; ptbin<=9; ptbin++){
+ double SysVar = fabs(FinalSpectrum->GetBinContent(ptbin) - Spectra_CutVar[i]->GetBinContent(ptbin))/FinalSpectrum->GetBinContent(ptbin);
+ double sigma = sqrt(fabs(pow(FinalSpectrum->GetBinError(ptbin),2) - pow(Spectra_CutVar[i]->GetBinError(ptbin),2)));// stat part
+ PureSysFromCutVar[i-1]->SetBinContent(ptbin-1, SysVar);
+ PureSysFromCutVar[i-1]->SetBinError(ptbin-1, sigma/FinalSpectrum->GetBinContent(ptbin));
+ }
+ PureSysFromCutVar[i-1]->Fit("pol0","IMEQ","",0,9);
+ MeanSysCutVar[i-1] = pol0->GetParameter(0);
+ sigmaSysCutVar[i-1] = pol0->GetParError(0);
+ cout<<"Cut Var "<<i<<" Mean Sys = "<<MeanSysCutVar[i-1]<<" += "<<pol0->GetParError(0)<<endl;
+ //cout<<"Chi2/NDF = "<<pol0->GetChisquare()/pol0->GetNDF()<<endl;
+ //TPaveStats *sb2 = (TPaveStats*)(PureSysFromCutVar[i-1]->GetListOfFunctions()->FindObject("stats"));
+ //sb2->SetX1NDC(.2);
+ }
+
+
+ double SysError[9]={0};
+ // Cut Variations
+ //int CutVarN=12;
+ for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13
+ for(int ptbin=1; ptbin<=9; ptbin++){
+ if(MeanSysCutVar[i-1] > 0 && MeanSysCutVar[i-1]>sigmaSysCutVar[i-1]) {
+ SysError[ptbin-1] += pow(MeanSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2) - pow(sigmaSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2);
+ }
+ }
+ }
+ // Fit Variations
+ for(int i=2; i<3; i++){// now only use i=2 (max deviation). was 0 to 4
+ for(int ptbin=1; ptbin<=9; ptbin++){
+ SysError[ptbin-1] += pow(FinalSpectrum->GetBinContent(ptbin) - Spectra_FitVar[i]->GetBinContent(ptbin),2);
+ }
+ }
+
+ // Add Stat errors and set Final Error
+ for(int ptbin=1; ptbin<=9; ptbin++){
+ double TotError = pow(FinalSpectrum->GetBinError(ptbin),2);// Stat
+ TotError += SysError[ptbin-1];// Add on Sys
+ FinalSpectrum->SetBinError(ptbin, sqrt(TotError));
+ cout<<"Spectrum ptbin:"<<ptbin-1<<" "<<FinalSpectrum->GetBinContent(ptbin)<<" +- "<<FinalSpectrum->GetBinError(ptbin)<<endl;
+ }
+
+ for(int ptbin=1; ptbin<=9; ptbin++) SysError[ptbin-1] = sqrt(SysError[ptbin-1]);
+
+ // print efficiencies
+ for(int ptbin=1; ptbin<=8; ptbin++){
+ cout<<"Efficiency ptbin:"<<ptbin<<" "<<Eff_CutVar[0]->GetBinContent(ptbin)<<" +- "<<Eff_CutVar[0]->GetBinError(ptbin)<<endl;
+ }
+
+ TCanvas *can = new TCanvas("can","can",13,34,700,500);
+ can->Range(-1.25,-0.2625,11.25,2.3625);
+ can->SetFillColor(10);
+ can->SetBorderMode(0);
+ can->SetBorderSize(2);
+ //can->SetGridx();
+ //can->SetGridy();
+ can->SetFrameFillColor(0);
+ can->SetFrameBorderMode(0);
+ can->SetFrameBorderMode(0);
+
+ TLegend *legend = new TLegend(.1,.7,.35,.9,NULL,"brNDC");
+ legend->SetBorderSize(1);
+ legend->SetTextSize(.04);// small .03; large .036
+ //legend->SetLineColor(0);
+ legend->SetFillColor(0);
+ TLegend *legend2 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
+ legend2->SetBorderSize(1);
+ legend2->SetTextSize(.04);// small .03; large .036
+ //legend2->SetLineColor(0);
+ legend2->SetFillColor(0);
+
+
+ double minFitpoint = 0.8; double maxFitpoint = 5.6;
+
+
+
+ TF1 *myLevy=new TF1("myLevy",myLevyPtFunc,0,8.5,3);
+ myLevy->SetParName(0,"dN/dy");
+ myLevy->SetParName(1,"C");
+ myLevy->SetParName(2,"n");
+ myLevy->SetParameter(0,.008);
+ myLevy->SetParameter(1,.3);
+ myLevy->SetParameter(2,15);
+ myLevy->SetParLimits(0,.001,.2);
+ myLevy->SetParLimits(1,.1,1);
+ myLevy->SetParLimits(2,1,500);
+ myLevy->SetLineColor(1);
+ TF1 *myPtExp = new TF1("myPtExp",myPtExpFunc,0,8.5,2);
+ myPtExp->SetParName(0,"dN/dy");
+ myPtExp->SetParName(1,"C");
+ myPtExp->SetParameter(0,.003);
+ myPtExp->SetParameter(1,.3);
+ myPtExp->SetParLimits(0,0.0001,0.1);
+ myPtExp->SetParLimits(1,.1,1);
+ myPtExp->SetLineColor(2);
+ TF1 *myMtExp = new TF1("myMtExp",myMtExpFunc,0,8.5,2);
+ myMtExp->SetParName(0,"dN/dy");
+ myMtExp->SetParName(1,"C");
+ myMtExp->SetParameter(0,.003);
+ myMtExp->SetParameter(1,.3);
+ myMtExp->SetParLimits(0,0.0001,0.1);
+ myMtExp->SetParLimits(1,.1,1);
+ myMtExp->SetLineColor(4);
+ TF1 *myBoltz = new TF1("myBoltz",myBoltzFunc,0,8.5,2);
+ myBoltz->SetParName(0,"dN/dy");
+ myBoltz->SetParName(1,"C");
+ myBoltz->SetParameter(0,.003);
+ myBoltz->SetParameter(1,.3);
+ myBoltz->SetParLimits(0,0.0001,0.1);
+ myBoltz->SetParLimits(1,.1,1);
+ myBoltz->SetLineColor(6);
+ TF1 *myPowerLawPt=new TF1("myPowerLawPt",myPowerLawPtFunc,0,8.5,3);
+ myPowerLawPt->SetParName(0,"dN/dy");
+ myPowerLawPt->SetParName(1,"pt0");
+ myPowerLawPt->SetParName(2,"n");
+ myPowerLawPt->SetParameter(0,.008);
+ myPowerLawPt->SetParameter(1,50);
+ myPowerLawPt->SetParameter(2,50);
+ myPowerLawPt->SetParLimits(0,.001,.01);
+ myPowerLawPt->SetParLimits(1,1,500);
+ myPowerLawPt->SetParLimits(2,1,500);
+ myPowerLawPt->SetLineColor(7);
+ TF1 *myPowerLaw=new TF1("myPowerLaw",myPowerLawFunc,0,8.5,3);
+ myPowerLaw->SetParName(0,"dN/dy");
+ myPowerLaw->SetParName(1,"pt0");
+ myPowerLaw->SetParName(2,"n");
+ myPowerLaw->SetParameter(0,.008);
+ myPowerLaw->SetParameter(1,50);
+ myPowerLaw->SetParameter(2,50);
+ myPowerLaw->SetParLimits(0,.001,.01);
+ myPowerLaw->SetParLimits(1,1,500);
+ myPowerLaw->SetParLimits(2,1,500);
+ myPowerLaw->SetLineColor(8);
+ TF1 *myBGBlastWave=new TF1("myBGBlastWave",myBoltsBWFunc,0,8.5,5);
+ myBGBlastWave->SetParName(0,"mass");
+ myBGBlastWave->SetParName(1,"beta");
+ myBGBlastWave->SetParName(2,"C");
+ myBGBlastWave->SetParName(3,"n");
+ myBGBlastWave->SetParName(4,"Norm");
+ myBGBlastWave->FixParameter(0,1.5318);
+ myBGBlastWave->SetParameter(1,.3);
+ myBGBlastWave->SetParameter(2,.3);
+ myBGBlastWave->SetParameter(3,5);
+ myBGBlastWave->SetParameter(4,.3);
+ myBGBlastWave->SetParLimits(1,.001,.9);
+ myBGBlastWave->SetParLimits(2,.1,200);
+ myBGBlastWave->SetParLimits(3,.1,200);
+ myBGBlastWave->SetLineColor(9);
+
+
+ FinalSpectrum->Draw();
+ legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2","p");
+ //legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2: Current Results","p");
+ FinalSpectrum->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
+ FinalSpectrum->Fit(myPtExp,"IMEN","",minFitpoint,maxFitpoint);
+ FinalSpectrum->Fit(myMtExp,"IMEN","",minFitpoint,maxFitpoint);
+ FinalSpectrum->Fit(myBoltz,"IMEN","",minFitpoint,maxFitpoint);
+ FinalSpectrum->Fit(myPowerLaw,"IMEN","",minFitpoint,maxFitpoint);
+ FinalSpectrum->Fit(myPowerLawPt,"IMEN","",minFitpoint,maxFitpoint);
+ //FinalSpectrum->Fit(myBGBlastWave,"IME","",minFitpoint,maxFitpoint);
+ myLevy->Draw("same");
+ /*myPtExp->Draw("same");
+ myMtExp->Draw("same");
+ myBoltz->Draw("same");
+ myPowerLaw->Draw("same");
+ myPowerLawPt->Draw("same");*/
+ //myBGBlastWave->Draw("same");
+ //
+ legend->AddEntry(myLevy,"Levy","l");
+ /*legend->AddEntry(myPtExp,"Pt-exp","l");
+ legend->AddEntry(myMtExp,"Mt-exp","l");
+ legend->AddEntry(myBoltz,"Pt-Boltzmann","l");
+ legend->AddEntry(myPowerLaw,"PowerLaw","l");
+ legend->AddEntry(myPowerLawPt,"Pt-PowerLaw","l");
+ */
+ double AvgLowPtExtrap = myLevy->Integral(0,0.8) + myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8) + myMtExp->Integral(0,0.8);
+ AvgLowPtExtrap /= 4.;
+ cout<<"Levy Extrap % = "<<myLevy->Integral(0,0.8)/myLevy->Integral(0,10.)<<" Avg Extrap % = "<<AvgLowPtExtrap/myLevy->Integral(0,10.)<<endl;
+ //double LowPtExtrapUncertaintyPercentage = 0.17;// Levy underestimation of Pythia
+ //LowPtExtrapUncertaintyPercentage -= (AvgLowPtExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.);
+ double High_lowptExtrap = (myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8))/2.;
+ double Low_lowptExtrap = myMtExp->Integral(0,0.8);
+ double High_LowPtExtrapUncertaintyPercentage = (High_lowptExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.);
+ double Low_LowPtExtrapUncertaintyPercentage = (myLevy->Integral(0,0.8) - Low_lowptExtrap)/myLevy->Integral(0,10.);
+ cout<<"high low-pt uncertainty = "<<High_LowPtExtrapUncertaintyPercentage<<endl;
+ cout<<"low low-pt uncertainty = "<<Low_LowPtExtrapUncertaintyPercentage<<endl;
+
+
+ //PureSysFromCutVar->Draw();
+
+
+ double dNdY = myLevy->Integral(0,10);
+ double dNdY_exp = myMtExp->Integral(0,10);
+ double dNdY_Stat = LevyFits_CutVar[0]->GetParError(0);
+ double dNdY_Covered = myLevy->Integral(0.8,5.6);
+ double dNdY_Covered_exp = myMtExp->Integral(0.8,5.6);
+ double dNdY_DataPoints = 0, dNdY_DataPoints_e=0;
+ for(int ii=0; ii<8; ii++){// to explicitly integrate dNdY without fit (Christina's suggestion)
+ dNdY_DataPoints += Spectra_CutVar[0]->GetBinContent(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2);
+ dNdY_DataPoints_e += pow(Spectra_CutVar[0]->GetBinError(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2),2);
+ }
+ dNdY_DataPoints_e = sqrt(dNdY_DataPoints_e);
+ double dNdYSysPlus = 0;
+ double dNdYSysMinus = 0;
+ double CSys = 0;
+ double nSys = 0;
+ //////////////////////////////////////
+ // Sys errors of Levy parameters
+ double sigma_dNdY_cv[12]={0};
+ double sigma_C_cv[12]={0};
+ double sigma_n_cv[12]={0};
+ double SysVar_dNdY_cv[12]={0};
+ double SysVar_C_cv[12]={0};
+ double SysVar_n_cv[12]={0};
+
+ for(int cv=1; cv<13; cv++){
+ if(MeanSysCutVar[cv-1] > 0 && MeanSysCutVar[cv-1] > sigmaSysCutVar[cv-1]) {// significant deviation
+ sigma_dNdY_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_CutVar[cv]->GetParError(0),2)));
+ sigma_C_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_CutVar[cv]->GetParError(1),2)));
+ sigma_n_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_CutVar[cv]->GetParError(2),2)));
+ SysVar_dNdY_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_CutVar[cv]->GetParameter(0));
+ SysVar_C_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_CutVar[cv]->GetParameter(1));
+ SysVar_n_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_CutVar[cv]->GetParameter(2));
+
+ if(SysVar_dNdY_cv[cv-1] > sigma_dNdY_cv[cv-1]){
+ dNdYSysPlus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2);
+ dNdYSysMinus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2);
+ }
+ if(SysVar_C_cv[cv-1] > sigma_C_cv[cv-1]){
+ CSys += pow(SysVar_C_cv[cv-1],2)-pow(sigma_C_cv[cv-1],2);
+ }
+ if(SysVar_n_cv[cv-1] > sigma_n_cv[cv-1]){
+ nSys += pow(SysVar_n_cv[cv-1],2)-pow(sigma_n_cv[cv-1],2);
+ }
+
+ }
+
+ }
+
+
+
+ double sigma_dNdY_fv[4]={0};
+ double sigma_C_fv[4]={0};
+ double sigma_n_fv[4]={0};
+ double SysVar_dNdY_fv[4]={0};
+ double SysVar_C_fv[4]={0};
+ double SysVar_n_fv[4]={0};
+ for(int fv=0; fv<4; fv++){
+ sigma_dNdY_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_FitVar[fv]->GetParError(0),2)));
+ SysVar_dNdY_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2));
+ sigma_C_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_FitVar[fv]->GetParError(1),2)));
+ SysVar_C_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2));
+ sigma_n_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_FitVar[fv]->GetParError(2),2)));
+ SysVar_n_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2));
+
+ if(fv==2){//max case for dN/dy
+ //dNdYSysPlus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2);
+ //dNdYSysMinus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2);
+ dNdYSysPlus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2);
+ dNdYSysMinus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2);
+ }
+ //if(fv==2) CSys += pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2) - pow(LevyFits_CutVar[0]->GetParError(1)-LevyFits_FitVar[fv]->GetParError(1),2);
+ //if(fv==1) nSys += pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2) - pow(LevyFits_CutVar[0]->GetParError(2)-LevyFits_FitVar[fv]->GetParError(2),2);
+ if(fv==2) CSys += pow(SysVar_C_fv[fv],2) - pow(sigma_C_fv[fv],2);
+ if(fv==1) nSys += pow(SysVar_n_fv[fv],2) - pow(sigma_n_fv[fv],2);
+ }
+
+
+ double dNdYSys_forRatio = dNdYSysPlus;
+ dNdYSysPlus += pow(0.062*dNdY,2);// INEL Norm uncertainty
+ dNdYSysPlus += pow(0.04*dNdY,2);// Material Budget uncertainty
+ dNdYSysPlus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty
+
+ dNdYSysPlus += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (was LowPtExtrapUncertaintyPercentage * fabs(dNdY-dNdY_Covered))
+ dNdYSys_forRatio += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty
+ dNdYSysPlus = sqrt(dNdYSysPlus);
+ dNdYSys_forRatio = sqrt(dNdYSys_forRatio);
+ //
+ dNdYSysMinus += pow(0.03*dNdY,2);// INEL Norm uncertainty
+ dNdYSysMinus += pow(0.04*dNdY,2);// Material Budget uncertainty
+ dNdYSysMinus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty
+ dNdYSysMinus += pow(Low_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (40% from Pythia study).
+ dNdYSysMinus = sqrt(dNdYSysMinus);
+
+ CSys = sqrt(CSys);
+ nSys = sqrt(nSys);
+
+ cout<<"dN/dy = "<<dNdY<<" , Covered range = "<<dNdY_Covered<<" , DataPoint Sum = "<<dNdY_DataPoints<<" , stat = "<<dNdY_DataPoints_e<<endl;
+ cout<<"dN/dy stat = "<<dNdY_Stat<<" +sys = "<<dNdYSysPlus<<" -sys = "<<dNdYSysMinus<<endl;
+
+ cout<<"C = "<<myLevy->GetParameter(1)<<" , stat = "<<LevyFits_CutVar[0]->GetParError(1)<<" , sys = "<<CSys<<endl;
+ cout<<"n = "<<myLevy->GetParameter(2)<<" , stat = "<<LevyFits_CutVar[0]->GetParError(2)<<" , sys = "<<nSys<<endl;
+
+ // mean pt
+ double ptEdges[9]={0.8, 1.2, 1.6, 2.0, 2.4, 3.2, 4.0, 4.8, 5.6};
+ double meanpt_def=0, meanpt_Sys1=0, meanpt_Sys2=0, meanpt_Sys3=0;// Sys1 = Spectrum only, Sys2 = bin centers, Sys3 = low-pt extrapolation
+ double yieldSum_def=0, yieldSum_Sys1=0, yieldSum_Sys2=0, yieldSum_Sys3=0;
+ double meanpt_statError=0;
+ meanpt_def += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8);// low-pt extrapolation
+ yieldSum_def += myLevy->Integral(0,0.8);
+ meanpt_Sys1 += meanpt_def; meanpt_Sys2 += meanpt_def;
+ yieldSum_Sys1 += yieldSum_def; yieldSum_Sys2 += yieldSum_def;
+ meanpt_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY) * myLevy->Moment(1,0.,0.8);// low-pt extrapolation uncertainty
+ yieldSum_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY);// low-pt extrapolation uncertainty
+ meanpt_statError += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8) * myLevy->GetParError(0)/myLevy->GetParameter(0);
+
+
+ double StatErrorPtBins=0;
+ for(int bin=1; bin<=8; bin++){
+ double bin_width = fabs(ptEdges[bin-1]-ptEdges[bin]);
+ meanpt_def += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ yieldSum_def += FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ meanpt_Sys1 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * myLevy->Integral(ptEdges[bin-1],ptEdges[bin]);
+ yieldSum_Sys1 += myLevy->Integral(ptEdges[bin-1],ptEdges[bin]);
+ meanpt_Sys2 += (ptEdges[bin]+ptEdges[bin-1])/2. * FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ yieldSum_Sys2 += FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ meanpt_Sys3 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ yieldSum_Sys3 += FinalSpectrum->GetBinContent(bin+1)*bin_width;
+ //
+ StatErrorPtBins += pow(myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * Spectra_CutVar[0]->GetBinError(bin+1)*bin_width,2);
+ }
+
+ meanpt_def += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
+ yieldSum_def += myLevy->Integral(5.6,10.);
+ meanpt_Sys1 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
+ yieldSum_Sys1 += myLevy->Integral(5.6,10.);
+ meanpt_Sys2 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
+ yieldSum_Sys2 += myLevy->Integral(5.6,10.);
+ meanpt_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.))*myLevy->Moment(1,5.6,10.);// high-pt extrapolation
+ yieldSum_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.));
+ //
+ meanpt_statError += sqrt(StatErrorPtBins);
+ meanpt_statError += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.) * myLevy->GetParError(0)/myLevy->GetParameter(0);// high-pt extrapolation
+ meanpt_statError /= yieldSum_def;
+
+ meanpt_def /= yieldSum_def;
+ meanpt_Sys1 /= yieldSum_Sys1;
+ meanpt_Sys2 /= yieldSum_Sys2;
+ meanpt_Sys3 /= yieldSum_Sys3;
+
+ double meanpt_Sys = sqrt(pow(meanpt_def-meanpt_Sys1,2)+pow(meanpt_def-meanpt_Sys2,2)+pow(meanpt_def-meanpt_Sys3,2));
+
+ cout<<"<pt> = "<<meanpt_def<<" , stat = "<<meanpt_statError<<" , sys = "<<meanpt_Sys<<endl;
+
+ /*double meanpt = myLevy->Moment(1,0.,10.);// 0.,30 for full range
+ double meanpt_CoveredLevy = myLevy->Moment(1,0.8,5.6);
+ double mPt_Stat = 0;
+ double mPt_Sys = 0;
+ double mPt_StatSys = 0;
+ double base_C = myLevy->GetParameter(1);
+ double base_n = myLevy->GetParameter(2);
+
+
+ // Sys
+ myLevy->SetParameter(1, base_C + myLevy->GetParError(1));
+ mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2);
+ //
+ //
+ myLevy->SetParameter(2, base_n + myLevy->GetParError(2));
+ mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2);
+ mPt_StatSys = sqrt(mPt_StatSys);
+ //
+ //
+ // Stat
+ myLevy->SetParameter(1, base_C + LevyFits_CutVar[0]->GetParError(1));
+ mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2);
+ //
+ //
+ myLevy->SetParameter(2, base_n + LevyFits_CutVar[0]->GetParError(2));
+ mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2);
+ mPt_Stat = sqrt(mPt_Stat);
+ //
+ myLevy->SetParameter(1, base_C); myLevy->SetParameter(1, base_n);
+
+ cout<<"Levy mean pt = "<<meanpt<<" , stat = "<<mPt_Stat<<" , sys = "<<sqrt( pow(mPt_StatSys,2)-pow(mPt_Stat,2) )<<endl;
+ cout<<"Levy mean pt in covered range = "<<meanpt_CoveredLevy<<endl;
+ */
+
+ double XidNdY = 0.0079;// (Xi^+ + Xi^-)/2. Published
+ double XiStarToXi = dNdY/XidNdY;
+ double XiStarToXi_Stat = sqrt(pow(dNdY_Stat/XidNdY,2) + pow(dNdY*0.0001/pow(XidNdY,2),2));// 0.0001 is stat error of Xi^{\pm}
+ double XiStartoXi_Sys = dNdYSys_forRatio/XidNdY;
+
+
+ cout<<"Xi*/Xi = "<<XiStarToXi<<" , stat = "<<XiStarToXi_Stat<<" , sys = "<<XiStartoXi_Sys<<endl;
+ cout<<endl;
+ cout<<endl;
+
+ //Old QM2011 values for (Xi* + cc)/2
+ double Old_XiStarPoints[8]={0.001259, 0.0008586, 0.0006316, 0.0003604, 1.856e-4, 5.4601e-5, 2.0187e-5, 7.6559e-6};
+ double Old_XiStarPoints_e[8]={2.44e-4, 1.224e-4, 7.136e-5, 3.85e-5, 1.59e-5, 7.49e-6, 3.51e-6, 1.822e-6};
+ double pt_points[8]={1.0, 1.4, 1.8, 2.2, 2.8, 3.6, 4.4, 5.2};
+ double pt_points_e[8]={.2, .2, .2, .2, .4, .4, .4, .4};
+ TGraphErrors *gr_OldXiStar = new TGraphErrors(8, pt_points, Old_XiStarPoints, pt_points_e, Old_XiStarPoints_e);
+ gr_OldXiStar->SetMarkerStyle(21);
+ gr_OldXiStar->SetMarkerColor(2);
+ gr_OldXiStar->SetLineColor(2);
+ double NewOldRatio[8]={0};
+ for(int ii=0; ii<8; ii++){
+ NewOldRatio[ii] = FinalSpectrum->GetBinContent(ii+2)/Old_XiStarPoints[ii];
+ //cout<<FinalSpectrum->GetBinError(ii+2)/FinalSpectrum->GetBinContent(ii+2)<<endl;
+ }
+ TGraph *gr_NewOldRatio = new TGraph(8, pt_points, NewOldRatio);
+ gr_NewOldRatio->SetMarkerStyle(20);
+ gr_NewOldRatio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_NewOldRatio->GetYaxis()->SetTitle("New/Old");
+ gr_NewOldRatio->SetTitle("New Spectrum Divided by QM11 Spectrum");
+ //gr_NewOldRatio->Draw("AP");
+
+ //gr_OldXiStar->Draw("P");
+ //legend->AddEntry(gr_OldXiStar,"(#Xi(1530) + cc)/2: QM11 Results","p");
+ // To Fit Old points
+ /*TH1D *h_OldXiStar = (TH1D*)h_Xispectrum->Clone();
+ h_OldXiStar->Add(h_OldXiStar,-1);
+ for(int i=0; i<8; i++){
+ h_OldXiStar->SetBinContent(i+2, Old_XiStarPoints[i]);
+ h_OldXiStar->SetBinError(i+2, Old_XiStarPoints_e[i]);
+ }
+ h_OldXiStar->SetBinContent(1,+100);// first bin has no data
+ h_OldXiStar->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
+ */
+ //legend->Draw("same");
+
+ double StatPercentage[8]={};
+ double SysPercentage[8]={};
+ double TotErrPercentage[8]={};
+ for(int ii=0; ii<8; ii++){
+ StatPercentage[ii] = StatError[ii]/FinalSpectrum->GetBinContent(ii+2);
+ SysPercentage[ii] = SysError[ii+1]/FinalSpectrum->GetBinContent(ii+2);
+ TotErrPercentage[ii] = sqrt(pow(StatError[ii],2)+pow(SysError[ii+1],2))/FinalSpectrum->GetBinContent(ii+2);
+ }
+ TGraph *gr_StatPercentage = new TGraph(8, pt_points, StatPercentage);
+ gr_StatPercentage->SetMarkerStyle(25);
+ gr_StatPercentage->SetMarkerColor(4);
+ gr_StatPercentage->SetLineColor(4);
+ gr_StatPercentage->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_StatPercentage->GetYaxis()->SetTitle("Uncertainty (fraction)");
+ gr_StatPercentage->SetTitle("");
+ gr_StatPercentage->SetMinimum(0);
+ gr_StatPercentage->SetMaximum(0.2);// 0.35 to compare with Enrico
+ //gr_StatPercentage->Draw("APC");
+ TGraph *gr_SysPercentage = new TGraph(8, pt_points, SysPercentage);
+ gr_SysPercentage->SetMarkerStyle(25);
+ gr_SysPercentage->SetMarkerColor(2);
+ gr_SysPercentage->SetLineColor(2);
+ //gr_SysPercentage->Draw("PC");
+ TGraph *gr_TotErrPercentage = new TGraph(8, pt_points, TotErrPercentage);
+ gr_TotErrPercentage->SetMarkerStyle(25);
+ gr_TotErrPercentage->SetMarkerColor(1);
+ gr_TotErrPercentage->SetLineColor(1);
+ //gr_TotErrPercentage->Draw("PC");
+ //legend->AddEntry(gr_StatPercentage,"Stat. uncertainty","p");
+ //legend->AddEntry(gr_SysPercentage,"Sys. uncertainty","p");
+ //legend->AddEntry(gr_TotErrPercentage,"#sqrt{Stat^{2}+Sys^{2}}","p");
+ legend->Draw("same");
+
+ ////////////////////////////////
+ // Cut Var
+ for(int var=1; var<13; var++){// 0 to 4 for FitVar, 1 to 13 for CutVar
+ TString *Varname = new TString();
+ if(var==1) Varname->Append("Nclusters TPC");
+ if(var==2) Varname->Append("DCA PV proton");
+ if(var==3) Varname->Append("DCA PV 1st pion");
+ if(var==4) Varname->Append("DCA PV 2nd pion");
+ if(var==5) Varname->Append("DCA PV 3rd pion");
+ if(var==6) Varname->Append("DCA PV 4th pion");
+ if(var==7) Varname->Append("DCA proton-pion");
+ if(var==8) Varname->Append("DCA Lambda-pion");
+ if(var==9) Varname->Append("Decay length xy Lambda");
+ if(var==10) Varname->Append("Decay length xy Xi");
+ if(var==11) Varname->Append("Cos PointingAngle Lambda");
+ if(var==12) Varname->Append("Cos PointingAngle Xi");
+
+ cout<<Varname->Data()<<" & ";
+ cout.precision(2);
+ for(int par=0; par<3; par++) {
+ cout<<100.*(LevyFits_CutVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/LevyFits_CutVar[0]->GetParameter(par)<<"(";
+ if(par==0) cout<<100.*sigma_dNdY_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+ if(par==1) cout<<100.*sigma_C_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+ if(par==2) cout<<100.*sigma_n_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+
+ if((par+1)!=3) cout<<"& ";
+ else cout<<" \\\\\ \\hline";
+ }
+ cout<<endl;
+ }
+ cout<<endl;
+
+ //////////////////////////////
+ // Fit Var
+ for(int var=0; var<4; var++){// 0 to 4 for FitVar, 1 to 13 for CutVar
+ TString *Varname = new TString();
+ if(var==0) Varname->Append("Pure Voigtian");
+ if(var==1) Varname->Append("Pure Voigtian w/o bkg subtraction");
+ if(var==2) Varname->Append("bkg norm to the right");
+ if(var==3) Varname->Append("1/2 bin-counting range");
+
+ cout<<Varname->Data()<<" & ";
+ cout.precision(2);
+ for(int par=0; par<3; par++) {
+ cout<<100.*(LevyFits_FitVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/(LevyFits_CutVar[0]->GetParameter(par))<<"(";
+ if(par==0) cout<<100.*sigma_dNdY_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+ if(par==1) cout<<100.*sigma_C_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+ if(par==2) cout<<100.*sigma_n_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
+
+ if((par+1)!=3) cout<<"& ";
+ else cout<<" \\\\\ \\hline";
+ }
+ cout<<endl;
+ }
+ cout<<endl;
+
+ ///////////////////////////////////////////
+
+ /*
+ cout<<"Yields"<<endl;
+ // print yields
+ for(int var=1; var<=13; var++){
+ TString *Varname = new TString();
+ if(var==1) Varname->Append("Standard ");
+ else {Varname->Append("Cut Var "); *Varname += var-1;}
+
+ cout<<Varname->Data()<<" & ";
+ cout.precision(4);// 4 or 6
+ for(int ptbin=1; ptbin<=8; ptbin++){
+ //cout<<Yields_CutVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[var-1]->GetBinContent(ptbin)<<" ";
+ cout<<Eff_CutVar[var-1]->GetBinContent(ptbin)<<" ";
+ if((ptbin)!=8) cout<<"& ";
+ else cout<<" \\\\\ \\hline";
+ }
+ cout<<endl;
+ }
+ */
+ /* for(int var=1; var<5; var++){
+ TString *Varname = new TString("Fit Var ");
+ *Varname += var;
+ cout<<Varname->Data()<<" & ";
+ cout.precision(6);
+ for(int ptbin=1; ptbin<=8; ptbin++){
+ cout<<Yields_FitVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[0]->GetBinContent(ptbin)<<" ";
+ if((ptbin)!=8) cout<<"& ";
+ else cout<<" \\\\\ \\hline";
+ }
+ cout<<endl;
+ }*/
+ cout<<endl;
+
+ //TString *name=new TString("../AnaNotes/XiStar/XiStar_yield_ErrorsCutVar");
+ //*name += CutVarN;
+ /*TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSys");
+ *name += CutVarLow;
+ name->Append(".png");
+ c1->SaveAs(name->Data());*/
+ //TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSummary");
+ //*name += CutVarLow;
+ //name->Append(".png");
+ //can->SaveAs(name->Data());
+
+}
+
+
+//________________________________________________________________________
+double myLevyPtFunc(Double_t *x, Double_t *par)
+{
+ Double_t lMass=0;
+ lMass = 1.5318; //Xi* mass
+
+
+ Double_t ldNdy = par[0]; // dN/dy
+ Double_t l2pi = 2*TMath::Pi(); // 2pi (cancels in return statement)
+ Double_t lTemp = par[1]; // Temperature
+ Double_t lPower = par[2]; // power=n
+
+ Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
+ Double_t lInPower = 1 + (TMath::Sqrt(x[0]*x[0]+lMass*lMass)-lMass) / (lPower*lTemp);
+
+ return l2pi * ldNdy * x[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
+}
+//________________________________________________________________________
+double myMtExpFunc(Double_t *x, Double_t *par)
+{
+ Double_t lMass=0;
+ lMass = 1.5318; //Xi* mass
+
+ return 2*TMath::Pi()*x[0]*par[0]*exp(-TMath::Sqrt(x[0]*x[0]+lMass*lMass)/par[1]);
+}
+//________________________________________________________________________
+double myPtExpFunc(Double_t *x, Double_t *par)
+{
+ return 2*TMath::Pi()*x[0]*par[0]*exp(-x[0]/par[1]);
+}
+//________________________________________________________________________
+double myBoltzFunc(Double_t *x, Double_t *par)
+{
+ Double_t lMass=0;
+ lMass = 1.5318; //Xi* mass
+ Double_t Mt = sqrt(pow(x[0],2)+pow(lMass,2));
+ return par[0]*x[0]*Mt*exp(-Mt/par[1]);
+
+}
+//________________________________________________________________________
+double myPowerLawPtFunc(Double_t *x, Double_t *par)
+{
+ return ( x[0]*par[0]*pow(1+x[0]/par[1], -par[2]) );
+}
+//________________________________________________________________________
+double myPowerLawFunc(Double_t *x, Double_t *par)
+{
+ return ( par[0]*pow(1+x[0]/par[1], -par[2]) );
+}
+//________________________________________________________________________
+double myBoltsBWFunc(Double_t *x, Double_t *par)
+{
+ double pT = x[0];
+ double mass = par[0];
+ double beta = par[1];
+ double temp = par[2];
+ double n = par[3];
+ double norm = par[4];
+
+ static TF1 * fIntBG = 0;
+ if(!fIntBG)
+ fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
+
+ fIntBG->SetParameters(mass, pT, beta, temp,n);
+ double result = fIntBG->Integral(0,1);
+ return result*norm;
+}
+//________________________________________________________________________
+double IntegrandBG(const double * x, const double* par){
+ double x0 = x[0];
+
+ double mass = par[0];
+ double pT = par[1];
+ double beta_max = par[2];
+ double temp = par[3];
+ Double_t n = par[4];
+
+ // Keep beta within reasonable limits
+ Double_t beta = beta_max * TMath::Power(x0, n);
+ if (beta > 0.9999999999999999) beta = 0.9999999999999999;
+
+ double mT = TMath::Sqrt(mass*mass+pT*pT);
+
+ double rho0 = TMath::ATanH(beta);
+ double arg00 = pT*TMath::SinH(rho0)/temp;
+ if (arg00 > 700.) arg00 = 700.; // avoid FPE
+ double arg01 = mT*TMath::CosH(rho0)/temp;
+ double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
+
+ // printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01);
+
+ return f0;
+}
--- /dev/null
+//gROOT->Reset();
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+#include "RooFit.h"
+#include "RooVoigtian.h"
+#include "RooAbsReal.h"
+#include "RooRealVar.h"
+#include "RooComplex.h"
+#include "RooMath.h"
+
+#include "TVector2.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TMath.h"
+#include "TText.h"
+#include "TRandom3.h"
+#include "TArray.h"
+#include "TLegend.h"
+#include "TStyle.h"
+#include "TMinuit.h"
+
+#include "TCanvas.h"
+#include "TLatex.h"
+#include "TImage.h"
+#include "TGaxis.h"
+#include "TGraphErrors.h"
+#include "TPaveStats.h"
+
+using namespace std;
+
+bool SaveToFile=kFALSE;
+int VariationType=0;// 0 = CutVariations, 1 = FitVariations
+//
+int CutVariation=0;// 0=Standard Cuts, 1-12 are Systematic deviations
+int FitVariation=1;// 1 for VoigtOnly with BkgSub, 2 for VoigtOnly without BkgSub, 3 for RightNorm, 4 Alternate BinCounting Range
+TString *outfilename = new TString();
+//
+// below are default (standard) settings
+int POLdegree=1;// 1 with BkgSub, 5 without
+bool XiStarCase=kTRUE;// Standard: kTRUE
+bool VoigtFitOnly=kFALSE;// Standard: kFALSE
+bool bkgSubtract=kTRUE;// Standard: kTRUE
+bool LeftNorm=kTRUE;// Standard: kTRUE
+bool MCAssocCase=kTRUE;// Standard: kTRUE
+short ParticleCase=0;// 0 for particle+anti-particle. 1 for particle. 2 for anti-particle
+
+
+const int max_ptbins = 8; // 18 for Xi, 8 for XiStar
+double raprange[2]={-0.5,0.5};// for both Xi and XiStar (bin size = 0.1)
+
+double XiyieldRange[2]={1.313, 1.332};
+double FitRangeXi[2]={1.27,1.36};
+
+double XiStaryieldRange[2]={1.51, 1.56};
+double FitRangeXiStar[2]={1.48,1.59};
+
+double XiCountBound[4]={1.304,1.313,1.331,1.340};// -9sigma, -4.5sigma, +4.5sigma, +9sigma
+double XiStarCountBound[6]={1.481, 1.490, 1.523, 1.541, 1.570, 1.579};// 1 gamma range (Standard)
+double XiStarCountBound_SysRange[6]={1.481, 1.490, 1.528, 1.537, 1.570, 1.579};// 1/2*gamma range (Systematic Variation)
+
+double yieldRange[2]={0};
+bool reject;
+
+
+double PolFunction(double *x, double *par);
+double PolFunctionSpecialXi(double *x, double *par);
+double BWFunction(double *x, double *par);
+double BWplusPol(double *x, double *par);
+double GausFunction(double *x, double *par);
+double GausplusPol(double *x, double *par);
+double GausplusPolSpecialXi(double *x, double *par);
+double myLevyPt(Double_t *x, Double_t *par);
+double myExpFit(Double_t *x, Double_t *par);
+double Voigt(Double_t *x, Double_t *par);
+double VoigtplusPol(double *x, double *par);
+
+
+
+void Plot_macroXiStar(){
+
+ if(VariationType==0) {
+ outfilename->Append("CutVariation_");
+ *outfilename += CutVariation;
+ }
+ if(VariationType==1) {
+ outfilename->Append("FitVariation_");
+ *outfilename += FitVariation;
+ // 1 for VoigtOnly with BkgSub, 2 for VoigtOnly without BkgSub, 3 for RightNorm, 4 Alternate BinCounting Range
+ if(FitVariation==1) {VoigtFitOnly=kTRUE;}
+ if(FitVariation==2) {VoigtFitOnly=kTRUE; bkgSubtract=kFALSE; POLdegree=5;}
+ if(FitVariation==3) {LeftNorm=kFALSE;}
+ if(FitVariation==4) {for(int ii=0; ii<6; ii++) XiStarCountBound[ii]=XiStarCountBound_SysRange[ii];}
+ }
+ outfilename->Append(".root");
+
+ double FitRange[2]={0};
+ if(XiStarCase) {
+ FitRange[0] = FitRangeXiStar[0]; FitRange[1] = FitRangeXiStar[1];
+ yieldRange[0] = XiStaryieldRange[0]; yieldRange[1] = XiStaryieldRange[1];
+ }
+ else {
+ FitRange[0] = FitRangeXi[0]; FitRange[1] = FitRangeXi[1];
+ yieldRange[0] = XiyieldRange[0]; yieldRange[1] = XiyieldRange[1];
+ }
+
+
+
+ TGaxis::SetMaxDigits(3);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(111);
+ //gStyle->SetOptFit(1111);
+
+
+ TCanvas *can = new TCanvas("can","can",13,34,700,500);
+ gStyle->SetOptFit(1111);
+ can->Range(-1.25,-0.2625,11.25,2.3625);
+ can->SetFillColor(10);
+ can->SetBorderMode(0);
+ can->SetBorderSize(2);
+ //can->SetGridx();
+ //can->SetGridy();
+ can->SetFrameFillColor(0);
+ can->SetFrameBorderMode(0);
+ can->SetFrameBorderMode(0);
+ can->Divide(4,2);
+
+ TLegend *legend = new TLegend(.1,.7,.35,.9,NULL,"brNDC");
+ legend->SetBorderSize(1);
+ legend->SetTextSize(.04);// small .03; large .036
+ //legend->SetLineColor(0);
+ legend->SetFillColor(0);
+ TLegend *legend2 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
+ legend2->SetBorderSize(1);
+ legend2->SetTextSize(.04);// small .03; large .036
+ //legend2->SetLineColor(0);
+ legend2->SetFillColor(0);
+
+
+
+ int PARBINS_temp;
+ PARBINS_temp = POLdegree+1+4;
+ const int PARBINS=PARBINS_temp;
+ int offset=0;
+ offset=4;
+
+ double bkg_params[POLdegree+1];
+
+ TString *polname=new TString();
+ if(POLdegree==1) polname->Append("pol1(0)");
+ if(POLdegree==2) polname->Append("pol2(0)");
+ if(POLdegree==3) polname->Append("pol3(0)");
+ if(POLdegree==4) polname->Append("pol4(0)");
+ if(POLdegree==5) polname->Append("pol5(0)");
+
+
+ TF1 *myPol[max_ptbins];
+ for(int i=0; i<max_ptbins; i++){
+ TString *name = new TString("myPol");
+ *name +=i+1;
+ myPol[i] = new TF1(name->Data(),polname->Data(),0,2);
+ myPol[i]->SetLineColor(3);
+ myPol[i]->SetRange(FitRange[0],FitRange[1]);
+ }
+
+
+ //
+
+
+ // use 10b + 10c + 10d for final result
+ //TFile *myfile = new TFile("Results/Real_LHC10b_lego.root","READ");
+ //TFile *myfile = new TFile("Results/Real_LHC10c_lego.root","READ");
+ //TFile *myfile = new TFile("Results/Real_LHC10d_lego.root","READ");
+ TFile *myfile = new TFile("Results/RealMerge_10b_10c_10d_lego.root","READ");// accepted version
+
+
+ // use 10d1 + 10d4 + 10f6a for final result
+ //TFile *myfile2 = new TFile("Results/MC_LHC10d1_lego.root","READ");
+ //TFile *myfile2 = new TFile("Results/MC_LHC10d4_lego.root","READ");
+ //TFile *myfile2 = new TFile("Results/MC_LHC10f6a_lego.root","READ");// Pythia
+ TFile *myfile2 = new TFile("Results/MCMerge_10d1_10d4_10f6a_lego.root","READ");// accepted version
+ //TFile *myfile2 = new TFile("Results/MC_LHC10f6_lego.root","READ");// Phojet
+
+ double David_XiMinus_eff[18]={0.012631, 0.0321763, 0.0473262, 0.059645, 0.0758325, 0.0926316, 0.109662, 0.120194, 0.145389, 0.17972, 0.218773, 0.248567, 0.287159, 0.298501, 0.287736, 0.278715, 0.241559, 0.168857};
+ double David_XiMinus_eff_e[18]={0.00038781, 0.000951899, 0.00123674, 0.00147704, 0.0017635, 0.00209177, 0.00243421, 0.00272509, 0.00236159, 0.0030364, 0.00329984, 0.00389785, 0.00508514, 0.00597671, 0.00825827, 0.0124455, 0.0167743, 0.0190201};
+ double David_XiPlus_eff[18]={0.010763, 0.0294655, 0.0429606, 0.0558855, 0.0722606, 0.0901825, 0.0982448, 0.114756, 0.140581, 0.172746, 0.208525, 0.247291, 0.28522, 0.29594, 0.28497, 0.276188, 0.249522, 0.210375};
+ double David_XiPlus_eff_e[18]={0.000382058, 0.000962961, 0.00122875, 0.00147933, 0.00178481, 0.00212508, 0.00236694, 0.00273767, 0.00238592, 0.00306165, 0.00329165, 0.0039834, 0.00516792, 0.00603728, 0.00836173, 0.0124634, 0.0168862, 0.0220455};
+
+ TDirectoryFile *tdir=(TDirectoryFile*)myfile->Get("PWGLF.outputXiStarAnalysis.root");
+ TList *List1=(TList*)tdir->Get("XiStarOutput");
+ //TList *List1=(TList*)myfile->Get("MyList");
+ myfile->Close();
+
+ TDirectoryFile *tdir2=(TDirectoryFile*)myfile2->Get("PWGLF.outputXiStarAnalysis.root");
+ TList *List2=(TList*)tdir2->Get("XiStarOutput");
+ //TList *List2=(TList*)myfile2->Get("MyList");
+ myfile2->Close();
+
+ TH1F *Events_PhysSel=(TH1F*)(List1->FindObject("fMultDist1"));
+ TH1F *Events_postPV=(TH1F*)(List1->FindObject("fMultDist3"));
+ TH1F *EventsMC_PhysSel=(TH1F*)(List2->FindObject("fMultDist1"));
+ TH1F *EventsMC_postPV=(TH1F*)(List2->FindObject("fMultDist3"));
+
+ cout<<"# Events passing PhysSelection Cuts = "<<Events_PhysSel->GetEntries()<<endl;
+ cout<<"# Events passing PhysSel + Vz Cut + PileUp Cut = "<<Events_postPV->GetEntries()<<endl;
+ cout<<"# MC Events passing Zvertex and PhysSelection Cuts = "<<EventsMC_postPV->GetEntries()<<endl;
+
+ TString *name = new TString("fXi_0");
+ TH3F *Xi3 = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXibar_0");
+ TH3F *Xi_bar3 = (TH3F*)(List1->FindObject(name->Data()));
+
+ name = new TString("fXiMinusPiPlus_"); *name += CutVariation;
+ TH3F *XiMinus_piPlus3 = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiMinusPiMinus_"); *name += CutVariation;
+ TH3F *XiMinus_piMinus3 = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiPlusPiPlus_"); *name += CutVariation;
+ TH3F *XiPlus_piPlus3 = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiPlusPiMinus_"); *name += CutVariation;
+ TH3F *XiPlus_piMinus3 = (TH3F*)(List1->FindObject(name->Data()));
+
+ name = new TString("fXiMinusPiPlusbkg_"); *name += CutVariation;
+ TH3F *XiMinus_piPlus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiMinusPiMinusbkg_"); *name += CutVariation;
+ TH3F *XiMinus_piMinus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiPlusPiPlusbkg_"); *name += CutVariation;
+ TH3F *XiPlus_piPlus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
+ name = new TString("fXiPlusPiMinusbkg_"); *name += CutVariation;
+ TH3F *XiPlus_piMinus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
+
+
+ TH3F *MCinput_Xi_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXi1"));
+ TH3F *MCinput_Xi_bar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXibar1"));
+ TH3F *MCinput_XiStar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXiStar1"));
+ TH3F *MCinput_XiStar_bar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXiStarbar1"));
+ TH3F *MCinput_XiStar = (TH3F*)(List2->FindObject("fMCinputTotalXiStar3"));
+ TH3F *MCinput_XiStar_bar = (TH3F*)(List2->FindObject("fMCinputTotalXiStarbar3"));
+ TH3F *MCinput_Xi = (TH3F*)(List2->FindObject("fMCinputTotalXi3"));
+ TH3F *MCinput_Xi_bar = (TH3F*)(List2->FindObject("fMCinputTotalXibar3"));
+
+ int binYL = MCinput_XiStar->GetYaxis()->FindBin(raprange[0]+0.01);
+ int binYH = MCinput_XiStar->GetYaxis()->FindBin(raprange[1]-0.01);
+ int binML = MCinput_XiStar->GetZaxis()->FindBin(1.51);
+ int binMH = MCinput_XiStar->GetZaxis()->FindBin(1.55);
+
+ TH1D *MCinput_Spectrum;
+ TH1D *MCinput_Spectrum_bar;
+ if(XiStarCase){
+ MCinput_Spectrum= (TH1D*)MCinput_XiStar->ProjectionX("MCinput_Spectrum",binYL,binYH, binML,binMH);
+ MCinput_Spectrum_bar = (TH1D*)MCinput_XiStar_bar->ProjectionX("MCinput_Spectrum_bar",binYL,binYH, binML,binMH);
+ }else{
+ MCinput_Spectrum= (TH1D*)MCinput_Xi->ProjectionX("MCinput_Spectrum",binYL,binYH, binML,binMH);
+ MCinput_Spectrum_bar = (TH1D*)MCinput_Xi_bar->ProjectionX("MCinput_Spectrum_bar",binYL,binYH, binML,binMH);
+ }
+ MCinput_Spectrum->Add(MCinput_Spectrum_bar);
+
+ name = new TString("fMCrecXiMinusPiPlus_"); *name += CutVariation;
+ TH3F *MCrec_XiStar = (TH3F*)(List2->FindObject(name->Data()));
+ name = new TString("fMCrecXiPlusPiMinus_"); *name += CutVariation;
+ TH3F *MCrec_XiStar_bar = (TH3F*)(List2->FindObject(name->Data()));
+
+
+ TH3F *MCrec_Xi;
+ TH3F *MCrec_Xi_bar;
+ if(MCAssocCase){
+ MCrec_Xi = (TH3F*)(List2->FindObject("fMCrecXi"));
+ MCrec_Xi_bar = (TH3F*)(List2->FindObject("fMCrecXibar"));
+ }else{
+ name = new TString("fXi_0");
+ MCrec_Xi = (TH3F*)(List2->FindObject(name->Data()));
+ name = new TString("fXibar_0");
+ MCrec_Xi_bar = (TH3F*)(List2->FindObject(name->Data()));
+ }
+
+
+
+
+ double ptedges_Dhevan[9]={0.8, 1.2, 1.6, 2.0, 2.4, 3.2, 4.0, 4.8, 5.6};// standard
+ double pt_points_Dhevan[8]={1.0, 1.4, 1.8, 2.2, 2.8, 3.6, 4.4, 5.2};// standard
+ double pt_points_e_Dhevan[8]={.2, .2, .2, .2, .4, .4, .4, .4};// standard
+ //double ptedges_Dhevan[17]={0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 4.8, 5.2, 5.6};
+ //double pt_points_Dhevan[16]={0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.6, 3.0, 3.4, 3.8, 4.2, 4.6, 5.0, 5.4};
+ //double pt_points_e_Dhevan[16]={.1, .1, .1, .1, .1, .1, .1, .1, .2, .2, .2, .2, .2, .2, .2, .2};
+
+ double ptedges_David[19]={0.6, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.7, 1.9, 2.2, 2.6, 3.1, 3.9, 4.9, 6, 7.2, 8.5};
+ double pt_points_David[18]={0.7, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.6, 1.8, 2.05, 2.4, 2.85, 3.5, 4.4, 5.45, 6.6, 7.85};
+ double pt_points_e_David[18]={0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.3, 0.4, 0.5, 0.8, 1, 1.1, 1.2, 1.3};
+
+
+ double ptedges[max_ptbins+1+1]={0};
+ double pt_points[max_ptbins]={0};
+ double pt_points_e[max_ptbins]={0};
+ TString *names[max_ptbins];
+
+ for(int i=0; i<max_ptbins; i++){
+ if(XiStarCase){
+ //ptedges[i] = ptedges_Dhevan[i];
+ //ptedges[max_ptbins] = ptedges_Dhevan[max_ptbins];
+ pt_points[i] = pt_points_Dhevan[i];
+ pt_points_e[i] = pt_points_e_Dhevan[i];
+ names[i] = new TString("pt bin ");
+ *names[i] += i+1;
+ }
+ else{
+ //ptedges[i] = ptedges_David[i];
+ //ptedges[max_ptbins] = ptedges_David[max_ptbins];
+ pt_points[i] = pt_points_David[i];
+ pt_points_e[i] = pt_points_e_David[i]/2.;
+ names[i] = new TString();
+ *names[i] += float(ptedges_David[i]);
+ names[i]->Append(" < pt < ");
+ *names[i] += float(ptedges_David[i+1]);
+ }
+ }
+ for(int i=0; i<max_ptbins+1; i++){
+ if(XiStarCase){
+ if(i==0) ptedges[i]=0;
+ else {
+ ptedges[i] = ptedges_Dhevan[i-1];
+ ptedges[max_ptbins+1] = ptedges_Dhevan[max_ptbins];
+ }
+ }else {
+ if(i==0) ptedges[i]=0;
+ else {
+ ptedges[i] = ptedges_David[i-1];
+ ptedges[max_ptbins+1] = ptedges_David[max_ptbins];
+ }
+ }
+ }
+
+ double Eff_corr_Xiplus[max_ptbins]={0};
+ double Eff_corr_Ximinus[max_ptbins]={0};
+
+ if(!XiStarCase){
+ // David binning
+ double geant3fluka_corr_Xiplus[18]={0.8646, 0.8958, 0.9063, 0.9135, 0.9190, 0.9241, 0.9285, 0.9321, 0.9370, 0.9428, 0.9485, 0.9548, 0.9612, 0.9677, 0.9745, 0.9802, 0.9847, 0.9884};
+ double geant3fluka_corr_Ximinus[18]={0.9843, 0.9903, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906};
+ for(int i=0; i<max_ptbins; i++) {
+ Eff_corr_Xiplus[i] = geant3fluka_corr_Xiplus[i];
+ Eff_corr_Ximinus[i] = geant3fluka_corr_Ximinus[i];
+ }
+ }
+ if(XiStarCase){
+ double geant3fluka_corr_Xiplus[8]={0.9118, 0.9305, 0.9427, 0.9514, 0.9601, 0.9686, 0.9746, 0.9792};// standard
+ double geant3fluka_corr_Ximinus[8]={0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906};// standard
+ //double geant3fluka_corr_Xiplus[16]={0.9025, 0.9166, 0.9264, 0.9339, 0.9400, 0.9452, 0.9496, 0.9534, 0.9579, 0.9631, 0.9672, 0.9706, 0.9737, 0.9760, 0.9784, 0.9803};
+ //double geant3fluka_corr_Ximinus[16]={0.9905, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906};
+ for(int i=0; i<max_ptbins; i++) {
+ Eff_corr_Xiplus[i] = geant3fluka_corr_Xiplus[i];
+ Eff_corr_Ximinus[i] = geant3fluka_corr_Ximinus[i];
+ }
+ }
+
+ double Eff[max_ptbins]={0};
+ double Eff_e[max_ptbins]={0};
+
+
+ double spectrum[max_ptbins]={0};
+ double spectrum_e[max_ptbins]={0};
+ double chi2perNDF[max_ptbins]={0};
+ double yield[max_ptbins]={0};
+ double yield_e[max_ptbins]={0};
+ double width[max_ptbins]={0};
+ double width_e[max_ptbins]={0};
+ double mass[max_ptbins]={0};
+ double mass_e[max_ptbins]={0};
+ double MCwidth[max_ptbins]={0};
+ double MCwidth_e[max_ptbins]={0};
+ double MCmass[max_ptbins]={0};
+ double MCmass_e[max_ptbins]={0};
+ double MissedYield[max_ptbins]={0};
+ double MissedYield_e[max_ptbins]={0};
+
+ TH1F *MCrec_1d[max_ptbins];
+ TH1F *MCrec_1d_bar[max_ptbins];
+ TF1 *MCfit = new TF1("MCfit",Voigt,FitRange[0],FitRange[1],4);
+ double LostRatio=0, LostRatio_e=0;
+
+ for(int ptbin=0; ptbin<max_ptbins; ptbin++){
+ can->cd(ptbin+1);
+ cout<<"+++++++++++++++ ptbin = "<<ptbin<<" +++++++++++++++++"<<endl;
+ int rapbinFirst = Xi3->GetYaxis()->FindBin(raprange[0]+0.01);// bin 16
+ int rapbinLast = Xi3->GetYaxis()->FindBin(raprange[1]-0.01);// bin 25
+ int ptbinFirst = Xi3->GetXaxis()->FindBin(ptedges[ptbin+1]+0.01);
+ int ptbinLast = Xi3->GetXaxis()->FindBin(ptedges[ptbin+1+1]-0.01);
+
+ ////////////////////////////
+ // Efficiency calcualtion
+ double rec_num=0, input_num=0;
+ double rec_num1=0, input_num1=0;
+ double rec_num2=0, input_num2=0;
+ if(XiStarCase){
+ int massbinFirst = 1;
+ int massbinLast = MCrec_XiStar->GetNbinsZ();
+ rec_num1 = MCrec_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ input_num1 = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ rec_num2 = MCrec_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ input_num2 = MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+
+ // set ranges for projections later on
+ MCrec_XiStar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ MCrec_XiStar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ MCrec_XiStar->GetZaxis()->SetRange(massbinFirst,massbinLast);
+ MCrec_XiStar_bar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ MCrec_XiStar_bar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ MCrec_XiStar_bar->GetZaxis()->SetRange(massbinFirst,massbinLast);
+ // Estimate Lost number of Xi* from PV cut
+ double num=0, den=0;
+ if(ParticleCase==0){
+ num = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast) + MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ den = MCinput_XiStar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast) + MCinput_XiStar_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ }else if(ParticleCase==1){
+ num = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ den = MCinput_XiStar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ }else {
+ num = MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ den = MCinput_XiStar_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
+ }
+ LostRatio = num/den;
+ LostRatio_e = sqrt(pow(sqrt(num)/den,2)+pow(sqrt(den)*num/pow(den,2),2));
+
+ ///////////////////////////////////////
+ // MC fitting for mass and sigma values
+ TString *nameMCrec=new TString("MCrec_1d_");
+ TString *nameMCrecbar=new TString("MCrecbar_1d_");
+ *nameMCrec += ptbin;
+ *nameMCrecbar += ptbin;
+ MCrec_1d[ptbin] = (TH1F*)MCrec_XiStar->ProjectionZ(nameMCrec->Data(),ptbinFirst,ptbinLast, rapbinFirst,rapbinLast);
+ MCrec_1d_bar[ptbin] = (TH1F*)MCrec_XiStar->ProjectionZ(nameMCrecbar->Data(),ptbinFirst,ptbinLast, rapbinFirst,rapbinLast);
+ MCrec_1d[ptbin]->Add(MCrec_1d_bar[ptbin]);
+ MCrec_1d[ptbin]->GetXaxis()->SetRangeUser(1.46, 1.62);
+ MCrec_1d[ptbin]->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
+ MCrec_1d[ptbin]->SetTitle("p_{T} bin 8");
+ MCfit->SetParameter(0,.5);// mean
+ MCfit->SetParameter(1,1.532);// mean
+ MCfit->SetParameter(2,.002);// Gaussian sigma
+ MCfit->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
+ MCfit->SetParLimits(1,1.52,1.54);
+ MCfit->SetParLimits(2,0.00001,.01);
+ MCfit->FixParameter(3,.0091);// BW width
+ MCrec_1d[ptbin]->Fit(MCfit,"IMQ","",FitRange[0],FitRange[1]);
+ MCwidth[ptbin] = MCfit->GetParameter(2);
+ MCwidth_e[ptbin] = MCfit->GetParError(2);
+ MCmass[ptbin] = MCfit->GetParameter(1);
+ MCmass_e[ptbin] = MCfit->GetParError(1);
+ //continue;// use this only to plot MC yields one by one
+ ///////////////////////////////////////////
+
+ }else {// Xi
+ rec_num1 = MCrec_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ input_num1 = MCinput_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCinput_Xi->GetZaxis()->GetNbins());
+ rec_num2 = MCrec_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi_bar->GetZaxis()->GetNbins());
+ input_num2 = MCinput_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCinput_Xi_bar->GetZaxis()->GetNbins());
+ //
+ MCrec_Xi->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ MCrec_Xi->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ MCrec_Xi_bar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ MCrec_Xi_bar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ MCrec_1d[ptbin]=(TH1F*)MCrec_Xi->Project3D("z");
+ MCrec_1d_bar[ptbin]=(TH1F*)MCrec_Xi_bar->Project3D("z");
+
+ rec_num1 = MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[1]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[2]-.0001));
+ rec_num1 -= MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[0]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[1]-.0001));
+ rec_num1 -= MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[2]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[3]-.0001));
+
+ rec_num2 = MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[1]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[2]-.0001));
+ rec_num2 -= MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[0]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[1]-.0001));
+ rec_num2 -= MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[2]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[3]-.0001));
+
+ // Estimate Lost number of Xi from PV cut
+ if(ParticleCase==0){
+ LostRatio = MCinput_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins()) + MCinput_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ LostRatio /= MCinput_Xi_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins()) + MCinput_Xi_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ }else if(ParticleCase==1){
+ LostRatio = MCinput_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ LostRatio /= MCinput_Xi_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ }else {
+ LostRatio = MCinput_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ LostRatio /= MCinput_Xi_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
+ }
+
+
+ }// else Xi case
+ cout<<"LostRatio = "<<LostRatio<<" +- "<<LostRatio_e<<endl;
+
+
+ // Apply Geant3/fluka correction
+ rec_num1 /= Eff_corr_Ximinus[ptbin];
+ rec_num2 /= Eff_corr_Xiplus[ptbin];
+
+
+ if(ParticleCase == 0){
+ rec_num = rec_num1 + rec_num2;
+ input_num = input_num1 + input_num2;
+ }
+ if(ParticleCase == 1){
+ rec_num = rec_num1;
+ input_num = input_num1;
+ }
+ if(ParticleCase == 2){
+ rec_num = rec_num2;
+ input_num = input_num2;
+ }
+
+
+ Eff[ptbin] = rec_num/input_num;
+ Eff_e[ptbin] = sqrt( pow(sqrt(rec_num)/input_num,2) + pow(sqrt(input_num)*rec_num/pow(input_num,2),2));
+
+ cout<<"Efficiency = "<<Eff[ptbin]<<" +- "<<Eff_e[ptbin]<<" MC_rec_count = "<<rec_num<<endl;
+ //cout<<"Reassigning Eff in last pt bin to second to last!!!!!!!!!!!!!"<<endl; Eff[7] = Eff[6];
+
+
+
+ //cout<<"Re-assigning to David's Eff"<<endl; Eff[ptbin] = David_XiMinus_eff[ptbin];
+
+
+ reject=kTRUE;
+ TF1 *myBkg=new TF1("myBkg",PolFunction,FitRange[0],FitRange[1],POLdegree+1);
+ //TF1 *myBkg=new TF1("myBkg",PolFunctionSpecialXi,1.27,1.36,POLdegree+1);// special bkg for Xi
+ myBkg->SetParameters(6,0);
+ myBkg->SetLineColor(3);
+
+
+ Xi3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);// pt bins
+ Xi_bar3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ Xi3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ Xi_bar3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+
+ XiMinus_piPlus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiPlus_piMinus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiMinus_piMinus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiPlus_piPlus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+
+ XiMinus_piPlus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiPlus_piMinus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiMinus_piMinus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiPlus_piPlus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+
+
+ XiMinus_piPlus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiPlus_piMinus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiMinus_piMinus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+ XiPlus_piPlus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
+
+ XiMinus_piPlus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiPlus_piMinus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiMinus_piMinus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+ XiPlus_piPlus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
+
+
+ TH1F *Xi=(TH1F*)Xi3->Project3D("z");
+ TH1F *Xibar=(TH1F*)Xi_bar3->Project3D("z");
+ if(ParticleCase==0) Xi->Add(Xibar);// to add particle and anti-particle
+
+ TH1F *XiMinus_piPlus = (TH1F*)XiMinus_piPlus3->Project3D("z");
+ TH1F *XiPlus_piMinus = (TH1F*)XiPlus_piMinus3->Project3D("z");
+ TH1F *XiMinus_piMinus = (TH1F*)XiMinus_piMinus3->Project3D("z");
+ TH1F *XiPlus_piPlus = (TH1F*)XiPlus_piPlus3->Project3D("z");
+ if(ParticleCase==0) {
+ XiMinus_piPlus->Add(XiPlus_piMinus);
+ XiMinus_piMinus->Add(XiPlus_piPlus);
+ }
+
+
+ TH1F *XiMinus_piPlus_bkg = (TH1F*)XiMinus_piPlus3_bkg->Project3D("z");
+ TH1F *XiPlus_piMinus_bkg = (TH1F*)XiPlus_piMinus3_bkg->Project3D("z");
+ TH1F *XiMinus_piMinus_bkg = (TH1F*)XiMinus_piMinus3_bkg->Project3D("z");
+ TH1F *XiPlus_piPlus_bkg = (TH1F*)XiPlus_piPlus3_bkg->Project3D("z");
+ if(ParticleCase==0) {
+ XiMinus_piPlus_bkg->Add(XiPlus_piMinus_bkg);
+ XiMinus_piMinus_bkg->Add(XiPlus_piPlus_bkg);
+ }
+
+ TH1F *os=0x0;
+ TH1F *os_bkg=0x0;
+ if(XiStarCase) {
+ if(ParticleCase==0 || ParticleCase==1) {
+ os=(TH1F*)XiMinus_piPlus->Clone();
+ os_bkg=(TH1F*)XiMinus_piPlus_bkg->Clone();
+ }
+ if(ParticleCase==2) os=(TH1F*)XiPlus_piMinus->Clone();
+ }else {
+ if(ParticleCase==0 || ParticleCase==1) os=(TH1F*)Xi->Clone();
+ if(ParticleCase==2) os=(TH1F*)Xibar->Clone();
+ }
+
+
+ os->Sumw2();
+ TH1F *ss=(TH1F*)XiMinus_piMinus->Clone();
+ TH1F *ss_bkg=(TH1F*)XiMinus_piMinus_bkg->Clone();
+ ss_bkg->SetLineColor(2);
+
+ ss->SetLineColor(4);
+ ss->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
+ ss->GetYaxis()->SetTitle("Counts");
+ //ss->SetTitle(name[ptbin]);
+ ss->SetTitle(names[ptbin]->Data());
+
+ os->SetLineColor(4);
+ os->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
+ os->GetYaxis()->SetTitle("Counts");
+ //os->SetTitle(name[ptbin]);
+ os->SetTitle(names[ptbin]->Data());
+
+
+ double scalefactor=1;
+ if(XiStarCase) {
+ os_bkg->SetLineColor(2);
+ os_bkg->Sumw2();
+ double leftS=0, rightS=0;
+ if(LeftNorm) {leftS=1.49; rightS=1.51;}
+ else {leftS=1.56; rightS=1.58;}
+
+ scalefactor = os->Integral(os->GetXaxis()->FindBin(leftS),os->GetXaxis()->FindBin(rightS));
+ scalefactor /= os_bkg->Integral(os_bkg->GetXaxis()->FindBin(leftS),os_bkg->GetXaxis()->FindBin(rightS));
+ os_bkg->Scale(scalefactor);
+ if(bkgSubtract) os->Add(os_bkg,-1);
+ }
+
+
+ /////////////////////////////////////
+ // 1st Iteration
+ reject=kTRUE;
+ if(bkgSubtract) {
+ for(int polbin=0; polbin<POLdegree+1-2; polbin++) myBkg->FixParameter(polbin+2,0);
+ myBkg->SetParameter(0,0);
+ myBkg->SetParameter(1,0);
+ }
+ /////////////
+ os->Fit(myBkg,"NQ","",FitRange[0],FitRange[1]);// 1st fit
+ /////////////
+
+
+ //ss->Draw();
+ //ss_bkg->Scale(ss->Integral()/ss_bkg->Integral());
+ //ss_bkg->Draw("same");
+
+ TF1 *FullFit_t2;
+ TF1 *FullFit_t6;
+
+
+
+ // Voigt
+ FullFit_t2=new TF1("FullFit_t2",VoigtplusPol,FitRange[0],FitRange[1],PARBINS+1);
+ FullFit_t2->SetParameter(0,1.0);// mean
+ FullFit_t2->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
+ FullFit_t2->SetParameter(1,1.532);// mean
+ FullFit_t2->SetParLimits(1,1.5,1.56);
+ //FullFit_t2->FixParameter(2,0);// Gaussian sigma
+ FullFit_t2->SetParameter(2,.003);// Gaussian sigma
+ FullFit_t2->SetParLimits(2,0.0001,.015);
+ FullFit_t2->FixParameter(3,.0091);// BW width
+ //FullFit_t2->SetParameter(3,.0091);// BW width
+ //FullFit_t2->SetParLimits(3,.001,.015);
+ FullFit_t2->SetParName(0,"Norm");
+ FullFit_t2->SetParName(1,"Mass");
+ FullFit_t2->SetParName(2,"#sigma");
+ FullFit_t2->SetParName(3,"#Gamma");
+
+ FullFit_t2->SetParameter(4,0);
+ FullFit_t2->SetParameter(5,0);
+ FullFit_t2->SetParameter(6,0);
+ FullFit_t2->SetParameter(7,0);
+ FullFit_t2->SetParameter(8,0);
+ FullFit_t2->SetParameter(9,0);
+
+
+
+
+ if(XiStarCase) os->GetXaxis()->SetRange(1000*(FitRangeXiStar[0]-1.4),1000*(FitRangeXiStar[1]-1.4));
+ else os->GetXaxis()->SetRange(1000*(FitRangeXi[0]-1.2),1000*(FitRangeXi[1]+-1.2));
+
+ /////////////////////////////////////////////
+ // 2nd interation
+ reject=kFALSE;
+ double pars[PARBINS];
+ double pars_e[PARBINS];
+
+ for(int polbin=0; polbin<PARBINS-offset; polbin++) bkg_params[polbin] = myBkg->GetParameter(polbin);
+
+
+ for(int polbin=0; polbin<PARBINS-offset; polbin++) FullFit_t2->FixParameter(offset+polbin, myBkg->GetParameter(polbin));
+ //////////////////
+ os->Fit(FullFit_t2,"IMQ","",FitRange[0],FitRange[1]);// 2nd fit
+ //////////////////
+ for(int polbin=0; polbin<PARBINS; polbin++) pars[polbin] = FullFit_t2->GetParameter(polbin);
+ ////////////////////////////////////////
+
+
+
+ ///////////////////////////////////////
+ // 3rd iteration
+ FullFit_t6=new TF1("FullFit_t6",VoigtplusPol,FitRange[0],FitRange[1],PARBINS);
+ FullFit_t6->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
+ FullFit_t6->SetParLimits(1,1.5,1.56);
+ FullFit_t6->SetParLimits(2,0.0001,.015);
+ FullFit_t6->FixParameter(3,.0091);// BW width
+ FullFit_t6->SetParName(0,"Norm");
+ FullFit_t6->SetParName(1,"Mass");
+ FullFit_t6->SetParName(2,"#sigma");
+ FullFit_t6->SetParName(3,"#Gamma");
+
+
+ for(int polbin=0; polbin<PARBINS; polbin++) FullFit_t6->SetParameter(polbin, pars[polbin]);
+ if(bkgSubtract) {
+ for(int polbin=0; polbin<PARBINS-offset-2; polbin++) FullFit_t6->FixParameter(polbin+offset+2,0);
+ FullFit_t6->SetParameter(offset,0);
+ FullFit_t6->SetParameter(offset+1,0);
+ }
+ //////////////////
+ os->Fit(FullFit_t6,"IMEQ","",FitRange[0],FitRange[1]);// 3rd fit
+ //////////////////
+ for(int polbin=0; polbin<PARBINS-offset; polbin++) bkg_params[polbin] = FullFit_t6->GetParameter(polbin+offset);
+ for(int polbin=0; polbin<PARBINS; polbin++) pars[polbin] = FullFit_t6->GetParameter(polbin);
+ for(int polbin=0; polbin<PARBINS; polbin++) pars_e[polbin] = FullFit_t6->GetParError(polbin);
+
+
+
+ /////////////////////////////////////////
+ if(!bkgSubtract) os_bkg->Draw("same");// XiStar
+ for(int polbin=0; polbin<POLdegree+1; polbin++) myBkg->SetParameter(polbin, bkg_params[polbin]);
+
+ for(int polbin=0; polbin<POLdegree+1; polbin++) myPol[ptbin]->FixParameter(polbin, bkg_params[polbin]);
+ myPol[ptbin]->DrawCopy("same");
+ if((ptbin+1)==max_ptbins){
+ legend->AddEntry(os,"opp-charge","l");
+ legend->AddEntry(os_bkg,"opp-charge event mixing","l");
+ //legend->Draw("same");
+ }
+
+
+ double temp_yield=0, temp_yield_e=0;
+ MissedYield[ptbin] = 1000.*(FullFit_t6->Integral(0.,XiStarCountBound[2]) + FullFit_t6->Integral(XiStarCountBound[3],100.) - myBkg->Integral(0.,XiStarCountBound[2]) - myBkg->Integral(XiStarCountBound[3],100.));
+ MissedYield_e[ptbin] = MissedYield[ptbin]*FullFit_t6->GetParError(0)/FullFit_t6->GetParameter(0);
+ //cout<<"Signal/Bkg = "<<temp_yield/(1000.*myBkg->Integral(XiStarYieldCountRange[0],XiStarYieldCountRange[1]))<<endl;
+ //cout<<"Significance = "<<temp_yield/sqrt(temp_yield)<<endl;
+
+ if(VoigtFitOnly) {
+ if(!bkgSubtract) temp_yield = 1000.*(FullFit_t6->Integral(1.48,1.59) - myBkg->Integral(1.48,1.59));// must have fit range limits for high degree polynomial bkg
+ else temp_yield = 1000.*(FullFit_t6->Integral(0.,100.) - myBkg->Integral(0.,100.));// yield from Voigt fit (1.48-1.59 is function fit range)
+ temp_yield_e = temp_yield*FullFit_t6->GetParError(0)/FullFit_t6->GetParameter(0);
+ }else {// yield from bin counting
+ if(XiStarCase){
+ temp_yield = os->Integral(os->GetXaxis()->FindBin(XiStarCountBound[2]),os->GetXaxis()->FindBin(XiStarCountBound[3]-.0001));
+ temp_yield -= myBkg->Integral(XiStarCountBound[2],XiStarCountBound[3]);
+ cout<<"Missed yield fraction = "<<MissedYield[ptbin]/(temp_yield+MissedYield[ptbin])<<endl;
+ temp_yield += MissedYield[ptbin];
+ for(int massBin=os->GetXaxis()->FindBin(XiStarCountBound[2]); massBin<=os->GetXaxis()->FindBin(XiStarCountBound[3]); massBin++){
+ temp_yield_e += pow(os->GetBinError(massBin),2);
+ }
+ temp_yield_e += pow(MissedYield_e[ptbin],2);
+ temp_yield_e = sqrt(temp_yield_e);
+ cout<<"Included/Total = "<<(temp_yield-MissedYield[ptbin])/temp_yield<<endl;
+ }else {// Xi case
+ temp_yield = os->Integral(os->GetXaxis()->FindBin(XiCountBound[1]),os->GetXaxis()->FindBin(XiCountBound[2]-.0001));
+ temp_yield -= os->Integral(os->GetXaxis()->FindBin(XiCountBound[0]),os->GetXaxis()->FindBin(XiCountBound[1]-.0001));
+ temp_yield -= os->Integral(os->GetXaxis()->FindBin(XiCountBound[2]),os->GetXaxis()->FindBin(XiCountBound[3]-.0001));
+ }
+ }
+ cout<<"yield = "<<temp_yield<<" +- "<<temp_yield_e<<endl;
+
+
+
+ double base = Eff[ptbin]*LostRatio*(raprange[1]-raprange[0])*(2*pt_points_e[ptbin])*Events_PhysSel->GetEntries()/0.852;
+
+ spectrum[ptbin] = temp_yield;
+ spectrum[ptbin] /= base;
+ //
+ spectrum_e[ptbin] = pow(temp_yield_e/base,2);
+ spectrum_e[ptbin] += pow(temp_yield/base*Eff_e[ptbin]/Eff[ptbin],2);
+ spectrum_e[ptbin] += pow(temp_yield/base*LostRatio_e/LostRatio,2);
+ spectrum_e[ptbin] += pow(temp_yield/base*sqrt(Events_PhysSel->GetEntries())/Events_PhysSel->GetEntries(),2);
+ spectrum_e[ptbin] = sqrt(spectrum_e[ptbin]);
+ //spectrum_e[ptbin] = temp_yield/base*Eff_e[ptbin]/Eff[ptbin];// for eff errors only
+ if(ParticleCase==0) spectrum[ptbin] /= 2.;// (particle + cc)/2
+ if(ParticleCase==0) spectrum_e[ptbin] /= 2.;// (particle + cc)/2
+ //
+
+
+ cout<<"spectrum point = "<<spectrum[ptbin]<<" +- "<<spectrum_e[ptbin]<<endl;
+
+ yield[ptbin] = temp_yield;
+ yield_e[ptbin] = temp_yield_e;
+ mass[ptbin] = pars[1];
+ mass_e[ptbin] = pars_e[1];
+ width[ptbin] = pars[2];
+ width_e[ptbin] = pars_e[2];
+
+ chi2perNDF[ptbin] = FullFit_t6->GetChisquare()/double(FullFit_t6->GetNDF());
+
+ }
+
+
+
+ TH1D *RawYields = new TH1D("RawYields","",8,0.5,8.5);
+ for(int ii=0; ii<8; ii++){ RawYields->SetBinContent(ii+1, yield[ii]); RawYields->SetBinError(ii+1, yield_e[ii]);}
+ TH1D *Efficiency = new TH1D("Efficiency","",8,0.5,8.5);
+ for(int ii=0; ii<8; ii++){ Efficiency->SetBinContent(ii+1, Eff[ii]); Efficiency->SetBinError(ii+1, Eff_e[ii]);}
+
+ /*
+ TLatex *tex2 = new TLatex(1.6,1.5,"#Xi(1530)^{0} #rightarrow #Xi^{+-} + #pi^{-+}");
+ tex2->SetTextSize(.06);
+ tex2->Draw();
+
+ TLatex *tex3 = new TLatex(1.6,1.5,"ALICE Performance");
+ tex3->SetTextSize(.06);
+ tex3->SetTextColor(2);
+ tex3->Draw();
+
+ TLatex *tex4 = new TLatex(1.6,1.5,"12/09/2011");
+ tex4->SetTextSize(.04);
+ tex4->Draw();
+
+ TLatex *tex5 = new TLatex(1.6,1.5,"pp, #sqrt{s} = 7 TeV");
+ tex5->SetTextSize(.06);
+ tex5->Draw();
+
+ TPad *c1_1 = new TPad("c1_1", "c1_1",0.25,0.3,0.4,0.5);
+ c1_1->SetBorderMode(0);
+ c1_1->SetLineColor(0);
+ c1_1->SetFillColor(0);
+ c1_1->Draw();
+
+ c1_1->cd();
+ TImage *alicelogo = TImage::Open("../../Pictures/alice_logo_transparent2.png");
+ //alicelogo->SetConstRatio(kTRUE);
+ //alicelogo->SetImageQuality(TAttImage::kImgBest);
+ alicelogo->Draw("same");
+ */
+
+
+ TCanvas *can2 = new TCanvas("can2","can2",13,34,700,500);
+ gStyle->SetOptFit(111);
+ can2->Range(-1.25,-0.2625,11.25,2.3625);
+ can2->SetFillColor(10);
+ can2->SetBorderMode(0);
+ can2->SetBorderSize(2);
+ can2->SetGridx();
+ can2->SetGridy();
+ //can2->SetLogy();
+ can2->SetFrameFillColor(0);
+ can2->SetFrameBorderMode(0);
+ can2->SetFrameBorderMode(0);
+
+
+
+ TH1D *h_Xispectrum = new TH1D("h_Xispectrum","Xi spectrum",max_ptbins+1,ptedges);
+ h_Xispectrum->SetMarkerStyle(20);
+ h_Xispectrum->SetMarkerSize(1.0);
+ if(XiStarCase) h_Xispectrum->SetTitle("#Xi(1530) spectrum");
+ else h_Xispectrum->SetTitle("#Xi spectrum");
+ h_Xispectrum->SetTitle("pp #sqrt{s}= 7 TeV");
+ h_Xispectrum->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ h_Xispectrum->GetYaxis()->SetTitle("1/N_{E}d^{2}N/dydp_{T} (GeV/c)^{-1}");
+ //h_Xispectrum->GetYaxis()->SetTitle("(Real N_Xi-)/(MC N_Xi-)");
+ h_Xispectrum->SetMinimum(3.0e-6);
+ h_Xispectrum->SetMaximum(0.02);
+ h_Xispectrum->GetXaxis()->SetLimits(0,5.6);
+
+ h_Xispectrum->SetBinContent(1,+100);// first bin is just there to visualize pt=0
+ for(int i=0; i<max_ptbins; i++) {
+ h_Xispectrum->SetBinContent(i+1+1, spectrum[i]);
+ h_Xispectrum->SetBinError(i+1+1, spectrum_e[i]);
+ }
+
+ //h_Xispectrum->Draw();
+
+
+ double minFitpoint=0, maxFitpoint=0;
+ if(XiStarCase) {minFitpoint = 0.8; maxFitpoint = 5.6;}// 0.8 to 5.6
+ else {minFitpoint = 0.6; maxFitpoint = 8.5;}
+
+
+ TF1 *myLevy=new TF1("myLevy",myLevyPt,0,8.5,3);
+ myLevy->SetParName(0,"dN/dy");
+ myLevy->SetParName(1,"C");
+ myLevy->SetParName(2,"n");
+ myLevy->SetParameter(0,.008);
+ myLevy->SetParameter(1,.3);
+ myLevy->SetParameter(2,15);
+ myLevy->SetParLimits(0,.001,.2);
+ myLevy->SetParLimits(1,.1,1);
+ myLevy->SetParLimits(2,1,500);
+
+
+ h_Xispectrum->Draw();
+ h_Xispectrum->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
+ //h_Xispectrum->Fit(myExp,"IME","",minFitpoint,maxFitpoint);
+ //myLevy->Draw("same");
+ //TH1D *temp = (TH1D*)myLevy->GetHistogram();
+
+ //legend2->AddEntry(h_Xispectrum,"#Xi^{+}, This analysis","p");
+ //legend2->AddEntry(h_Xispectrum,"10b+10c+10d","p");
+ legend2->AddEntry(h_Xispectrum,"(#Xi(1530) + cc)/2","p");
+
+ cout<<"dN/dy = "<<myLevy->Integral(0,10)<<" , Covered range = "<<myLevy->Integral(0.8,5.6)<<endl;
+
+
+ double Data_Levy_Ratio[8]={0};
+ double Data_Levy_Ratio_e[8]={0};
+ for(int i=0; i<8; i++){
+ Data_Levy_Ratio[i] = (ptedges[i+2]-ptedges[i+1])*spectrum[i]/myLevy->Integral(ptedges[i+1],ptedges[i+2]);
+ Data_Levy_Ratio_e[i] = (ptedges[i+2]-ptedges[i+1])*spectrum_e[i]/myLevy->Integral(ptedges[i+1],ptedges[i+2]);
+ }
+ TGraphErrors *gr_Data_Levy_Ratio = new TGraphErrors(8, pt_points, Data_Levy_Ratio, pt_points_e, Data_Levy_Ratio_e);
+ gr_Data_Levy_Ratio->SetMarkerStyle(20);
+ gr_Data_Levy_Ratio->SetTitle("");
+ gr_Data_Levy_Ratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_Data_Levy_Ratio->GetYaxis()->SetTitle("(Measured Yield)/(Levy fit integral)");
+ //gr_Data_Levy_Ratio->Draw("AP");
+
+ // Reference Spectra
+ double LHC10b_XiStar[8]={0.00175909, 0.000893716, 0.000957646, 0.000399652, 0.000217706, 6.86327e-05, 3.24599e-05, 1.28201e-05};
+ double LHC10b_XiStar_e[8]={0.000460783, 0.000152526, 0.000136327, 6.06295e-05, 2.79773e-05, 1.2783e-05, 9.14534e-06, 4.461e-06};
+ double LHC10c_XiStar[8]={0.00153422, 0.00106629, 0.000752749, 0.000495777, 0.000255897, 6.88883e-05, 2.41039e-05, 1.09414e-05};
+ double LHC10c_XiStar_e[8]={0.000262538, 0.000106795, 6.6372e-05, 4.51931e-05, 2.07905e-05, 8.01244e-06, 4.03196e-06, 2.75332e-06};
+ double LHC10d_XiStar[8]={0.00137741, 0.00118096, 0.000813771, 0.000465871, 0.000223755, 6.5902e-05, 2.2662e-05, 7.11972e-06};
+ double LHC10d_XiStar_e[8]={7.81918e-05, 4.79133e-05, 2.92694e-05, 1.16203e-05, 5.00591e-06, 2.60709e-06, 1.18255e-06};
+ TGraphErrors *gr_10bXiStar = new TGraphErrors(8, pt_points, LHC10b_XiStar, pt_points_e, LHC10b_XiStar_e);
+ gr_10bXiStar->SetMarkerStyle(21);
+ gr_10bXiStar->SetMarkerColor(2);
+ gr_10bXiStar->SetLineColor(2);
+ TGraphErrors *gr_10cXiStar = new TGraphErrors(8, pt_points, LHC10c_XiStar, pt_points_e, LHC10c_XiStar_e);
+ gr_10cXiStar->SetMarkerStyle(22);
+ gr_10cXiStar->SetMarkerColor(4);
+ gr_10cXiStar->SetLineColor(4);
+ TGraphErrors *gr_10dXiStar = new TGraphErrors(8, pt_points, LHC10d_XiStar, pt_points_e, LHC10d_XiStar_e);
+ gr_10dXiStar->SetMarkerStyle(23);
+ gr_10dXiStar->SetMarkerColor(6);
+ gr_10dXiStar->SetLineColor(6);
+ //gr_10bXiStar->Draw("P");
+ //gr_10cXiStar->Draw("P");
+ //gr_10dXiStar->Draw("P");
+ //legend2->AddEntry(gr_10bXiStar,"LHC10b","p");
+ //legend2->AddEntry(gr_10cXiStar,"LHC10c","p");
+ //legend2->AddEntry(gr_10dXiStar,"LHC10d","p");
+
+ // Old QM2011 values for (Xi* + cc)/2
+ double Old_XiStarPoints[8]={0.001259, 0.0008586, 0.0006316, 0.0003604, 1.856e-4, 5.4601e-5, 2.0187e-5, 7.6559e-6};
+ double Old_XiStarPoints_e[8]={2.44e-4, 1.224e-4, 7.136e-5, 3.85e-5, 1.59e-5, 7.49e-6, 3.51e-6, 1.822e-6};
+ TGraphErrors *gr_OldXiStar = new TGraphErrors(8, pt_points, Old_XiStarPoints, pt_points_e, Old_XiStarPoints_e);
+ gr_OldXiStar->SetMarkerStyle(21);
+ gr_OldXiStar->SetMarkerColor(2);
+ gr_OldXiStar->SetLineColor(2);
+ //gr_OldXiStar->Draw("P");
+ //legend2->AddEntry(gr_OldXiStar,"(#Xi(1530) + cc)/2: QM11 Results","p");
+ /*
+ // To fit the Old QM results
+ TH1D *h_OldXiStar = (TH1D*)h_Xispectrum->Clone();
+ h_OldXiStar->Add(h_OldXiStar,-1);
+ for(int i=0; i<8; i++){
+ h_OldXiStar->SetBinContent(i+2, Old_XiStarPoints[i]);
+ h_OldXiStar->SetBinError(i+2, Old_XiStarPoints_e[i]);
+ }
+ h_OldXiStar->SetBinContent(1,+100);// first bin has no data
+ h_OldXiStar->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
+ cout<<"old mean pt = "<<myLevy->Moment(1,0.,10)<<endl;
+ */
+ //double pt_points_shifted[max_ptbins]={0};
+ //TF1 *AvgX = new TF1("AvgX","x*myExp",0.6,5.6);
+
+ //AvgX->SetParName(0,"A");
+ //AvgX->SetParName(1,"T");
+ //for(int i=0; i<8; i++){
+ //cout<<AvgX->Integral(ptedges[i],ptedges[i+1])/myExp->Integral(ptedges[i],ptedges[i+1])<<endl;
+ //}
+
+
+
+
+
+ double David_points_pt[18]={0.7, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.6, 1.8, 2.05, 2.4, 2.85, 3.5, 4.4, 5.45, 6.6, 7.85};
+
+
+ // David's new points
+ // David's Xi-
+ //double David_points_y[18]={0.00470846, 0.00494042, 0.00474707, 0.0046308, 0.00426987, 0.00390345, 0.00345295, 0.0032779, 0.00269152, 0.00202204, 0.00145725, 0.00088381, 0.000460706, 0.000196437, 6.31383e-05, 1.69709e-05, 5.72763e-06, 2.22189e-06};
+ //double David_points_y_e[18]={0.000313742, 0.000316252, 0.000294082, 0.000292834, 0.000255195, 0.000233124, 0.00020646, 0.000207611, 0.000160407, 0.000116738, 8.30453e-05, 5.03115e-05, 2.65867e-05, 1.14573e-05, 4.13348e-06, 1.4395e-06, 6.63351e-07, 4.01286e-07};
+ // David's Xi+
+ //double David_points_y[18]={0.00476135, 0.00480974, 0.00469503, 0.0046272, 0.0040248, 0.00364938, 0.00356068, 0.00306262, 0.0026221, 0.00196555, 0.00139214, 0.000856611, 0.000443817, 0.000192649, 6.06012e-05, 1.77073e-05, 5.04921e-06, 2.01156e-06};
+ //double David_points_y_e[18]={0.000324969, 0.000323376, 0.000320961, 0.00029042, 0.000248456, 0.000219039, 0.000231559, 0.000187812, 0.000153252, 0.000114861, 7.95886e-05, 4.92256e-05, 2.57394e-05, 1.15276e-05, 4.04731e-06, 1.43029e-06, 6.10334e-07, 3.50592e-07};
+ // Published ALICE points
+ // Xi-
+ double David_points_y[18]={4.734067e-03, 4.967289e-03, 4.772892e-03, 4.655989e-03, 4.293091e-03, 3.924680e-03, 3.471732e-03, 3.295728e-03, 2.706158e-03, 2.033039e-03, 1.465172e-03, 8.886171e-04, 4.632116e-04, 1.975056e-04, 6.348164e-05, 1.706320e-05, 5.758779e-06, 2.233978e-06};
+ double David_points_y_e[18]={3.203305e-04, 3.230481e-04, 3.006813e-04, 2.992160e-04, 2.610787e-04, 2.384269e-04, 2.111132e-04, 2.118811e-04, 1.638975e-04, 1.193448e-04, 8.491450e-05, 5.149418e-05, 2.718932e-05, 1.174048e-05, 4.233300e-06, 1.475196e-06, 6.953602e-07, 4.110991e-07};
+ // Xi+
+ //double David_points_y[18]={4.787242e-03, 4.835899e-03, 4.720567e-03, 4.652362e-03, 4.046691e-03, 3.669227e-03, 3.580047e-03, 3.079273e-03, 2.636364e-03, 1.976242e-03, 1.399711e-03, 8.612696e-04, 4.462308e-04, 1.936968e-04, 6.093081e-05, 1.780359e-05, 5.076670e-06, 2.022504e-06};
+ //double David_points_y_e[18]={3.315195e-04, 3.301755e-04, 3.272160e-04, 2.971086e-04, 2.538620e-04, 2.236666e-04, 2.363366e-04, 1.920731e-04, 1.566324e-04, 1.174641e-04, 8.139283e-05, 5.033151e-05, 2.632723e-05, 1.179350e-05, 4.145917e-06, 1.471487e-06, 6.243488e-07, 3.638050e-07};
+
+ TGraphErrors *gr_DavidXi = new TGraphErrors(18, David_points_pt, David_points_y, 0, David_points_y_e);
+ gr_DavidXi->SetMarkerStyle(24);
+ gr_DavidXi->SetMarkerSize(2);
+ gr_DavidXi->SetMarkerColor(4);
+ //gr_DavidXi->Draw("P");
+ //legend2->AddEntry(gr_DavidXi,"ALICE Published #Xi^{+}","p");
+ //legend2->AddEntry(gr_DavidXi,"#Xi^{+-}","p");
+
+
+
+ //cout<<"Efficiencies:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<Eff[i]<<", ";
+ //cout<<endl;
+ //cout<<"Corresponding errors:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<Eff_e[i]<<", ";
+ //cout<<endl;
+ //cout<<"Raw Yields:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<yield[i]<<", ";
+ //cout<<endl;
+ //cout<<"Ratio of My spectrum to Davids:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<spectrum[i]/David_points_y[i]<<", ";
+ //cout<<endl;
+ //cout<<"Corresponding ratio errors:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<sqrt( pow(spectrum_e[i]/David_points_y[i],2) + pow(spectrum[i]*David_points_y_e[i]/(David_points_y[i]*David_points_y[i]),2))<<", ";
+ //cout<<endl;
+ //cout<<"Ratio of My Eff to Davids:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<Eff[i]/David_XiMinus_eff[i]<<", ";
+ //cout<<endl;
+ //cout<<"Corresponding ratio errors:"<<endl;
+ //for(int i=0; i<max_ptbins; i++) cout<<sqrt( pow(Eff_e[i]/David_XiMinus_eff[i],2) + pow(Eff_e[i]*David_XiMinus_eff_e[i]/(David_XiMinus_eff[i]*David_XiMinus_eff[i]),2))<<", ";
+ //cout<<endl;
+
+
+ //Xi 10b+10c+10d with first eff bug fix, second bug fix is a very tiny change
+ //double Eff_ref[18]={0.0155246, 0.0366893, 0.0522083, 0.0653027, 0.0795129, 0.0942177, 0.109902, 0.125146, 0.150523, 0.184772, 0.225384, 0.27233, 0.324995, 0.371011, 0.382669, 0.354245, 0.344473, 0.299919};
+ // My eff with David's cuts XiMinus
+ //double Eff_refXiMinus[18]={0.0105995, 0.0312385, 0.0477496, 0.0645244, 0.0812671, 0.0996255, 0.117602, 0.13177, 0.158334, 0.19623, 0.230584, 0.280132, 0.324593, 0.351765, 0.353174, 0.321661, 0.285389, 0.20558};
+ //double Eff_refXiMinus_e[18]={0.000430506, 0.00115503, 0.00152918, 0.00189702, 0.00227815, 0.00268649, 0.0031259, 0.00357281, 0.00307359, 0.00400903, 0.00422553, 0.0052229, 0.00685152, 0.00835609, 0.0118039, 0.0168519, 0.0234869, 0.0274467};
+ //double Eff_refXiPlus[18]={0.0103308, 0.0291685, 0.0473052, 0.0586698, 0.0741367, 0.0936288, 0.110488, 0.129917, 0.15185, 0.192956, 0.226617, 0.286435, 0.318134, 0.355376, 0.358586, 0.325878, 0.249803, 0.256949};
+ //double Eff_refXiPlus_e[18]={0.000429664, 0.00112569, 0.00153199, 0.00180478, 0.00218061, 0.00261963, 0.00307867, 0.00358387, 0.00304106, 0.00398421, 0.00422447, 0.00531824, 0.00683356, 0.00842619, 0.012285, 0.0171912, 0.0224039, 0.0320204};
+
+ // 10d1 eff
+ double Eff_ref1[8]={0.00585138, 0.0271111, 0.0370528, 0.0767083, 0.0970953, 0.129362, 0.144039, 0.130414};
+ double Eff_ref1_e[8]={0.00110308, 0.00301515, 0.00433263, 0.00845994, 0.00972208, 0.0175748, 0.0318935, 0.0367763};
+ // 10d4 eff
+ double Eff_ref2[8]={0.00646719, 0.0236728, 0.0454735, 0.0692855, 0.0922256, 0.13542, 0.161154, 0.126502};
+ double Eff_ref2_e[8]={0.000740606, 0.00176465, 0.00315074, 0.00498191, 0.00600009, 0.0117801, 0.0203019, 0.0247307};
+ // 10f6a eff
+ double Eff_ref3[8]={0.00589799, 0.0228898, 0.043247, 0.0692314, 0.104212, 0.14611, 0.1587, 0.176539};
+ double Eff_ref3_e[8]={0.000465501, 0.00114491, 0.00199253, 0.00333088, 0.00421532, 0.00822058, 0.0133099, 0.0202604};
+ // 10f6a eff doubled binning
+ //double Eff_ref3[16]={0.00288254, 0.00957083, 0.0177512, 0.0293962, 0.0382567, 0.0496768, 0.0650758, 0.074635, 0.096423, 0.116841, 0.142147, 0.15214, 0.169098, 0.143671, 0.185977, 0.162155};
+ //double Eff_ref3_e[16]={0.000438375, 0.000885026, 0.00134525, 0.00195975, 0.00249222, 0.00323944, 0.00428832, 0.00525628, 0.00513653, 0.00726998, 0.0104245, 0.0133468, 0.0179437, 0.0196858, 0.0268477, 0.030696};
+
+ // Pythia to Phojet eff ratio
+ for(int ii=0; ii<max_ptbins; ii++){
+ Eff_e[ii] = sqrt(pow(Eff_e[ii]/Eff_ref3[ii],2) + pow(Eff[ii]*Eff_ref3_e[ii]/pow(Eff_ref3[ii],2),2));
+ Eff[ii] /= Eff_ref3[ii];
+ }
+
+ TGraphErrors *gr_Eff = new TGraphErrors(max_ptbins,pt_points, Eff, pt_points_e, Eff_e);
+ gr_Eff->SetMarkerStyle(20);
+ gr_Eff->SetMinimum(0);
+ gr_Eff->SetMaximum(.4);
+ //gr_Eff->GetXaxis()->SetLimits(0,6.1);
+ gr_Eff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_Eff->GetYaxis()->SetTitle("Efficiency");
+ //gr_Eff->SetTitle("#Xi(1530)^{0}+cc Efficiency");
+ //gr_Eff->Draw("AP");
+ //legend2->AddEntry(gr_Eff,"10d1+10d4+10f6a","p");
+ //legend2->AddEntry(gr_Eff,"10f6 (Phojet)","p");
+ //gr_Eff->Fit("pol0","IME","",1.6,4.8);
+
+ TGraphErrors *gr_Eff_ref1 = new TGraphErrors(max_ptbins,pt_points, Eff_ref1, pt_points_e, Eff_ref1_e);
+ gr_Eff_ref1->SetMarkerStyle(20);
+ gr_Eff_ref1->SetMarkerColor(2);
+ gr_Eff_ref1->SetLineColor(2);
+ //gr_Eff_ref1->Draw("P");
+ //legend2->AddEntry(gr_Eff_ref1,"10d1 MC","p");
+ TGraphErrors *gr_Eff_ref2 = new TGraphErrors(max_ptbins,pt_points, Eff_ref2, pt_points_e, Eff_ref2_e);
+ gr_Eff_ref2->SetMarkerStyle(20);
+ gr_Eff_ref2->SetMarkerColor(4);
+ gr_Eff_ref2->SetLineColor(4);
+ //gr_Eff_ref2->Draw("P");
+ //legend2->AddEntry(gr_Eff_ref2,"10d4 MC","p");
+ TGraphErrors *gr_Eff_ref3 = new TGraphErrors(max_ptbins,pt_points, Eff_ref3, pt_points_e, Eff_ref3_e);
+ gr_Eff_ref3->SetMarkerStyle(20);
+ gr_Eff_ref3->SetMarkerColor(3);
+ gr_Eff_ref3->SetLineColor(3);
+ //gr_Eff_ref3->Draw("P");
+ //legend2->AddEntry(gr_Eff_ref3,"10f6a (Pythia Perugia)","p");
+
+
+ //TGraphErrors *gr_Eff_xiplus = new TGraphErrors(max_ptbins,pt_points, Eff_refXiPlus, 0, Eff_refXiPlus_e);
+ //gr_Eff_xiplus->SetMarkerStyle(21);
+ //gr_Eff_xiplus->SetMarkerColor(4);
+ //gr_Eff_xiplus->GetXaxis()->SetTitle("pt (GeV/c)");
+ //gr_Eff_xiplus->GetYaxis()->SetTitle("Efficiency");
+ //gr_Eff_xiplus->SetTitle("#Xi(1530)^0 Efficiency");
+ //gr_Eff_xiplus->Draw("P");
+
+
+
+ TGraphErrors *gr_David_Eff_XiMinus = new TGraphErrors(18,pt_points_David,David_XiMinus_eff,0,David_XiMinus_eff_e);
+ gr_David_Eff_XiMinus->SetMarkerStyle(24);
+ gr_David_Eff_XiMinus->SetMarkerColor(2);
+ //gr_David_Eff_XiMinus->Draw("P");
+ //legend2->AddEntry(gr_David_Eff_XiMinus,"David/Antonin Xi- Eff","p");
+
+ TGraphErrors *gr_David_Eff_XiPlus = new TGraphErrors(18,pt_points_David,David_XiPlus_eff,0,David_XiPlus_eff_e);
+ gr_David_Eff_XiPlus->SetMarkerStyle(25);
+ gr_David_Eff_XiPlus->SetMarkerColor(4);
+ //gr_David_Eff_XiPlus->Draw("P");
+
+
+ legend2->Draw("same");
+
+
+ TLegend *legend3 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
+ legend3->SetBorderSize(1);
+ legend3->SetTextSize(.04);// small .03; large .036
+ //legend3->SetLineColor(0);
+ legend3->SetFillColor(0);
+
+
+ TGraphErrors *gr_RealWidth = new TGraphErrors(8, pt_points, width, pt_points_e, width_e);
+ gr_RealWidth->SetMarkerStyle(20);
+ gr_RealWidth->SetMinimum(0);
+ gr_RealWidth->SetMaximum(0.0035);
+ gr_RealWidth->GetYaxis()->SetTitle("#sigma (GeV/c)");
+ gr_RealWidth->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_RealWidth->SetTitle("");
+ //gr_RealWidth->Draw("AP");
+ //legend3->AddEntry(gr_RealWidth,"Real data (10b+10c+10d)","p");
+
+
+ TGraphErrors *gr_MCWidth = new TGraphErrors(8, pt_points, MCwidth, pt_points_e, MCwidth_e);
+ gr_MCWidth->SetMarkerStyle(20);
+ gr_MCWidth->SetMarkerColor(2);
+ gr_MCWidth->SetLineColor(2);
+ gr_MCWidth->GetYaxis()->SetTitle("#sigma (GeV/c)");
+ gr_MCWidth->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_MCWidth->SetTitle("");
+ //gr_MCWidth->Draw("P");
+ //legend3->AddEntry(gr_MCWidth,"Pythia+ALICE (10d1+10d4+10f6a)","p");
+
+
+
+
+ TGraphErrors *gr_RealMass = new TGraphErrors(8, pt_points, mass, pt_points_e, mass_e);
+ gr_RealMass->SetMarkerStyle(20);
+ gr_RealMass->SetMinimum(1.529);
+ gr_RealMass->SetMaximum(1.5345);
+ gr_RealMass->SetTitle("");
+ gr_RealMass->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ gr_RealMass->GetYaxis()->SetTitleOffset(1.3);
+ gr_RealMass->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
+ //gr_RealMass->Draw("AP");
+ //legend3->AddEntry(gr_RealMass,"Real data (10b+10c+10d)","p");
+
+ TGraphErrors *gr_MCMass = new TGraphErrors(8, pt_points, MCmass, pt_points_e, MCmass_e);
+ gr_MCMass->SetMarkerStyle(20);
+ gr_MCMass->SetMarkerColor(2);
+ gr_MCMass->SetTitle("");
+ gr_MCMass->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
+ gr_MCMass->GetYaxis()->SetTitleOffset(1.3);
+ gr_MCMass->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ //gr_MCMass->Draw("P");
+ //legend3->AddEntry(gr_MCMass,"Pythia+ALICE (10d1+10d4+10f6a)","p");
+ //TF1 *pdg_massline=new TF1("pdg_massline","1.53178",0,10);
+ //pdg_massline->Draw("same");
+
+ //legend3->Draw("same");
+ TF1 *myLevy2=new TF1("myLevy2",myLevyPt,0,8.5,3);
+ myLevy2->SetLineColor(4);
+ myLevy2->SetParName(0,"dN/dy");
+ myLevy2->SetParName(1,"C");
+ myLevy2->SetParName(2,"n");
+ myLevy2->SetParameter(0,.003);
+ myLevy2->SetParameter(1,.3);
+ myLevy2->SetParameter(2,15);
+ myLevy2->SetParLimits(0,0.0001,0.01);
+ myLevy2->SetParLimits(1,.1,1);
+ myLevy2->SetParLimits(2,1,500);
+ TF1 *myExp2 = new TF1("myExp2",myExpFit,0,8.5,2);
+ myExp2->SetLineColor(4);
+ myExp2->SetParName(0,"dN/dy");
+ myExp2->SetParName(1,"C");
+ myExp2->SetParameter(0,.003);
+ myExp2->SetParameter(1,.3);
+ myExp2->SetParLimits(0,0.0001,0.1);
+ myExp2->SetParLimits(1,.1,1);
+
+
+ MCinput_Spectrum->SetMarkerStyle(20);
+ MCinput_Spectrum->SetMarkerColor(2);
+ MCinput_Spectrum->SetLineColor(2);
+ MCinput_Spectrum->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+ MCinput_Spectrum->GetYaxis()->SetTitle("1/N_{E}d^{2}N/dydp_{T} (GeV/c)^{-1}");
+ //MCinput_Spectrum->SetTitle("Pythia Perugia");
+ //MCinput_Spectrum->SetTitle("Phojet");
+
+ MCinput_Spectrum->Sumw2();
+ MCinput_Spectrum->Scale(1/MCinput_Spectrum->GetXaxis()->GetBinWidth(1));
+ MCinput_Spectrum->Scale(1/EventsMC_postPV->GetEntries());
+
+ // Scale Pythia to Data fit Tsallis
+ MCinput_Spectrum->Scale(myLevy->Integral(0.8,1.2)/(MCinput_Spectrum->Integral(9,12)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)));
+ //
+ cout<<"low-pt extrapolation from Levy = "<<myLevy->Integral(0,0.8)<<" Pythia Histogram Integral = "<<MCinput_Spectrum->Integral(1,8)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
+ cout<<"high pt extrapolation over total Levy = "<<myLevy->Integral(5.6,10.)/myLevy->Integral(0,10.)<<endl;
+ cout<<"high pt Pythia over total Pythia = "<<MCinput_Spectrum->Integral(56,100)/MCinput_Spectrum->Integral(1,100)<<endl;
+ cout<<"Total yield Pythia over total Levy = "<<MCinput_Spectrum->Integral(1,100)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)/myLevy->Integral(0,10.)<<endl;
+
+ //MCinput_Spectrum->Draw();
+ //minFitpoint=0.; maxFitpoint=10.;
+ //MCinput_Spectrum->Fit(myLevy2,"IME","",minFitpoint,maxFitpoint);
+ //MCinput_Spectrum->Fit(myExp,"IME","",minFitpoint,maxFitpoint);
+ //myExp->Draw("same");
+ //myLevy2->Draw("same");
+ //myLevy->Draw("same");
+ //legend3->AddEntry(MCinput_Spectrum,"Pythia Perugia (Scaled to Real Data)","p");
+ //legend3->AddEntry(myLevy2,"Phojet Levy fit (fit range from pt=0.8 to 5.6 GeV/c","l");
+ //legend3->AddEntry(myLevy,"Real data Levy fit","l");
+
+ double xaxis[60]={0};
+ double yaxis[60]={0};
+ for(int ii=0; ii<60; ii++){
+ xaxis[ii] = MCinput_Spectrum->GetXaxis()->GetBinCenter(ii+1);
+ yaxis[ii] = myLevy->Eval(xaxis[ii])/myLevy2->Eval(xaxis[ii]);
+ }
+ TGraph *LevyRatio = new TGraph(60, xaxis, yaxis);
+ LevyRatio->SetTitle("Levy Ratio");
+ LevyRatio->GetXaxis()->SetTitle("p_{T}");
+ LevyRatio->GetYaxis()->SetTitle("(Real Data Levy)/(Pythia fit Levy)");
+ LevyRatio->SetMarkerStyle(20);
+ //LevyRatio->Draw("AP");
+
+
+ //legend3->Draw("same");
+
+ //cout<<"low-pt extrapolation from Levy = "<<myLevy2->Integral(0,0.8)<<" Histogram Integral = "<<MCinput_Spectrum->Integral(1,8)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
+ //cout<<"Total Bin Count yield = "<<MCinput_Spectrum->Integral(1,MCinput_Spectrum->GetNbinsX())*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
+ //cout<<"low-pt extrapolation from Exp = "<<myExp->Integral(0,0.8)<<" Histogram Integral = "<<MCinput_Spectrum->Integral(1,8)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
+
+
+ if(SaveToFile){
+ TFile *outputfile = new TFile(outfilename->Data(),"RECREATE");
+ can->Write("can");
+ can2->Write("can2");
+ myLevy->Write("myLevy");
+ h_Xispectrum->Write("h_Spectrum");
+ RawYields->Write("RawYields");
+ Efficiency->Write("Efficiency");
+
+ outputfile->Close();
+ }
+
+ cout<<endl;
+
+}
+
+//________________________________________________________________________
+double PolFunction(double *x, double *par){
+
+ if (reject && x[0] > yieldRange[0] && x[0] < yieldRange[1]) {
+ TF1::RejectPoint();
+ return 0;
+ }
+
+ if(POLdegree==1) return par[0] + par[1]*x[0];
+ else if(POLdegree==2) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2);
+ else if(POLdegree==3) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3);
+ else if(POLdegree==4) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3) + par[4]*pow(x[0],4);
+ else if(POLdegree==5) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3) + par[4]*pow(x[0],4) + par[5]*pow(x[0],5);
+ else return 0;
+
+}
+//________________________________________________________________________
+double PolFunctionSpecialXi(double *x, double *par){
+
+ if (x[0] > 1.28 && x[0] < 1.3) {
+ TF1::RejectPoint();
+ return 0;
+ }
+ if (reject && x[0] > 1.31 && x[0] < 1.334) {
+ TF1::RejectPoint();
+ return 0;
+ }
+
+ if(POLdegree==1) return par[0] + par[1]*x[0];
+ else if(POLdegree==2) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2);
+ else if(POLdegree==3) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3);
+ else if(POLdegree==4) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3) + par[4]*pow(x[0],4);
+ else if(POLdegree==5) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3) + par[4]*pow(x[0],4) + par[5]*pow(x[0],5);
+ else return 0;
+
+}
+//________________________________________________________________________
+double BWFunction(double *x, double *par){
+ return (par[0]*par[2]/(1000.*2.*3.1415926))/( pow(x[0]-par[1],2) + par[2]*par[2]/4.);
+ //return (par[0]*par[1]/(1000.*2.*3.1415926))/( pow(x[0]-par[2],2) + par[1]*par[1]/4.);
+}
+//________________________________________________________________________
+double BWplusPol(double *x, double *par){
+ return BWFunction(x,par) + PolFunction(x,&par[3]);
+}
+//________________________________________________________________________
+double GausFunction(double *x, double *par){
+ return (par[0]*.001)/(sqrt(2*3.14159*par[2]*par[2]))*exp(-pow((x[0]-par[1])/(sqrt(2)*par[2]),2));
+}
+//________________________________________________________________________
+double GausplusPol(double *x, double *par){
+ return GausFunction(x,par) + PolFunction(x,&par[3]);
+}
+//________________________________________________________________________
+double GausplusPolSpecialXi(double *x, double *par){
+ return GausFunction(x,par) + PolFunctionSpecialXi(x,&par[3]);
+}
+//________________________________________________________________________
+double myLevyPt(Double_t *x, Double_t *par)
+{
+ Double_t lMass=0;
+ if(XiStarCase) lMass = 1.5318; //Xi mass
+ else lMass = 1.32171; //Xi mass
+
+ Double_t ldNdy = par[0]; // dN/dy
+ Double_t l2pi = 2*TMath::Pi(); // 2pi
+ Double_t lTemp = par[1]; // Temperature
+ Double_t lPower = par[2]; // power=n
+
+ Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
+ Double_t lInPower = 1 + (TMath::Sqrt(x[0]*x[0]+lMass*lMass)-lMass) / (lPower*lTemp);
+
+ return l2pi * ldNdy * x[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
+}
+//________________________________________________________________________
+double myExpFit(Double_t *x, Double_t *par)
+{
+ Double_t lMass=0;
+ if(XiStarCase) lMass = 1.5318; //Xi mass
+ else lMass = 1.32171; //Xi mass
+
+ return 2*TMath::Pi()*x[0]*par[0]*exp(-TMath::Sqrt(x[0]*x[0]+lMass*lMass)/par[1]);
+}
+//________________________________________________________________________
+double Voigt(Double_t *x, Double_t *par)
+{// code taken directly from RooVoigtian in ROOT
+
+ Double_t Norm = par[0];
+ Double_t mean = par[1];
+ Double_t s = fabs(par[2]);// sigma; Gaussian sigma
+ Double_t w = fabs(par[3]);// Width; BW width
+
+ Double_t arg = x[0] - mean;
+ Double_t _invRootPi = 1./sqrt(atan2(0.,-1.));
+
+ // return constant for zero width and sigma
+ if (s==0. && w==0.) return 1.;
+
+ // Breit-Wigner for zero sigma
+ if (s==0.) return (Norm*1./(arg*arg+0.25*w*w));
+ Double_t coef= -0.5/(s*s);
+
+ // Gauss for zero width
+ if (w==0.) return Norm*exp(coef*arg*arg);
+
+ // actual Voigtian for non-trivial width and sigma
+ Double_t c = 1./(sqrt(2.)*s);
+ Double_t a = 0.5*c*w;
+ Double_t u = c*arg;
+ RooComplex z(u,a);
+ RooComplex v(0.);
+
+
+ v = RooMath::ComplexErrFunc(z);
+
+ return Norm*c*_invRootPi*v.re();
+
+}
+//________________________________________________________________________
+ double VoigtplusPol(double *x, double *par){
+ return Voigt(x,par) + PolFunction(x,&par[4]);
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff