PWGLF/RESONANCES/extra/Plot_macroXiStar.C

   1 //gROOT->Reset();
   2 #include <math.h>
   3 #include <time.h>
   4 #include <stdio.h>
   5 #include <stdlib.h>
   6 #include <Riostream.h>
   7 #include "RooFit.h"
   8 #include "RooVoigtian.h"
   9 #include "RooAbsReal.h"
  10 #include "RooRealVar.h"
  11 #include "RooComplex.h"
  12 #include "RooMath.h"
  13
  14 #include "TVector2.h"
  15 #include "TFile.h"
  16 #include "TString.h"
  17 #include "TF1.h"
  18 #include "TH1.h"
  19 #include "TH2.h"
  20 #include "TH3.h"
  21 #include "TProfile.h"
  22 #include "TProfile2D.h"
  23 #include "TMath.h"
  24 #include "TText.h"
  25 #include "TRandom3.h"
  26 #include "TArray.h"
  27 #include "TLegend.h"
  28 #include "TStyle.h"
  29 #include "TMinuit.h"
  30
  31 #include "TCanvas.h"
  32 #include "TLatex.h"
  33 #include "TImage.h"
  34 #include "TGaxis.h"
  35 #include "TGraphErrors.h"
  36 #include "TPaveStats.h"
  37
  38 using namespace std;
  39
  40 bool SaveToFile=kFALSE;
  41 int VariationType=0;// 0 = CutVariations, 1 = FitVariations
  42 //
  43 int CutVariation=0;// 0=Standard Cuts, 1-12 are Systematic deviations
  44 int FitVariation=1;// 1 for VoigtOnly with BkgSub, 2 for VoigtOnly without BkgSub, 3 for RightNorm, 4 Alternate BinCounting Range
  45 TString *outfilename = new TString();
  46 //
  47 // below are default (standard) settings
  48 int POLdegree=1;// 1 with BkgSub, 5 without
  49 bool XiStarCase=kTRUE;// Standard: kTRUE
  50 bool VoigtFitOnly=kFALSE;// Standard: kFALSE
  51 bool bkgSubtract=kTRUE;// Standard: kTRUE
  52 bool LeftNorm=kTRUE;// Standard: kTRUE
  53 bool MCAssocCase=kTRUE;// Standard: kTRUE
  54 short ParticleCase=0;// 0 for particle+anti-particle. 1 for particle. 2 for anti-particle
  55
  56
  57 const int max_ptbins = 8; // 18 for Xi, 8 for XiStar
  58 double raprange[2]={-0.5,0.5};// for both Xi and XiStar (bin size = 0.1)
  59
  60 double XiyieldRange[2]={1.313, 1.332};
  61 double FitRangeXi[2]={1.27,1.36};
  62
  63 double XiStaryieldRange[2]={1.51, 1.56};
  64 double FitRangeXiStar[2]={1.48,1.59};
  65
  66 double XiCountBound[4]={1.304,1.313,1.331,1.340};// -9sigma, -4.5sigma, +4.5sigma, +9sigma
  67 double XiStarCountBound[6]={1.481, 1.490, 1.523, 1.541, 1.570, 1.579};// 1 gamma range (Standard)
  68 double XiStarCountBound_SysRange[6]={1.481, 1.490, 1.528, 1.537, 1.570, 1.579};// 1/2*gamma range (Systematic Variation)
  69
  70 double yieldRange[2]={0};
  71 bool reject;
  72
  73
  74 double PolFunction(double *x, double *par);
  75 double PolFunctionSpecialXi(double *x, double *par);
  76 double BWFunction(double *x, double *par);
  77 double BWplusPol(double *x, double *par);
  78 double GausFunction(double *x, double *par);
  79 double GausplusPol(double *x, double *par);
  80 double GausplusPolSpecialXi(double *x, double *par);
  81 double myLevyPt(Double_t *x, Double_t *par);
  82 double myExpFit(Double_t *x, Double_t *par);
  83 double Voigt(Double_t *x, Double_t *par);
  84 double VoigtplusPol(double *x, double *par);
  85
  86
  87
  88 void Plot_macroXiStar(){
  89
  90   if(VariationType==0) {
  91     outfilename->Append("CutVariation_");
  92     *outfilename += CutVariation;
  93   }
  94   if(VariationType==1) {
  95     outfilename->Append("FitVariation_");
  96     *outfilename += FitVariation;
  97     // 1 for VoigtOnly with BkgSub, 2 for VoigtOnly without BkgSub, 3 for RightNorm, 4 Alternate BinCounting Range
  98     if(FitVariation==1) {VoigtFitOnly=kTRUE;}
  99     if(FitVariation==2) {VoigtFitOnly=kTRUE; bkgSubtract=kFALSE; POLdegree=5;}
 100     if(FitVariation==3) {LeftNorm=kFALSE;}
 101     if(FitVariation==4) {for(int ii=0; ii<6; ii++) XiStarCountBound[ii]=XiStarCountBound_SysRange[ii];}
 102   }
 103   outfilename->Append(".root");
 104
 105   double FitRange[2]={0};
 106   if(XiStarCase) {
 107     FitRange[0] = FitRangeXiStar[0]; FitRange[1] = FitRangeXiStar[1];
 108     yieldRange[0] = XiStaryieldRange[0]; yieldRange[1] = XiStaryieldRange[1];
 109   }
 110   else {
 111     FitRange[0] = FitRangeXi[0]; FitRange[1] = FitRangeXi[1];
 112     yieldRange[0] = XiyieldRange[0]; yieldRange[1] = XiyieldRange[1];
 113   }
 114
 115
 116
 117   TGaxis::SetMaxDigits(3);
 118   gStyle->SetOptStat(0);
 119   gStyle->SetOptFit(111);
 120   //gStyle->SetOptFit(1111);
 121
 122
 123   TCanvas *can = new TCanvas("can","can",13,34,700,500);
 124   gStyle->SetOptFit(1111);
 125   can->Range(-1.25,-0.2625,11.25,2.3625);
 126   can->SetFillColor(10);
 127   can->SetBorderMode(0);
 128   can->SetBorderSize(2);
 129   //can->SetGridx();
 130   //can->SetGridy();
 131   can->SetFrameFillColor(0);
 132   can->SetFrameBorderMode(0);
 133   can->SetFrameBorderMode(0);
 134   can->Divide(4,2);
 135
 136   TLegend *legend = new TLegend(.1,.7,.35,.9,NULL,"brNDC");
 137   legend->SetBorderSize(1);
 138   legend->SetTextSize(.04);// small .03; large .036
 139   //legend->SetLineColor(0);
 140   legend->SetFillColor(0);
 141   TLegend *legend2 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
 142   legend2->SetBorderSize(1);
 143   legend2->SetTextSize(.04);// small .03; large .036
 144   //legend2->SetLineColor(0);
 145   legend2->SetFillColor(0);
 146
 147
 148
 149   int PARBINS_temp;
 150   PARBINS_temp = POLdegree+1+4;
 151   const int PARBINS=PARBINS_temp;
 152   int offset=0;
 153   offset=4;
 154
 155   double bkg_params[POLdegree+1];
 156
 157   TString *polname=new TString();
 158   if(POLdegree==1) polname->Append("pol1(0)");
 159   if(POLdegree==2) polname->Append("pol2(0)");
 160   if(POLdegree==3) polname->Append("pol3(0)");
 161   if(POLdegree==4) polname->Append("pol4(0)");
 162   if(POLdegree==5) polname->Append("pol5(0)");
 163
 164
 165   TF1 *myPol[max_ptbins];
 166   for(int i=0; i<max_ptbins; i++){
 167     TString *name = new TString("myPol");
 168     *name +=i+1;
 169     myPol[i] = new TF1(name->Data(),polname->Data(),0,2);
 170     myPol[i]->SetLineColor(3);
 171     myPol[i]->SetRange(FitRange[0],FitRange[1]);
 172   }
 173
 174
 175   //
 176
 177
 178   // use 10b + 10c + 10d for final result
 179   //TFile *myfile = new TFile("Results/Real_LHC10b_lego.root","READ");
 180   //TFile *myfile = new TFile("Results/Real_LHC10c_lego.root","READ");
 181   //TFile *myfile = new TFile("Results/Real_LHC10d_lego.root","READ");
 182   TFile *myfile = new TFile("Results/RealMerge_10b_10c_10d_lego.root","READ");// accepted version
 183
 184
 185   // use 10d1 + 10d4 + 10f6a for final result
 186   //TFile *myfile2 = new TFile("Results/MC_LHC10d1_lego.root","READ");
 187   //TFile *myfile2 = new TFile("Results/MC_LHC10d4_lego.root","READ");
 188   //TFile *myfile2 = new TFile("Results/MC_LHC10f6a_lego.root","READ");// Pythia
 189   TFile *myfile2 = new TFile("Results/MCMerge_10d1_10d4_10f6a_lego.root","READ");// accepted version
 190   //TFile *myfile2 = new TFile("Results/MC_LHC10f6_lego.root","READ");// Phojet
 191
 192   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};
 193   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};
 194   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};
 195   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};
 196
 197   TDirectoryFile *tdir=(TDirectoryFile*)myfile->Get("PWGLF.outputXiStarAnalysis.root");
 198   TList *List1=(TList*)tdir->Get("XiStarOutput");
 199   //TList *List1=(TList*)myfile->Get("MyList");
 200   myfile->Close();
 201
 202   TDirectoryFile *tdir2=(TDirectoryFile*)myfile2->Get("PWGLF.outputXiStarAnalysis.root");
 203   TList *List2=(TList*)tdir2->Get("XiStarOutput");
 204   //TList *List2=(TList*)myfile2->Get("MyList");
 205   myfile2->Close();
 206
 207   TH1F *Events_PhysSel=(TH1F*)(List1->FindObject("fMultDist1"));
 208   TH1F *Events_postPV=(TH1F*)(List1->FindObject("fMultDist3"));
 209   TH1F *EventsMC_PhysSel=(TH1F*)(List2->FindObject("fMultDist1"));
 210   TH1F *EventsMC_postPV=(TH1F*)(List2->FindObject("fMultDist3"));
 211
 212   cout<<"# Events passing PhysSelection Cuts = "<<Events_PhysSel->GetEntries()<<endl;
 213   cout<<"# Events passing PhysSel + Vz Cut + PileUp Cut = "<<Events_postPV->GetEntries()<<endl;
 214   cout<<"# MC Events passing Zvertex and PhysSelection Cuts = "<<EventsMC_postPV->GetEntries()<<endl;
 215
 216   TString *name = new TString("fXi_0");
 217   TH3F *Xi3 = (TH3F*)(List1->FindObject(name->Data()));
 218   name = new TString("fXibar_0");
 219   TH3F *Xi_bar3 = (TH3F*)(List1->FindObject(name->Data()));
 220
 221   name = new TString("fXiMinusPiPlus_"); *name += CutVariation;
 222   TH3F *XiMinus_piPlus3 = (TH3F*)(List1->FindObject(name->Data()));
 223   name = new TString("fXiMinusPiMinus_"); *name += CutVariation;
 224   TH3F *XiMinus_piMinus3 = (TH3F*)(List1->FindObject(name->Data()));
 225   name = new TString("fXiPlusPiPlus_"); *name += CutVariation;
 226   TH3F *XiPlus_piPlus3 = (TH3F*)(List1->FindObject(name->Data()));
 227   name = new TString("fXiPlusPiMinus_"); *name += CutVariation;
 228   TH3F *XiPlus_piMinus3 = (TH3F*)(List1->FindObject(name->Data()));
 229
 230   name = new TString("fXiMinusPiPlusbkg_"); *name += CutVariation;
 231   TH3F *XiMinus_piPlus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
 232   name = new TString("fXiMinusPiMinusbkg_"); *name += CutVariation;
 233   TH3F *XiMinus_piMinus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
 234   name = new TString("fXiPlusPiPlusbkg_"); *name += CutVariation;
 235   TH3F *XiPlus_piPlus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
 236   name = new TString("fXiPlusPiMinusbkg_"); *name += CutVariation;
 237   TH3F *XiPlus_piMinus3_bkg = (TH3F*)(List1->FindObject(name->Data()));
 238
 239
 240   TH3F *MCinput_Xi_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXi1"));
 241   TH3F *MCinput_Xi_bar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXibar1"));
 242   TH3F *MCinput_XiStar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXiStar1"));
 243   TH3F *MCinput_XiStar_bar_prePV = (TH3F*)(List2->FindObject("fMCinputTotalXiStarbar1"));
 244   TH3F *MCinput_XiStar = (TH3F*)(List2->FindObject("fMCinputTotalXiStar3"));
 245   TH3F *MCinput_XiStar_bar = (TH3F*)(List2->FindObject("fMCinputTotalXiStarbar3"));
 246   TH3F *MCinput_Xi = (TH3F*)(List2->FindObject("fMCinputTotalXi3"));
 247   TH3F *MCinput_Xi_bar = (TH3F*)(List2->FindObject("fMCinputTotalXibar3"));
 248
 249   int binYL = MCinput_XiStar->GetYaxis()->FindBin(raprange[0]+0.01);
 250   int binYH = MCinput_XiStar->GetYaxis()->FindBin(raprange[1]-0.01);
 251   int binML = MCinput_XiStar->GetZaxis()->FindBin(1.51);
 252   int binMH = MCinput_XiStar->GetZaxis()->FindBin(1.55);
 253
 254   TH1D *MCinput_Spectrum;
 255   TH1D *MCinput_Spectrum_bar;
 256   if(XiStarCase){
 257     MCinput_Spectrum= (TH1D*)MCinput_XiStar->ProjectionX("MCinput_Spectrum",binYL,binYH, binML,binMH);
 258     MCinput_Spectrum_bar = (TH1D*)MCinput_XiStar_bar->ProjectionX("MCinput_Spectrum_bar",binYL,binYH, binML,binMH);
 259   }else{
 260     MCinput_Spectrum= (TH1D*)MCinput_Xi->ProjectionX("MCinput_Spectrum",binYL,binYH, binML,binMH);
 261     MCinput_Spectrum_bar = (TH1D*)MCinput_Xi_bar->ProjectionX("MCinput_Spectrum_bar",binYL,binYH, binML,binMH);
 262   }
 263   MCinput_Spectrum->Add(MCinput_Spectrum_bar);
 264
 265   name = new TString("fMCrecXiMinusPiPlus_"); *name += CutVariation;
 266   TH3F *MCrec_XiStar = (TH3F*)(List2->FindObject(name->Data()));
 267   name = new TString("fMCrecXiPlusPiMinus_"); *name += CutVariation;
 268   TH3F *MCrec_XiStar_bar = (TH3F*)(List2->FindObject(name->Data()));
 269
 270
 271   TH3F *MCrec_Xi;
 272   TH3F *MCrec_Xi_bar;
 273   if(MCAssocCase){
 274     MCrec_Xi = (TH3F*)(List2->FindObject("fMCrecXi"));
 275     MCrec_Xi_bar = (TH3F*)(List2->FindObject("fMCrecXibar"));
 276   }else{
 277     name = new TString("fXi_0");
 278     MCrec_Xi = (TH3F*)(List2->FindObject(name->Data()));
 279     name = new TString("fXibar_0");
 280     MCrec_Xi_bar = (TH3F*)(List2->FindObject(name->Data()));
 281   }
 282
 283
 284
 285
 286   double ptedges_Dhevan[9]={0.8, 1.2, 1.6, 2.0, 2.4, 3.2, 4.0, 4.8, 5.6};// standard
 287   double pt_points_Dhevan[8]={1.0, 1.4, 1.8, 2.2, 2.8, 3.6, 4.4, 5.2};// standard
 288   double pt_points_e_Dhevan[8]={.2, .2, .2, .2, .4, .4, .4, .4};// standard
 289   //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};
 290   //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};
 291   //double pt_points_e_Dhevan[16]={.1, .1, .1, .1, .1, .1, .1, .1, .2, .2, .2, .2, .2, .2, .2, .2};
 292
 293   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};
 294   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};
 295   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};
 296
 297
 298   double ptedges[max_ptbins+1+1]={0};
 299   double pt_points[max_ptbins]={0};
 300   double pt_points_e[max_ptbins]={0};
 301   TString *names[max_ptbins];
 302
 303   for(int i=0; i<max_ptbins; i++){
 304     if(XiStarCase){
 305       //ptedges[i] = ptedges_Dhevan[i];
 306       //ptedges[max_ptbins] = ptedges_Dhevan[max_ptbins];
 307       pt_points[i] = pt_points_Dhevan[i];
 308       pt_points_e[i] = pt_points_e_Dhevan[i];
 309       names[i] = new TString("pt bin ");
 310       *names[i] += i+1;
 311     }
 312     else{
 313       //ptedges[i] = ptedges_David[i];
 314       //ptedges[max_ptbins] = ptedges_David[max_ptbins];
 315       pt_points[i] = pt_points_David[i];
 316       pt_points_e[i] = pt_points_e_David[i]/2.;
 317       names[i] = new TString();
 318       *names[i] += float(ptedges_David[i]);
 319       names[i]->Append(" < pt < ");
 320       *names[i] += float(ptedges_David[i+1]);
 321     }
 322   }
 323   for(int i=0; i<max_ptbins+1; i++){
 324     if(XiStarCase){
 325       if(i==0) ptedges[i]=0;
 326       else {
 327         ptedges[i] = ptedges_Dhevan[i-1];
 328         ptedges[max_ptbins+1] = ptedges_Dhevan[max_ptbins];
 329       }
 330     }else {
 331       if(i==0) ptedges[i]=0;
 332       else {
 333         ptedges[i] = ptedges_David[i-1];
 334         ptedges[max_ptbins+1] = ptedges_David[max_ptbins];
 335       }
 336     }
 337   }
 338
 339   double Eff_corr_Xiplus[max_ptbins]={0};
 340   double Eff_corr_Ximinus[max_ptbins]={0};
 341
 342   if(!XiStarCase){
 343     // David binning
 344     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};
 345     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};
 346     for(int i=0; i<max_ptbins; i++) {
 347       Eff_corr_Xiplus[i] = geant3fluka_corr_Xiplus[i];
 348       Eff_corr_Ximinus[i] = geant3fluka_corr_Ximinus[i];
 349     }
 350   }
 351   if(XiStarCase){
 352     double geant3fluka_corr_Xiplus[8]={0.9118, 0.9305, 0.9427, 0.9514, 0.9601, 0.9686, 0.9746, 0.9792};// standard
 353     double geant3fluka_corr_Ximinus[8]={0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906, 0.9906};// standard
 354     //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};
 355     //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};
 356     for(int i=0; i<max_ptbins; i++) {
 357       Eff_corr_Xiplus[i] = geant3fluka_corr_Xiplus[i];
 358       Eff_corr_Ximinus[i] = geant3fluka_corr_Ximinus[i];
 359     }
 360   }
 361
 362   double Eff[max_ptbins]={0};
 363   double Eff_e[max_ptbins]={0};
 364
 365
 366   double spectrum[max_ptbins]={0};
 367   double spectrum_e[max_ptbins]={0};
 368   double chi2perNDF[max_ptbins]={0};
 369   double yield[max_ptbins]={0};
 370   double yield_e[max_ptbins]={0};
 371   double width[max_ptbins]={0};
 372   double width_e[max_ptbins]={0};
 373   double mass[max_ptbins]={0};
 374   double mass_e[max_ptbins]={0};
 375   double MCwidth[max_ptbins]={0};
 376   double MCwidth_e[max_ptbins]={0};
 377   double MCmass[max_ptbins]={0};
 378   double MCmass_e[max_ptbins]={0};
 379   double MissedYield[max_ptbins]={0};
 380   double MissedYield_e[max_ptbins]={0};
 381
 382   TH1F *MCrec_1d[max_ptbins];
 383   TH1F *MCrec_1d_bar[max_ptbins];
 384   TF1 *MCfit = new TF1("MCfit",Voigt,FitRange[0],FitRange[1],4);
 385   double LostRatio=0, LostRatio_e=0;
 386
 387   for(int ptbin=0; ptbin<max_ptbins; ptbin++){
 388     can->cd(ptbin+1);
 389     cout<<"+++++++++++++++ ptbin = "<<ptbin<<"  +++++++++++++++++"<<endl;
 390     int rapbinFirst = Xi3->GetYaxis()->FindBin(raprange[0]+0.01);// bin 16
 391     int rapbinLast = Xi3->GetYaxis()->FindBin(raprange[1]-0.01);// bin 25
 392     int ptbinFirst = Xi3->GetXaxis()->FindBin(ptedges[ptbin+1]+0.01);
 393     int ptbinLast = Xi3->GetXaxis()->FindBin(ptedges[ptbin+1+1]-0.01);
 394
 395     ////////////////////////////
 396     // Efficiency calcualtion
 397     double rec_num=0, input_num=0;
 398     double rec_num1=0, input_num1=0;
 399     double rec_num2=0, input_num2=0;
 400     if(XiStarCase){
 401       int massbinFirst = 1;
 402       int massbinLast = MCrec_XiStar->GetNbinsZ();
 403       rec_num1 = MCrec_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 404       input_num1 = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 405       rec_num2 = MCrec_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 406       input_num2 = MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 407
 408       // set ranges for projections later on
 409       MCrec_XiStar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 410       MCrec_XiStar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 411       MCrec_XiStar->GetZaxis()->SetRange(massbinFirst,massbinLast);
 412       MCrec_XiStar_bar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 413       MCrec_XiStar_bar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 414       MCrec_XiStar_bar->GetZaxis()->SetRange(massbinFirst,massbinLast);
 415       // Estimate Lost number of Xi* from PV cut
 416       double num=0, den=0;
 417       if(ParticleCase==0){
 418         num = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast) + MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 419         den = MCinput_XiStar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast) + MCinput_XiStar_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 420       }else if(ParticleCase==1){
 421         num = MCinput_XiStar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 422         den = MCinput_XiStar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 423       }else {
 424         num = MCinput_XiStar_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 425         den = MCinput_XiStar_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, massbinFirst,massbinLast);
 426       }
 427       LostRatio = num/den;
 428       LostRatio_e = sqrt(pow(sqrt(num)/den,2)+pow(sqrt(den)*num/pow(den,2),2));
 429
 430       ///////////////////////////////////////
 431       // MC fitting for mass and sigma values
 432       TString *nameMCrec=new TString("MCrec_1d_");
 433       TString *nameMCrecbar=new TString("MCrecbar_1d_");
 434       *nameMCrec += ptbin;
 435       *nameMCrecbar += ptbin;
 436       MCrec_1d[ptbin] = (TH1F*)MCrec_XiStar->ProjectionZ(nameMCrec->Data(),ptbinFirst,ptbinLast, rapbinFirst,rapbinLast);
 437       MCrec_1d_bar[ptbin] = (TH1F*)MCrec_XiStar->ProjectionZ(nameMCrecbar->Data(),ptbinFirst,ptbinLast, rapbinFirst,rapbinLast);
 438       MCrec_1d[ptbin]->Add(MCrec_1d_bar[ptbin]);
 439       MCrec_1d[ptbin]->GetXaxis()->SetRangeUser(1.46, 1.62);
 440       MCrec_1d[ptbin]->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
 441       MCrec_1d[ptbin]->SetTitle("p_{T} bin 8");
 442       MCfit->SetParameter(0,.5);// mean
 443       MCfit->SetParameter(1,1.532);// mean
 444       MCfit->SetParameter(2,.002);// Gaussian sigma
 445       MCfit->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
 446       MCfit->SetParLimits(1,1.52,1.54);
 447       MCfit->SetParLimits(2,0.00001,.01);
 448       MCfit->FixParameter(3,.0091);// BW width
 449       MCrec_1d[ptbin]->Fit(MCfit,"IMQ","",FitRange[0],FitRange[1]);
 450       MCwidth[ptbin] = MCfit->GetParameter(2);
 451       MCwidth_e[ptbin] = MCfit->GetParError(2);
 452       MCmass[ptbin] = MCfit->GetParameter(1);
 453       MCmass_e[ptbin] = MCfit->GetParError(1);
 454       //continue;// use this only to plot MC yields one by one
 455       ///////////////////////////////////////////
 456
 457     }else {// Xi
 458       rec_num1 = MCrec_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
 459       input_num1 = MCinput_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCinput_Xi->GetZaxis()->GetNbins());
 460       rec_num2 = MCrec_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi_bar->GetZaxis()->GetNbins());
 461       input_num2 = MCinput_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCinput_Xi_bar->GetZaxis()->GetNbins());
 462       //
 463       MCrec_Xi->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 464       MCrec_Xi->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 465       MCrec_Xi_bar->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 466       MCrec_Xi_bar->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 467       MCrec_1d[ptbin]=(TH1F*)MCrec_Xi->Project3D("z");
 468       MCrec_1d_bar[ptbin]=(TH1F*)MCrec_Xi_bar->Project3D("z");
 469
 470       rec_num1 = MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[1]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[2]-.0001));
 471       rec_num1 -= MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[0]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[1]-.0001));
 472       rec_num1 -= MCrec_1d[ptbin]->Integral(MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[2]),MCrec_1d[ptbin]->GetXaxis()->FindBin(XiCountBound[3]-.0001));
 473
 474       rec_num2 = MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[1]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[2]-.0001));
 475       rec_num2 -= MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[0]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[1]-.0001));
 476       rec_num2 -= MCrec_1d_bar[ptbin]->Integral(MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[2]),MCrec_1d_bar[ptbin]->GetXaxis()->FindBin(XiCountBound[3]-.0001));
 477
 478       // Estimate Lost number of Xi from PV cut
 479       if(ParticleCase==0){
 480         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());
 481         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());
 482       }else if(ParticleCase==1){
 483         LostRatio = MCinput_Xi->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
 484         LostRatio /= MCinput_Xi_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
 485       }else {
 486         LostRatio = MCinput_Xi_bar->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
 487         LostRatio /= MCinput_Xi_bar_prePV->Integral(ptbinFirst,ptbinLast, rapbinFirst,rapbinLast, 1,MCrec_Xi->GetZaxis()->GetNbins());
 488       }
 489
 490
 491     }// else Xi case
 492     cout<<"LostRatio = "<<LostRatio<<" +- "<<LostRatio_e<<endl;
 493
 494
 495     // Apply Geant3/fluka correction
 496     rec_num1 /= Eff_corr_Ximinus[ptbin];
 497     rec_num2 /= Eff_corr_Xiplus[ptbin];
 498
 499
 500     if(ParticleCase == 0){
 501       rec_num = rec_num1 + rec_num2;
 502       input_num = input_num1 + input_num2;
 503     }
 504     if(ParticleCase == 1){
 505       rec_num = rec_num1;
 506       input_num = input_num1;
 507     }
 508     if(ParticleCase == 2){
 509       rec_num = rec_num2;
 510       input_num = input_num2;
 511     }
 512
 513
 514     Eff[ptbin] = rec_num/input_num;
 515     Eff_e[ptbin] = sqrt( pow(sqrt(rec_num)/input_num,2) + pow(sqrt(input_num)*rec_num/pow(input_num,2),2));
 516
 517     cout<<"Efficiency = "<<Eff[ptbin]<<"  +- "<<Eff_e[ptbin]<<"   MC_rec_count = "<<rec_num<<endl;
 518     //cout<<"Reassigning Eff in last pt bin to second to last!!!!!!!!!!!!!"<<endl; Eff[7] = Eff[6];
 519
 520
 521
 522     //cout<<"Re-assigning to David's Eff"<<endl; Eff[ptbin] = David_XiMinus_eff[ptbin];
 523
 524
 525     reject=kTRUE;
 526     TF1 *myBkg=new TF1("myBkg",PolFunction,FitRange[0],FitRange[1],POLdegree+1);
 527     //TF1 *myBkg=new TF1("myBkg",PolFunctionSpecialXi,1.27,1.36,POLdegree+1);// special bkg for Xi
 528     myBkg->SetParameters(6,0);
 529     myBkg->SetLineColor(3);
 530
 531
 532     Xi3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);// pt bins
 533     Xi_bar3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 534     Xi3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 535     Xi_bar3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 536
 537     XiMinus_piPlus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 538     XiPlus_piMinus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 539     XiMinus_piMinus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 540     XiPlus_piPlus3->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 541
 542     XiMinus_piPlus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 543     XiPlus_piMinus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 544     XiMinus_piMinus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 545     XiPlus_piPlus3->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 546
 547
 548     XiMinus_piPlus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 549     XiPlus_piMinus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 550     XiMinus_piMinus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 551     XiPlus_piPlus3_bkg->GetXaxis()->SetRange(ptbinFirst,ptbinLast);
 552
 553     XiMinus_piPlus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 554     XiPlus_piMinus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 555     XiMinus_piMinus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 556     XiPlus_piPlus3_bkg->GetYaxis()->SetRange(rapbinFirst,rapbinLast);
 557
 558
 559     TH1F *Xi=(TH1F*)Xi3->Project3D("z");
 560     TH1F *Xibar=(TH1F*)Xi_bar3->Project3D("z");
 561     if(ParticleCase==0) Xi->Add(Xibar);// to add particle and anti-particle
 562
 563     TH1F *XiMinus_piPlus = (TH1F*)XiMinus_piPlus3->Project3D("z");
 564     TH1F *XiPlus_piMinus = (TH1F*)XiPlus_piMinus3->Project3D("z");
 565     TH1F *XiMinus_piMinus = (TH1F*)XiMinus_piMinus3->Project3D("z");
 566     TH1F *XiPlus_piPlus = (TH1F*)XiPlus_piPlus3->Project3D("z");
 567     if(ParticleCase==0) {
 568       XiMinus_piPlus->Add(XiPlus_piMinus);
 569       XiMinus_piMinus->Add(XiPlus_piPlus);
 570     }
 571
 572
 573     TH1F *XiMinus_piPlus_bkg = (TH1F*)XiMinus_piPlus3_bkg->Project3D("z");
 574     TH1F *XiPlus_piMinus_bkg = (TH1F*)XiPlus_piMinus3_bkg->Project3D("z");
 575     TH1F *XiMinus_piMinus_bkg = (TH1F*)XiMinus_piMinus3_bkg->Project3D("z");
 576     TH1F *XiPlus_piPlus_bkg = (TH1F*)XiPlus_piPlus3_bkg->Project3D("z");
 577     if(ParticleCase==0) {
 578       XiMinus_piPlus_bkg->Add(XiPlus_piMinus_bkg);
 579       XiMinus_piMinus_bkg->Add(XiPlus_piPlus_bkg);
 580     }
 581
 582     TH1F *os=0x0;
 583     TH1F *os_bkg=0x0;
 584     if(XiStarCase) {
 585       if(ParticleCase==0 || ParticleCase==1) {
 586         os=(TH1F*)XiMinus_piPlus->Clone();
 587         os_bkg=(TH1F*)XiMinus_piPlus_bkg->Clone();
 588       }
 589       if(ParticleCase==2) os=(TH1F*)XiPlus_piMinus->Clone();
 590     }else {
 591       if(ParticleCase==0 || ParticleCase==1) os=(TH1F*)Xi->Clone();
 592       if(ParticleCase==2) os=(TH1F*)Xibar->Clone();
 593     }
 594
 595
 596     os->Sumw2();
 597     TH1F *ss=(TH1F*)XiMinus_piMinus->Clone();
 598     TH1F *ss_bkg=(TH1F*)XiMinus_piMinus_bkg->Clone();
 599     ss_bkg->SetLineColor(2);
 600
 601     ss->SetLineColor(4);
 602     ss->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
 603     ss->GetYaxis()->SetTitle("Counts");
 604     //ss->SetTitle(name[ptbin]);
 605     ss->SetTitle(names[ptbin]->Data());
 606
 607     os->SetLineColor(4);
 608     os->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
 609     os->GetYaxis()->SetTitle("Counts");
 610     //os->SetTitle(name[ptbin]);
 611     os->SetTitle(names[ptbin]->Data());
 612
 613
 614     double scalefactor=1;
 615     if(XiStarCase) {
 616       os_bkg->SetLineColor(2);
 617       os_bkg->Sumw2();
 618       double leftS=0, rightS=0;
 619       if(LeftNorm) {leftS=1.49; rightS=1.51;}
 620       else {leftS=1.56; rightS=1.58;}
 621
 622       scalefactor = os->Integral(os->GetXaxis()->FindBin(leftS),os->GetXaxis()->FindBin(rightS));
 623       scalefactor /= os_bkg->Integral(os_bkg->GetXaxis()->FindBin(leftS),os_bkg->GetXaxis()->FindBin(rightS));
 624       os_bkg->Scale(scalefactor);
 625       if(bkgSubtract) os->Add(os_bkg,-1);
 626     }
 627
 628
 629     /////////////////////////////////////
 630     // 1st Iteration
 631     reject=kTRUE;
 632     if(bkgSubtract) {
 633       for(int polbin=0; polbin<POLdegree+1-2; polbin++) myBkg->FixParameter(polbin+2,0);
 634       myBkg->SetParameter(0,0);
 635       myBkg->SetParameter(1,0);
 636     }
 637     /////////////
 638     os->Fit(myBkg,"NQ","",FitRange[0],FitRange[1]);// 1st fit
 639     /////////////
 640
 641
 642     //ss->Draw();
 643     //ss_bkg->Scale(ss->Integral()/ss_bkg->Integral());
 644     //ss_bkg->Draw("same");
 645
 646     TF1 *FullFit_t2;
 647     TF1 *FullFit_t6;
 648
 649
 650
 651     // Voigt
 652     FullFit_t2=new TF1("FullFit_t2",VoigtplusPol,FitRange[0],FitRange[1],PARBINS+1);
 653     FullFit_t2->SetParameter(0,1.0);// mean
 654     FullFit_t2->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
 655     FullFit_t2->SetParameter(1,1.532);// mean
 656     FullFit_t2->SetParLimits(1,1.5,1.56);
 657     //FullFit_t2->FixParameter(2,0);// Gaussian sigma
 658     FullFit_t2->SetParameter(2,.003);// Gaussian sigma
 659     FullFit_t2->SetParLimits(2,0.0001,.015);
 660     FullFit_t2->FixParameter(3,.0091);// BW width
 661     //FullFit_t2->SetParameter(3,.0091);// BW width
 662     //FullFit_t2->SetParLimits(3,.001,.015);
 663     FullFit_t2->SetParName(0,"Norm");
 664     FullFit_t2->SetParName(1,"Mass");
 665     FullFit_t2->SetParName(2,"#sigma");
 666     FullFit_t2->SetParName(3,"#Gamma");
 667
 668     FullFit_t2->SetParameter(4,0);
 669     FullFit_t2->SetParameter(5,0);
 670     FullFit_t2->SetParameter(6,0);
 671     FullFit_t2->SetParameter(7,0);
 672     FullFit_t2->SetParameter(8,0);
 673     FullFit_t2->SetParameter(9,0);
 674
 675
 676
 677
 678     if(XiStarCase) os->GetXaxis()->SetRange(1000*(FitRangeXiStar[0]-1.4),1000*(FitRangeXiStar[1]-1.4));
 679     else os->GetXaxis()->SetRange(1000*(FitRangeXi[0]-1.2),1000*(FitRangeXi[1]+-1.2));
 680
 681     /////////////////////////////////////////////
 682     // 2nd interation
 683     reject=kFALSE;
 684     double pars[PARBINS];
 685     double pars_e[PARBINS];
 686
 687     for(int polbin=0; polbin<PARBINS-offset; polbin++) bkg_params[polbin] = myBkg->GetParameter(polbin);
 688
 689
 690     for(int polbin=0; polbin<PARBINS-offset; polbin++) FullFit_t2->FixParameter(offset+polbin, myBkg->GetParameter(polbin));
 691     //////////////////
 692     os->Fit(FullFit_t2,"IMQ","",FitRange[0],FitRange[1]);// 2nd fit
 693     //////////////////
 694     for(int polbin=0; polbin<PARBINS; polbin++) pars[polbin] = FullFit_t2->GetParameter(polbin);
 695     ////////////////////////////////////////
 696
 697
 698
 699     ///////////////////////////////////////
 700     // 3rd iteration
 701     FullFit_t6=new TF1("FullFit_t6",VoigtplusPol,FitRange[0],FitRange[1],PARBINS);
 702     FullFit_t6->SetParLimits(0,0.01,10.);// 0.0001 to 0.01 for BW only; 0.01 to 10 for true Voigtian
 703     FullFit_t6->SetParLimits(1,1.5,1.56);
 704     FullFit_t6->SetParLimits(2,0.0001,.015);
 705     FullFit_t6->FixParameter(3,.0091);// BW width
 706     FullFit_t6->SetParName(0,"Norm");
 707     FullFit_t6->SetParName(1,"Mass");
 708     FullFit_t6->SetParName(2,"#sigma");
 709     FullFit_t6->SetParName(3,"#Gamma");
 710
 711
 712     for(int polbin=0; polbin<PARBINS; polbin++) FullFit_t6->SetParameter(polbin, pars[polbin]);
 713     if(bkgSubtract) {
 714       for(int polbin=0; polbin<PARBINS-offset-2; polbin++) FullFit_t6->FixParameter(polbin+offset+2,0);
 715       FullFit_t6->SetParameter(offset,0);
 716       FullFit_t6->SetParameter(offset+1,0);
 717     }
 718     //////////////////
 719     os->Fit(FullFit_t6,"IMEQ","",FitRange[0],FitRange[1]);// 3rd fit
 720     //////////////////
 721     for(int polbin=0; polbin<PARBINS-offset; polbin++) bkg_params[polbin] = FullFit_t6->GetParameter(polbin+offset);
 722     for(int polbin=0; polbin<PARBINS; polbin++) pars[polbin] = FullFit_t6->GetParameter(polbin);
 723     for(int polbin=0; polbin<PARBINS; polbin++) pars_e[polbin] = FullFit_t6->GetParError(polbin);
 724
 725
 726
 727     /////////////////////////////////////////
 728     if(!bkgSubtract) os_bkg->Draw("same");// XiStar
 729     for(int polbin=0; polbin<POLdegree+1; polbin++) myBkg->SetParameter(polbin, bkg_params[polbin]);
 730
 731     for(int polbin=0; polbin<POLdegree+1; polbin++) myPol[ptbin]->FixParameter(polbin, bkg_params[polbin]);
 732     myPol[ptbin]->DrawCopy("same");
 733     if((ptbin+1)==max_ptbins){
 734       legend->AddEntry(os,"opp-charge","l");
 735       legend->AddEntry(os_bkg,"opp-charge event mixing","l");
 736       //legend->Draw("same");
 737     }
 738
 739
 740     double temp_yield=0, temp_yield_e=0;
 741     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.));
 742     MissedYield_e[ptbin] = MissedYield[ptbin]*FullFit_t6->GetParError(0)/FullFit_t6->GetParameter(0);
 743     //cout<<"Signal/Bkg = "<<temp_yield/(1000.*myBkg->Integral(XiStarYieldCountRange[0],XiStarYieldCountRange[1]))<<endl;
 744     //cout<<"Significance = "<<temp_yield/sqrt(temp_yield)<<endl;
 745
 746     if(VoigtFitOnly) {
 747       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
 748       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)
 749       temp_yield_e = temp_yield*FullFit_t6->GetParError(0)/FullFit_t6->GetParameter(0);
 750     }else {// yield from bin counting
 751       if(XiStarCase){
 752         temp_yield = os->Integral(os->GetXaxis()->FindBin(XiStarCountBound[2]),os->GetXaxis()->FindBin(XiStarCountBound[3]-.0001));
 753         temp_yield -= myBkg->Integral(XiStarCountBound[2],XiStarCountBound[3]);
 754         cout<<"Missed yield fraction = "<<MissedYield[ptbin]/(temp_yield+MissedYield[ptbin])<<endl;
 755         temp_yield += MissedYield[ptbin];
 756         for(int massBin=os->GetXaxis()->FindBin(XiStarCountBound[2]); massBin<=os->GetXaxis()->FindBin(XiStarCountBound[3]); massBin++){
 757           temp_yield_e += pow(os->GetBinError(massBin),2);
 758         }
 759         temp_yield_e += pow(MissedYield_e[ptbin],2);
 760         temp_yield_e = sqrt(temp_yield_e);
 761         cout<<"Included/Total = "<<(temp_yield-MissedYield[ptbin])/temp_yield<<endl;
 762       }else {// Xi case
 763         temp_yield = os->Integral(os->GetXaxis()->FindBin(XiCountBound[1]),os->GetXaxis()->FindBin(XiCountBound[2]-.0001));
 764         temp_yield -= os->Integral(os->GetXaxis()->FindBin(XiCountBound[0]),os->GetXaxis()->FindBin(XiCountBound[1]-.0001));
 765         temp_yield -= os->Integral(os->GetXaxis()->FindBin(XiCountBound[2]),os->GetXaxis()->FindBin(XiCountBound[3]-.0001));
 766       }
 767     }
 768     cout<<"yield = "<<temp_yield<<" +- "<<temp_yield_e<<endl;
 769
 770
 771
 772     double base = Eff[ptbin]*LostRatio*(raprange[1]-raprange[0])*(2*pt_points_e[ptbin])*Events_PhysSel->GetEntries()/0.852;
 773
 774     spectrum[ptbin] = temp_yield;
 775     spectrum[ptbin] /= base;
 776     //
 777     spectrum_e[ptbin] = pow(temp_yield_e/base,2);
 778     spectrum_e[ptbin] += pow(temp_yield/base*Eff_e[ptbin]/Eff[ptbin],2);
 779     spectrum_e[ptbin] += pow(temp_yield/base*LostRatio_e/LostRatio,2);
 780     spectrum_e[ptbin] += pow(temp_yield/base*sqrt(Events_PhysSel->GetEntries())/Events_PhysSel->GetEntries(),2);
 781     spectrum_e[ptbin] = sqrt(spectrum_e[ptbin]);
 782     //spectrum_e[ptbin] = temp_yield/base*Eff_e[ptbin]/Eff[ptbin];// for eff errors only
 783     if(ParticleCase==0) spectrum[ptbin] /= 2.;// (particle + cc)/2
 784     if(ParticleCase==0) spectrum_e[ptbin] /= 2.;// (particle + cc)/2
 785     //
 786
 787
 788     cout<<"spectrum point = "<<spectrum[ptbin]<<"  +- "<<spectrum_e[ptbin]<<endl;
 789
 790     yield[ptbin] = temp_yield;
 791     yield_e[ptbin] = temp_yield_e;
 792     mass[ptbin] = pars[1];
 793     mass_e[ptbin] = pars_e[1];
 794     width[ptbin] = pars[2];
 795     width_e[ptbin] = pars_e[2];
 796
 797     chi2perNDF[ptbin] = FullFit_t6->GetChisquare()/double(FullFit_t6->GetNDF());
 798
 799   }
 800
 801
 802
 803   TH1D *RawYields = new TH1D("RawYields","",8,0.5,8.5);
 804   for(int ii=0; ii<8; ii++){ RawYields->SetBinContent(ii+1, yield[ii]); RawYields->SetBinError(ii+1, yield_e[ii]);}
 805   TH1D *Efficiency = new TH1D("Efficiency","",8,0.5,8.5);
 806   for(int ii=0; ii<8; ii++){ Efficiency->SetBinContent(ii+1, Eff[ii]); Efficiency->SetBinError(ii+1, Eff_e[ii]);}
 807
 808   /*
 809   TLatex *tex2 = new TLatex(1.6,1.5,"#Xi(1530)^{0} #rightarrow #Xi^{+-} + #pi^{-+}");
 810   tex2->SetTextSize(.06);
 811   tex2->Draw();
 812
 813   TLatex *tex3 = new TLatex(1.6,1.5,"ALICE Performance");
 814   tex3->SetTextSize(.06);
 815   tex3->SetTextColor(2);
 816   tex3->Draw();
 817
 818   TLatex *tex4 = new TLatex(1.6,1.5,"12/09/2011");
 819   tex4->SetTextSize(.04);
 820   tex4->Draw();
 821
 822   TLatex *tex5 = new TLatex(1.6,1.5,"pp, #sqrt{s} = 7 TeV");
 823   tex5->SetTextSize(.06);
 824   tex5->Draw();
 825
 826   TPad *c1_1 = new TPad("c1_1", "c1_1",0.25,0.3,0.4,0.5);
 827   c1_1->SetBorderMode(0);
 828   c1_1->SetLineColor(0);
 829   c1_1->SetFillColor(0);
 830   c1_1->Draw();
 831
 832   c1_1->cd();
 833   TImage *alicelogo = TImage::Open("../../Pictures/alice_logo_transparent2.png");
 834   //alicelogo->SetConstRatio(kTRUE);
 835   //alicelogo->SetImageQuality(TAttImage::kImgBest);
 836   alicelogo->Draw("same");
 837   */
 838
 839
 840   TCanvas *can2 = new TCanvas("can2","can2",13,34,700,500);
 841   gStyle->SetOptFit(111);
 842   can2->Range(-1.25,-0.2625,11.25,2.3625);
 843   can2->SetFillColor(10);
 844   can2->SetBorderMode(0);
 845   can2->SetBorderSize(2);
 846   can2->SetGridx();
 847   can2->SetGridy();
 848   //can2->SetLogy();
 849   can2->SetFrameFillColor(0);
 850   can2->SetFrameBorderMode(0);
 851   can2->SetFrameBorderMode(0);
 852
 853
 854
 855   TH1D *h_Xispectrum = new TH1D("h_Xispectrum","Xi spectrum",max_ptbins+1,ptedges);
 856   h_Xispectrum->SetMarkerStyle(20);
 857   h_Xispectrum->SetMarkerSize(1.0);
 858   if(XiStarCase) h_Xispectrum->SetTitle("#Xi(1530) spectrum");
 859   else h_Xispectrum->SetTitle("#Xi spectrum");
 860   h_Xispectrum->SetTitle("pp #sqrt{s}= 7 TeV");
 861   h_Xispectrum->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 862   h_Xispectrum->GetYaxis()->SetTitle("1/N_{E}d^{2}N/dydp_{T} (GeV/c)^{-1}");
 863   //h_Xispectrum->GetYaxis()->SetTitle("(Real N_Xi-)/(MC N_Xi-)");
 864   h_Xispectrum->SetMinimum(3.0e-6);
 865   h_Xispectrum->SetMaximum(0.02);
 866   h_Xispectrum->GetXaxis()->SetLimits(0,5.6);
 867
 868   h_Xispectrum->SetBinContent(1,+100);// first bin is just there to visualize pt=0
 869   for(int i=0; i<max_ptbins; i++) {
 870     h_Xispectrum->SetBinContent(i+1+1, spectrum[i]);
 871     h_Xispectrum->SetBinError(i+1+1, spectrum_e[i]);
 872   }
 873
 874   //h_Xispectrum->Draw();
 875
 876
 877   double minFitpoint=0, maxFitpoint=0;
 878   if(XiStarCase) {minFitpoint = 0.8; maxFitpoint = 5.6;}// 0.8 to 5.6
 879   else {minFitpoint = 0.6; maxFitpoint = 8.5;}
 880
 881
 882   TF1 *myLevy=new TF1("myLevy",myLevyPt,0,8.5,3);
 883   myLevy->SetParName(0,"dN/dy");
 884   myLevy->SetParName(1,"C");
 885   myLevy->SetParName(2,"n");
 886   myLevy->SetParameter(0,.008);
 887   myLevy->SetParameter(1,.3);
 888   myLevy->SetParameter(2,15);
 889   myLevy->SetParLimits(0,.001,.2);
 890   myLevy->SetParLimits(1,.1,1);
 891   myLevy->SetParLimits(2,1,500);
 892
 893
 894   h_Xispectrum->Draw();
 895   h_Xispectrum->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
 896   //h_Xispectrum->Fit(myExp,"IME","",minFitpoint,maxFitpoint);
 897   //myLevy->Draw("same");
 898   //TH1D *temp = (TH1D*)myLevy->GetHistogram();
 899
 900   //legend2->AddEntry(h_Xispectrum,"#Xi^{+}, This analysis","p");
 901   //legend2->AddEntry(h_Xispectrum,"10b+10c+10d","p");
 902   legend2->AddEntry(h_Xispectrum,"(#Xi(1530) + cc)/2","p");
 903
 904   cout<<"dN/dy = "<<myLevy->Integral(0,10)<<"  , Covered range = "<<myLevy->Integral(0.8,5.6)<<endl;
 905
 906
 907   double Data_Levy_Ratio[8]={0};
 908   double Data_Levy_Ratio_e[8]={0};
 909   for(int i=0; i<8; i++){
 910     Data_Levy_Ratio[i] = (ptedges[i+2]-ptedges[i+1])*spectrum[i]/myLevy->Integral(ptedges[i+1],ptedges[i+2]);
 911     Data_Levy_Ratio_e[i] = (ptedges[i+2]-ptedges[i+1])*spectrum_e[i]/myLevy->Integral(ptedges[i+1],ptedges[i+2]);
 912   }
 913   TGraphErrors *gr_Data_Levy_Ratio = new TGraphErrors(8, pt_points, Data_Levy_Ratio, pt_points_e, Data_Levy_Ratio_e);
 914   gr_Data_Levy_Ratio->SetMarkerStyle(20);
 915   gr_Data_Levy_Ratio->SetTitle("");
 916   gr_Data_Levy_Ratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 917   gr_Data_Levy_Ratio->GetYaxis()->SetTitle("(Measured Yield)/(Levy fit integral)");
 918   //gr_Data_Levy_Ratio->Draw("AP");
 919
 920   // Reference Spectra
 921   double LHC10b_XiStar[8]={0.00175909, 0.000893716, 0.000957646, 0.000399652, 0.000217706, 6.86327e-05, 3.24599e-05, 1.28201e-05};
 922   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};
 923   double LHC10c_XiStar[8]={0.00153422, 0.00106629, 0.000752749, 0.000495777, 0.000255897, 6.88883e-05, 2.41039e-05, 1.09414e-05};
 924   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};
 925   double LHC10d_XiStar[8]={0.00137741, 0.00118096, 0.000813771, 0.000465871, 0.000223755, 6.5902e-05, 2.2662e-05, 7.11972e-06};
 926   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};
 927   TGraphErrors *gr_10bXiStar = new TGraphErrors(8, pt_points, LHC10b_XiStar, pt_points_e, LHC10b_XiStar_e);
 928   gr_10bXiStar->SetMarkerStyle(21);
 929   gr_10bXiStar->SetMarkerColor(2);
 930   gr_10bXiStar->SetLineColor(2);
 931   TGraphErrors *gr_10cXiStar = new TGraphErrors(8, pt_points, LHC10c_XiStar, pt_points_e, LHC10c_XiStar_e);
 932   gr_10cXiStar->SetMarkerStyle(22);
 933   gr_10cXiStar->SetMarkerColor(4);
 934   gr_10cXiStar->SetLineColor(4);
 935   TGraphErrors *gr_10dXiStar = new TGraphErrors(8, pt_points, LHC10d_XiStar, pt_points_e, LHC10d_XiStar_e);
 936   gr_10dXiStar->SetMarkerStyle(23);
 937   gr_10dXiStar->SetMarkerColor(6);
 938   gr_10dXiStar->SetLineColor(6);
 939   //gr_10bXiStar->Draw("P");
 940   //gr_10cXiStar->Draw("P");
 941   //gr_10dXiStar->Draw("P");
 942   //legend2->AddEntry(gr_10bXiStar,"LHC10b","p");
 943   //legend2->AddEntry(gr_10cXiStar,"LHC10c","p");
 944   //legend2->AddEntry(gr_10dXiStar,"LHC10d","p");
 945
 946   // Old QM2011 values for (Xi* + cc)/2
 947   double Old_XiStarPoints[8]={0.001259, 0.0008586, 0.0006316, 0.0003604, 1.856e-4, 5.4601e-5, 2.0187e-5, 7.6559e-6};
 948   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};
 949   TGraphErrors *gr_OldXiStar = new TGraphErrors(8, pt_points, Old_XiStarPoints, pt_points_e, Old_XiStarPoints_e);
 950   gr_OldXiStar->SetMarkerStyle(21);
 951   gr_OldXiStar->SetMarkerColor(2);
 952   gr_OldXiStar->SetLineColor(2);
 953   //gr_OldXiStar->Draw("P");
 954   //legend2->AddEntry(gr_OldXiStar,"(#Xi(1530) + cc)/2: QM11 Results","p");
 955   /*
 956   // To fit the Old QM results
 957   TH1D *h_OldXiStar = (TH1D*)h_Xispectrum->Clone();
 958   h_OldXiStar->Add(h_OldXiStar,-1);
 959   for(int i=0; i<8; i++){
 960     h_OldXiStar->SetBinContent(i+2, Old_XiStarPoints[i]);
 961     h_OldXiStar->SetBinError(i+2, Old_XiStarPoints_e[i]);
 962   }
 963   h_OldXiStar->SetBinContent(1,+100);// first bin has no data
 964   h_OldXiStar->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
 965   cout<<"old mean pt = "<<myLevy->Moment(1,0.,10)<<endl;
 966   */
 967   //double pt_points_shifted[max_ptbins]={0};
 968   //TF1 *AvgX = new TF1("AvgX","x*myExp",0.6,5.6);
 969
 970   //AvgX->SetParName(0,"A");
 971   //AvgX->SetParName(1,"T");
 972   //for(int i=0; i<8; i++){
 973   //cout<<AvgX->Integral(ptedges[i],ptedges[i+1])/myExp->Integral(ptedges[i],ptedges[i+1])<<endl;
 974   //}
 975
 976
 977
 978
 979
 980   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};
 981
 982
 983   // David's new points
 984   // David's Xi-
 985   //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};
 986   //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};
 987   // David's Xi+
 988   //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};
 989   //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};
 990   // Published ALICE points
 991   // Xi-
 992   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};
 993   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};
 994   // Xi+
 995   //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};
 996   //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};
 997
 998   TGraphErrors *gr_DavidXi = new TGraphErrors(18, David_points_pt, David_points_y, 0, David_points_y_e);
 999   gr_DavidXi->SetMarkerStyle(24);
1000   gr_DavidXi->SetMarkerSize(2);
1001   gr_DavidXi->SetMarkerColor(4);
1002   //gr_DavidXi->Draw("P");
1003   //legend2->AddEntry(gr_DavidXi,"ALICE Published #Xi^{+}","p");
1004   //legend2->AddEntry(gr_DavidXi,"#Xi^{+-}","p");
1005
1006
1007
1008   //cout<<"Efficiencies:"<<endl;
1009   //for(int i=0; i<max_ptbins; i++) cout<<Eff[i]<<", ";
1010   //cout<<endl;
1011   //cout<<"Corresponding errors:"<<endl;
1012   //for(int i=0; i<max_ptbins; i++) cout<<Eff_e[i]<<", ";
1013   //cout<<endl;
1014   //cout<<"Raw Yields:"<<endl;
1015   //for(int i=0; i<max_ptbins; i++) cout<<yield[i]<<", ";
1016   //cout<<endl;
1017   //cout<<"Ratio of My spectrum to Davids:"<<endl;
1018   //for(int i=0; i<max_ptbins; i++) cout<<spectrum[i]/David_points_y[i]<<", ";
1019   //cout<<endl;
1020   //cout<<"Corresponding ratio errors:"<<endl;
1021   //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))<<", ";
1022   //cout<<endl;
1023   //cout<<"Ratio of My Eff to Davids:"<<endl;
1024   //for(int i=0; i<max_ptbins; i++) cout<<Eff[i]/David_XiMinus_eff[i]<<", ";
1025   //cout<<endl;
1026   //cout<<"Corresponding ratio errors:"<<endl;
1027   //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))<<", ";
1028   //cout<<endl;
1029
1030
1031   //Xi 10b+10c+10d with first eff bug fix, second bug fix is a very tiny change
1032   //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};
1033   // My eff with David's cuts XiMinus
1034   //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};
1035   //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};
1036   //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};
1037   //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};
1038
1039   // 10d1 eff
1040   double Eff_ref1[8]={0.00585138, 0.0271111, 0.0370528, 0.0767083, 0.0970953, 0.129362, 0.144039, 0.130414};
1041   double Eff_ref1_e[8]={0.00110308, 0.00301515, 0.00433263, 0.00845994, 0.00972208, 0.0175748, 0.0318935, 0.0367763};
1042   // 10d4 eff
1043   double Eff_ref2[8]={0.00646719, 0.0236728, 0.0454735, 0.0692855, 0.0922256, 0.13542, 0.161154, 0.126502};
1044   double Eff_ref2_e[8]={0.000740606, 0.00176465, 0.00315074, 0.00498191, 0.00600009, 0.0117801, 0.0203019, 0.0247307};
1045   // 10f6a eff
1046   double Eff_ref3[8]={0.00589799, 0.0228898, 0.043247, 0.0692314, 0.104212, 0.14611, 0.1587, 0.176539};
1047   double Eff_ref3_e[8]={0.000465501, 0.00114491, 0.00199253, 0.00333088, 0.00421532, 0.00822058, 0.0133099, 0.0202604};
1048   // 10f6a eff doubled binning
1049   //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};
1050   //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};
1051
1052   // Pythia to Phojet eff ratio
1053   for(int ii=0; ii<max_ptbins; ii++){
1054     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));
1055     Eff[ii] /= Eff_ref3[ii];
1056   }
1057
1058   TGraphErrors *gr_Eff = new TGraphErrors(max_ptbins,pt_points, Eff, pt_points_e, Eff_e);
1059   gr_Eff->SetMarkerStyle(20);
1060   gr_Eff->SetMinimum(0);
1061   gr_Eff->SetMaximum(.4);
1062   //gr_Eff->GetXaxis()->SetLimits(0,6.1);
1063   gr_Eff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1064   gr_Eff->GetYaxis()->SetTitle("Efficiency");
1065   //gr_Eff->SetTitle("#Xi(1530)^{0}+cc Efficiency");
1066   //gr_Eff->Draw("AP");
1067   //legend2->AddEntry(gr_Eff,"10d1+10d4+10f6a","p");
1068   //legend2->AddEntry(gr_Eff,"10f6 (Phojet)","p");
1069   //gr_Eff->Fit("pol0","IME","",1.6,4.8);
1070
1071   TGraphErrors *gr_Eff_ref1 = new TGraphErrors(max_ptbins,pt_points, Eff_ref1, pt_points_e, Eff_ref1_e);
1072   gr_Eff_ref1->SetMarkerStyle(20);
1073   gr_Eff_ref1->SetMarkerColor(2);
1074   gr_Eff_ref1->SetLineColor(2);
1075   //gr_Eff_ref1->Draw("P");
1076   //legend2->AddEntry(gr_Eff_ref1,"10d1 MC","p");
1077   TGraphErrors *gr_Eff_ref2 = new TGraphErrors(max_ptbins,pt_points, Eff_ref2, pt_points_e, Eff_ref2_e);
1078   gr_Eff_ref2->SetMarkerStyle(20);
1079   gr_Eff_ref2->SetMarkerColor(4);
1080   gr_Eff_ref2->SetLineColor(4);
1081   //gr_Eff_ref2->Draw("P");
1082   //legend2->AddEntry(gr_Eff_ref2,"10d4 MC","p");
1083   TGraphErrors *gr_Eff_ref3 = new TGraphErrors(max_ptbins,pt_points, Eff_ref3, pt_points_e, Eff_ref3_e);
1084   gr_Eff_ref3->SetMarkerStyle(20);
1085   gr_Eff_ref3->SetMarkerColor(3);
1086   gr_Eff_ref3->SetLineColor(3);
1087   //gr_Eff_ref3->Draw("P");
1088   //legend2->AddEntry(gr_Eff_ref3,"10f6a (Pythia Perugia)","p");
1089
1090
1091   //TGraphErrors *gr_Eff_xiplus = new TGraphErrors(max_ptbins,pt_points, Eff_refXiPlus, 0, Eff_refXiPlus_e);
1092   //gr_Eff_xiplus->SetMarkerStyle(21);
1093   //gr_Eff_xiplus->SetMarkerColor(4);
1094   //gr_Eff_xiplus->GetXaxis()->SetTitle("pt (GeV/c)");
1095   //gr_Eff_xiplus->GetYaxis()->SetTitle("Efficiency");
1096   //gr_Eff_xiplus->SetTitle("#Xi(1530)^0 Efficiency");
1097   //gr_Eff_xiplus->Draw("P");
1098
1099
1100
1101   TGraphErrors *gr_David_Eff_XiMinus = new TGraphErrors(18,pt_points_David,David_XiMinus_eff,0,David_XiMinus_eff_e);
1102   gr_David_Eff_XiMinus->SetMarkerStyle(24);
1103   gr_David_Eff_XiMinus->SetMarkerColor(2);
1104   //gr_David_Eff_XiMinus->Draw("P");
1105   //legend2->AddEntry(gr_David_Eff_XiMinus,"David/Antonin Xi- Eff","p");
1106
1107   TGraphErrors *gr_David_Eff_XiPlus = new TGraphErrors(18,pt_points_David,David_XiPlus_eff,0,David_XiPlus_eff_e);
1108   gr_David_Eff_XiPlus->SetMarkerStyle(25);
1109   gr_David_Eff_XiPlus->SetMarkerColor(4);
1110   //gr_David_Eff_XiPlus->Draw("P");
1111
1112
1113   legend2->Draw("same");
1114
1115
1116   TLegend *legend3 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
1117   legend3->SetBorderSize(1);
1118   legend3->SetTextSize(.04);// small .03; large .036
1119   //legend3->SetLineColor(0);
1120   legend3->SetFillColor(0);
1121
1122
1123   TGraphErrors *gr_RealWidth = new TGraphErrors(8, pt_points, width, pt_points_e, width_e);
1124   gr_RealWidth->SetMarkerStyle(20);
1125   gr_RealWidth->SetMinimum(0);
1126   gr_RealWidth->SetMaximum(0.0035);
1127   gr_RealWidth->GetYaxis()->SetTitle("#sigma (GeV/c)");
1128   gr_RealWidth->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1129   gr_RealWidth->SetTitle("");
1130   //gr_RealWidth->Draw("AP");
1131   //legend3->AddEntry(gr_RealWidth,"Real data (10b+10c+10d)","p");
1132
1133
1134   TGraphErrors *gr_MCWidth = new TGraphErrors(8, pt_points, MCwidth, pt_points_e, MCwidth_e);
1135   gr_MCWidth->SetMarkerStyle(20);
1136   gr_MCWidth->SetMarkerColor(2);
1137   gr_MCWidth->SetLineColor(2);
1138   gr_MCWidth->GetYaxis()->SetTitle("#sigma (GeV/c)");
1139   gr_MCWidth->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1140   gr_MCWidth->SetTitle("");
1141   //gr_MCWidth->Draw("P");
1142   //legend3->AddEntry(gr_MCWidth,"Pythia+ALICE (10d1+10d4+10f6a)","p");
1143
1144
1145
1146
1147   TGraphErrors *gr_RealMass = new TGraphErrors(8, pt_points, mass, pt_points_e, mass_e);
1148   gr_RealMass->SetMarkerStyle(20);
1149   gr_RealMass->SetMinimum(1.529);
1150   gr_RealMass->SetMaximum(1.5345);
1151   gr_RealMass->SetTitle("");
1152   gr_RealMass->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1153   gr_RealMass->GetYaxis()->SetTitleOffset(1.3);
1154   gr_RealMass->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
1155   //gr_RealMass->Draw("AP");
1156   //legend3->AddEntry(gr_RealMass,"Real data (10b+10c+10d)","p");
1157
1158   TGraphErrors *gr_MCMass = new TGraphErrors(8, pt_points, MCmass, pt_points_e, MCmass_e);
1159   gr_MCMass->SetMarkerStyle(20);
1160   gr_MCMass->SetMarkerColor(2);
1161   gr_MCMass->SetTitle("");
1162   gr_MCMass->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
1163   gr_MCMass->GetYaxis()->SetTitleOffset(1.3);
1164   gr_MCMass->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1165   //gr_MCMass->Draw("P");
1166   //legend3->AddEntry(gr_MCMass,"Pythia+ALICE (10d1+10d4+10f6a)","p");
1167   //TF1 *pdg_massline=new TF1("pdg_massline","1.53178",0,10);
1168   //pdg_massline->Draw("same");
1169
1170   //legend3->Draw("same");
1171   TF1 *myLevy2=new TF1("myLevy2",myLevyPt,0,8.5,3);
1172   myLevy2->SetLineColor(4);
1173   myLevy2->SetParName(0,"dN/dy");
1174   myLevy2->SetParName(1,"C");
1175   myLevy2->SetParName(2,"n");
1176   myLevy2->SetParameter(0,.003);
1177   myLevy2->SetParameter(1,.3);
1178   myLevy2->SetParameter(2,15);
1179   myLevy2->SetParLimits(0,0.0001,0.01);
1180   myLevy2->SetParLimits(1,.1,1);
1181   myLevy2->SetParLimits(2,1,500);
1182   TF1 *myExp2 = new TF1("myExp2",myExpFit,0,8.5,2);
1183   myExp2->SetLineColor(4);
1184   myExp2->SetParName(0,"dN/dy");
1185   myExp2->SetParName(1,"C");
1186   myExp2->SetParameter(0,.003);
1187   myExp2->SetParameter(1,.3);
1188   myExp2->SetParLimits(0,0.0001,0.1);
1189   myExp2->SetParLimits(1,.1,1);
1190
1191
1192   MCinput_Spectrum->SetMarkerStyle(20);
1193   MCinput_Spectrum->SetMarkerColor(2);
1194   MCinput_Spectrum->SetLineColor(2);
1195   MCinput_Spectrum->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1196   MCinput_Spectrum->GetYaxis()->SetTitle("1/N_{E}d^{2}N/dydp_{T} (GeV/c)^{-1}");
1197   //MCinput_Spectrum->SetTitle("Pythia Perugia");
1198   //MCinput_Spectrum->SetTitle("Phojet");
1199
1200   MCinput_Spectrum->Sumw2();
1201   MCinput_Spectrum->Scale(1/MCinput_Spectrum->GetXaxis()->GetBinWidth(1));
1202   MCinput_Spectrum->Scale(1/EventsMC_postPV->GetEntries());
1203
1204   // Scale Pythia to Data fit Tsallis
1205   MCinput_Spectrum->Scale(myLevy->Integral(0.8,1.2)/(MCinput_Spectrum->Integral(9,12)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)));
1206   //
1207   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;
1208   cout<<"high pt extrapolation over total Levy = "<<myLevy->Integral(5.6,10.)/myLevy->Integral(0,10.)<<endl;
1209   cout<<"high pt Pythia over total Pythia = "<<MCinput_Spectrum->Integral(56,100)/MCinput_Spectrum->Integral(1,100)<<endl;
1210   cout<<"Total yield Pythia over total Levy = "<<MCinput_Spectrum->Integral(1,100)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)/myLevy->Integral(0,10.)<<endl;
1211
1212   //MCinput_Spectrum->Draw();
1213   //minFitpoint=0.; maxFitpoint=10.;
1214   //MCinput_Spectrum->Fit(myLevy2,"IME","",minFitpoint,maxFitpoint);
1215   //MCinput_Spectrum->Fit(myExp,"IME","",minFitpoint,maxFitpoint);
1216   //myExp->Draw("same");
1217   //myLevy2->Draw("same");
1218   //myLevy->Draw("same");
1219   //legend3->AddEntry(MCinput_Spectrum,"Pythia Perugia (Scaled to Real Data)","p");
1220   //legend3->AddEntry(myLevy2,"Phojet Levy fit (fit range from pt=0.8 to 5.6 GeV/c","l");
1221   //legend3->AddEntry(myLevy,"Real data Levy fit","l");
1222
1223   double xaxis[60]={0};
1224   double yaxis[60]={0};
1225   for(int ii=0; ii<60; ii++){
1226     xaxis[ii] = MCinput_Spectrum->GetXaxis()->GetBinCenter(ii+1);
1227     yaxis[ii] = myLevy->Eval(xaxis[ii])/myLevy2->Eval(xaxis[ii]);
1228   }
1229   TGraph *LevyRatio = new TGraph(60, xaxis, yaxis);
1230   LevyRatio->SetTitle("Levy Ratio");
1231   LevyRatio->GetXaxis()->SetTitle("p_{T}");
1232   LevyRatio->GetYaxis()->SetTitle("(Real Data Levy)/(Pythia fit Levy)");
1233   LevyRatio->SetMarkerStyle(20);
1234   //LevyRatio->Draw("AP");
1235
1236
1237   //legend3->Draw("same");
1238
1239   //cout<<"low-pt extrapolation from Levy = "<<myLevy2->Integral(0,0.8)<<"   Histogram Integral = "<<MCinput_Spectrum->Integral(1,8)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
1240   //cout<<"Total Bin Count yield = "<<MCinput_Spectrum->Integral(1,MCinput_Spectrum->GetNbinsX())*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
1241   //cout<<"low-pt extrapolation from Exp = "<<myExp->Integral(0,0.8)<<"   Histogram Integral = "<<MCinput_Spectrum->Integral(1,8)*MCinput_Spectrum->GetXaxis()->GetBinWidth(1)<<endl;
1242
1243
1244   if(SaveToFile){
1245     TFile *outputfile = new TFile(outfilename->Data(),"RECREATE");
1246     can->Write("can");
1247     can2->Write("can2");
1248     myLevy->Write("myLevy");
1249     h_Xispectrum->Write("h_Spectrum");
1250     RawYields->Write("RawYields");
1251     Efficiency->Write("Efficiency");
1252
1253     outputfile->Close();
1254   }
1255
1256   cout<<endl;
1257
1258 }
1259
1260 //________________________________________________________________________
1261 double PolFunction(double *x, double *par){
1262
1263   if (reject && x[0] > yieldRange[0] && x[0] < yieldRange[1]) {
1264     TF1::RejectPoint();
1265     return 0;
1266   }
1267
1268   if(POLdegree==1) return par[0] + par[1]*x[0];
1269   else if(POLdegree==2) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2);
1270   else if(POLdegree==3) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3);
1271   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);
1272   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);
1273   else return 0;
1274
1275 }
1276 //________________________________________________________________________
1277 double PolFunctionSpecialXi(double *x, double *par){
1278
1279   if (x[0] > 1.28 && x[0] < 1.3) {
1280     TF1::RejectPoint();
1281     return 0;
1282   }
1283   if (reject && x[0] > 1.31 && x[0] < 1.334) {
1284     TF1::RejectPoint();
1285     return 0;
1286   }
1287
1288   if(POLdegree==1) return par[0] + par[1]*x[0];
1289   else if(POLdegree==2) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2);
1290   else if(POLdegree==3) return par[0] + par[1]*x[0] + par[2]*pow(x[0],2) + par[3]*pow(x[0],3);
1291   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);
1292   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);
1293   else return 0;
1294
1295 }
1296 //________________________________________________________________________
1297 double BWFunction(double *x, double *par){
1298   return (par[0]*par[2]/(1000.*2.*3.1415926))/( pow(x[0]-par[1],2) + par[2]*par[2]/4.);
1299   //return (par[0]*par[1]/(1000.*2.*3.1415926))/( pow(x[0]-par[2],2) + par[1]*par[1]/4.);
1300 }
1301 //________________________________________________________________________
1302 double BWplusPol(double *x, double *par){
1303   return BWFunction(x,par) + PolFunction(x,&par[3]);
1304 }
1305 //________________________________________________________________________
1306 double GausFunction(double *x, double *par){
1307   return (par[0]*.001)/(sqrt(2*3.14159*par[2]*par[2]))*exp(-pow((x[0]-par[1])/(sqrt(2)*par[2]),2));
1308 }
1309 //________________________________________________________________________
1310 double GausplusPol(double *x, double *par){
1311   return GausFunction(x,par) + PolFunction(x,&par[3]);
1312 }
1313 //________________________________________________________________________
1314 double GausplusPolSpecialXi(double *x, double *par){
1315   return GausFunction(x,par) + PolFunctionSpecialXi(x,&par[3]);
1316 }
1317 //________________________________________________________________________
1318 double myLevyPt(Double_t *x, Double_t *par)
1319 {
1320   Double_t lMass=0;
1321   if(XiStarCase) lMass = 1.5318; //Xi mass
1322   else lMass = 1.32171; //Xi mass
1323
1324   Double_t ldNdy  = par[0]; // dN/dy
1325   Double_t l2pi   = 2*TMath::Pi(); // 2pi
1326   Double_t lTemp = par[1]; // Temperature
1327   Double_t lPower = par[2]; // power=n
1328
1329   Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
1330   Double_t lInPower = 1 + (TMath::Sqrt(x[0]*x[0]+lMass*lMass)-lMass) / (lPower*lTemp);
1331
1332   return l2pi * ldNdy * x[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
1333 }
1334 //________________________________________________________________________
1335 double myExpFit(Double_t *x, Double_t *par)
1336 {
1337   Double_t lMass=0;
1338   if(XiStarCase) lMass = 1.5318; //Xi mass
1339   else lMass = 1.32171; //Xi mass
1340
1341   return 2*TMath::Pi()*x[0]*par[0]*exp(-TMath::Sqrt(x[0]*x[0]+lMass*lMass)/par[1]);
1342 }
1343 //________________________________________________________________________
1344 double Voigt(Double_t *x, Double_t *par)
1345 {// code taken directly from RooVoigtian in ROOT
1346
1347   Double_t Norm = par[0];
1348   Double_t mean = par[1];
1349   Double_t s = fabs(par[2]);// sigma; Gaussian sigma
1350   Double_t w = fabs(par[3]);// Width; BW width
1351
1352   Double_t arg = x[0] - mean;
1353   Double_t _invRootPi = 1./sqrt(atan2(0.,-1.));
1354
1355   // return constant for zero width and sigma
1356   if (s==0. && w==0.) return 1.;
1357
1358   // Breit-Wigner for zero sigma
1359   if (s==0.) return (Norm*1./(arg*arg+0.25*w*w));
1360   Double_t coef= -0.5/(s*s);
1361
1362   // Gauss for zero width
1363   if (w==0.) return Norm*exp(coef*arg*arg);
1364
1365   // actual Voigtian for non-trivial width and sigma
1366   Double_t c = 1./(sqrt(2.)*s);
1367   Double_t a = 0.5*c*w;
1368   Double_t u = c*arg;
1369   RooComplex z(u,a);
1370   RooComplex v(0.);
1371
1372
1373   v = RooMath::ComplexErrFunc(z);
1374
1375   return Norm*c*_invRootPi*v.re();
1376
1377 }
1378 //________________________________________________________________________
1379  double VoigtplusPol(double *x, double *par){
1380    return Voigt(x,par) + PolFunction(x,&par[4]);
1381 }