[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / SpectraV0CutVariations.C

#if !defined(__CINT__) || defined(__MAKECINT__)
   #include <TMath.h>
   #include <TROOT.h>
   #include <Riostream.h>
   #include <TCanvas.h>
   #include <TStyle.h>
   #include <TString.h>
   #include <TLegend.h>
   #include <TFile.h>
   #include <TList.h>
   #include <TH1F.h>
   #include <TH1D.h>
   #include <TF2.h>
   #include <TFitResult.h>
   #include <TFitResultPtr.h>
   #include <TH2F.h>
   #include <TH3F.h>
#endif

TFile *fAss;
TFile *fGen;
TFile *fRaw;

void Generated();
void Corrections(Float_t cMin, Float_t cMax, TString centr);
void RawYields(Float_t cMin, Float_t cMax, TString centr);
void Spectra(TString centr);

void SpectraV0CutVariations(Float_t cMin, Float_t cMax, TString centr) {

  fRaw=TFile::Open((centr+"/AliV0CutVariations.root").Data());
  fAss=TFile::Open((centr+"/AliV0CutVariationsMC.root").Data());
  fGen=TFile::Open("Generated.root");
  if (!fGen) {
     Generated();
     fGen=TFile::Open("Generated.root");
  }

  Corrections(cMin,cMax,centr);
  RawYields(cMin,cMax,centr);
  Spectra(centr);

  fAss->Close();
  fGen->Close();
  fRaw->Close();
}

TH1 *GetEfficiency(Float_t, Float_t, const Char_t *, const Char_t *);
TH1 *GetFeedDown(Float_t, Float_t, TString);

void Corrections(Float_t cmin, Float_t cmax, TString centr) {

  TString name;

  TH1 *effK0s=
  GetEfficiency(cmin, cmax, "fK0sAs","f3dHistPrimRawPtVsYVsMultK0Short");
  name="eff_K0s_";
  name+=centr;
  effK0s->SetName(name.Data());

  TH1 *effLambda=
  GetEfficiency(cmin, cmax, "fLambdaAs","f3dHistPrimRawPtVsYVsMultLambda");
  name="eff_Lambda_";
  name+=centr;
  effLambda->SetName(name.Data());
    TH1 *fdLambda=GetFeedDown(cmin, cmax, centr);
    name="fd_Lambda_";
    name+=centr;
    fdLambda->SetName(name.Data());


  TH1 *effLambdaBar=
  GetEfficiency(cmin,cmax,"fLambdaBarAs","f3dHistPrimRawPtVsYVsMultAntiLambda");
  name="eff_LambdaBar_";
  name+=centr;
  effLambdaBar->SetName(name.Data());

  TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
    effK0s->Write();
    effLambda->Write();    fdLambda->Write();
    effLambdaBar->Write();
  f->Close();
}

TH1 *GetEfficiency(Float_t cMin, Float_t cMax, 
		   const Char_t *chis, const Char_t *znam) {
  Double_t xbins[]={
   0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
   1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
   2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
   4.5,5.0,5.5,6.5,8.0,10.0,12.0
  };
  Int_t nb=sizeof(xbins)/sizeof(Double_t);
  Int_t nb1=nb-1;

  // Numerator
  fAss->cd();
  TH2F *f2d=(TH2F*)gDirectory->Get(chis); f2d->Sumw2();
  TH1D *hAs=f2d->ProjectionX("hAs",0,-1,"e"); 

  // Denominator
  fGen->cd();
  TH3F *f3d = (TH3F*)gDirectory->Get(znam);
  f3d->Sumw2();

  TH1D *fpt = f3d->ProjectionX("fpt",
                  f3d->GetYaxis()->FindBin(-0.5+1e-2),
                  f3d->GetYaxis()->FindBin(+0.5-1e-2),
                  f3d->GetZaxis()->FindBin(cMin+1e-2),
                  f3d->GetZaxis()->FindBin(cMax-1e-2)
              );
   TH1 *hMc=fpt->Rebin(nb1,"hMc",xbins);
  
  //Efficiency 
  TH1 *eff = (TH1*)hAs->Clone();
  eff->Divide(hAs,hMc,1,1,"b");

  return eff;
}

Bool_t GetBinContentError(const TH1 *h, Double_t x, Double_t &c, Double_t &e) {
   Int_t i1=h->GetXaxis()->FindFixBin(x);
   if (i1 <=0) return kFALSE;
   Int_t nb=h->GetNbinsX();
   if (i1 >nb) return kFALSE;

   Double_t x1=h->GetBinCenter(i1);

   if (TMath::Abs(x1-x) < 1e-13) {
     c=h->GetBinContent(i1);
     e=h->GetBinError(i1);
     return kTRUE;
   }

   Int_t i2 = (x1 < x) ? i1+1 : i1-1;
   if (i2 <=0) i2=i1+1;
   if (i2 >nb) i2=i1-1;
   Double_t x2=h->GetBinCenter(i2);

   Double_t c1=h->GetBinContent(i1);
   Double_t c2=h->GetBinContent(i2);
   c = c1 + (x - x1)*(c2 - c1)/(x2 - x1);
   e = 0.5*(h->GetBinError(i1) + h->GetBinError(i2));

   return kTRUE;
}

TH1 *GetFeedDown(Float_t cMin, Float_t cMax, TString cent) {
  // Get the FD matrix
  fAss->cd();
  TH3F *f3d = (TH3F*)gDirectory->Get("fLambdaFromXi");
  f3d->Sumw2();
  TH2D *fdMat = (TH2D*)f3d->Project3D("zxe");
  const Int_t nbx=fdMat->GetNbinsX();
  const Int_t nby=fdMat->GetNbinsY();
  TAxis *xiPtAxis=fdMat->GetYaxis();

  // Re-normalise the FD matrix with MC Xi spectrum
  fGen->cd();
  f3d = (TH3F*)gDirectory->Get("f3dHistGenPtVsYVsMultXiMinus");
  f3d->Sumw2();

  TH1D *hMcXi = f3d->ProjectionX("fpt",
                  f3d->GetYaxis()->FindBin(-0.5+1e-2),
                  f3d->GetYaxis()->FindBin(+0.5-1e-2),
                  f3d->GetZaxis()->FindBin(cMin+1e-2),
                  f3d->GetZaxis()->FindBin(cMax-1e-2)
              );

  TH2D *h2McXi=new TH2D(*fdMat); 
  h2McXi->Reset();
  h2McXi->Sumw2();

  for (Int_t i=1; i<=nbx; i++) {
      for (Int_t j=1; j<=nby; j++) {
	  Double_t pt=xiPtAxis->GetBinCenter(j);
          Double_t c=0.,e=0.;
          if (!GetBinContentError(hMcXi,pt,c,e)) continue;
          h2McXi->SetBinContent(i,j,c);
          h2McXi->SetBinError(i,j,e);
      }
  }
  fdMat->Divide(h2McXi);
  delete h2McXi;


  // Multiply the re-normalised matrix with the real Xi spectrum
  TString fileName("systcorrectedpt_cent");
     //*** "Dictionary" ***
     if (cent=="0005") cent="0010";
     if (cent=="6080") cent="6090";
     if (cent=="8090") cent="6090";
  TFile *f=TFile::Open((fileName+cent+".root").Data());
     TH1F *hReXi=(TH1F*)gDirectory->Get("correctedpt_0");

     TH2D *h2ReXi=new TH2D(*fdMat); 
     h2ReXi->Reset();
     h2ReXi->Sumw2();

     for (Int_t i=1; i<=nbx; i++) {
         for (Int_t j=1; j<=nby; j++) {
	    Double_t pt=xiPtAxis->GetBinCenter(j);
            Double_t c=0.,e=0.;
            if (!GetBinContentError(hReXi,pt,c,e)) continue;
            h2ReXi->SetBinContent(i,j,c);
            h2ReXi->SetBinError(i,j,e);
         }
     }
     fdMat->Multiply(h2ReXi);
     delete h2ReXi;
     f->Close();

  TH1 *fd=fdMat->ProjectionX("_px",0,-1,"e");
  fd->Scale(0.1,"width");
  return fd;
}


void RawYields(Float_t cMin, Float_t cMax, TString centr) {
  TString name;
  TH2F *f2d=0;

  //+++ Number of events for normalisation
  TFile *file=TFile::Open("LHC10h_pass2/Merged.root");
    TList *v0list=(TList *)gFile->Get("PWGLFExtractV0_PP/clistV0");
    TH1F *fMult=(TH1F*)v0list->FindObject("fHistMultiplicity");
    Int_t i1=fMult->GetXaxis()->FindBin(cMin+1e-2);
    Int_t i2=fMult->GetXaxis()->FindBin(cMax-1e-2);
    Float_t nEvents=fMult->Integral(i1,i2);
  file->Close();

  fRaw->cd();

  name="raw_K0s_";
  name+=centr;  
  f2d=(TH2F*)gDirectory->Get("fK0sSi"); f2d->Sumw2();
  TH1D *rawK0s=f2d->ProjectionX(name,0,-1,"e");
  rawK0s->Scale(1/nEvents,"width");
 
  name="raw_Lambda_";
  name+=centr;  
  f2d=(TH2F*)gDirectory->Get("fLambdaSi"); f2d->Sumw2();
  TH1D *rawLambda=f2d->ProjectionX(name,0,-1,"e");
  rawLambda->Scale(1/nEvents,"width");
 
  name="raw_LambdaBar_";
  name+=centr;  
  f2d=(TH2F*)gDirectory->Get("fLambdaBarSi"); f2d->Sumw2();
  TH1D *rawLambdaBar=f2d->ProjectionX(name,0,-1,"e");
  rawLambdaBar->Scale(1/nEvents,"width");
 
  TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
    rawK0s->Write();
    rawLambda->Write();
    rawLambdaBar->Write();
  f->Close();
}

Double_t fd(Double_t x) {
  //Effective FD correction
   return 0.219901 + 0.0322588*x - 0.0173934*x*x + 0.00179039*x*x*x - 5.67881e-05*x*x*x*x;
}

void Spectra(TString centr) {
  TString name;
  TH1 *eff=0;
  TH1D *raw=0;
  TH1D *spe=0;

  TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
    name="eff_K0s_";
    name+=centr;
    eff = (TH1*)gDirectory->Get(name.Data());
    name="raw_K0s_";
    name+=centr;
    raw = (TH1D*)gDirectory->Get(name.Data());
    spe = new TH1D(*raw);
    spe->Divide(eff);
    name="K0s_";
    name+=centr;
    spe->SetName(name.Data());
    spe->Write();  

    name="eff_Lambda_";
    name+=centr;
    eff = (TH1*)gDirectory->Get(name.Data());
    name="raw_Lambda_";
    name+=centr;
    raw = (TH1D*)gDirectory->Get(name.Data());
    spe = new TH1D(*raw);
    for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
        Double_t pt=spe->GetBinCenter(i);
        Double_t c=spe->GetBinContent(i);
        c -= (c*fd(pt));
        spe->SetBinContent(i,c);
    }
    spe->Divide(eff);
    name="Lambda_";
    name+=centr;
    spe->SetName(name.Data());
    spe->Write();  

    name="eff_LambdaBar_";
    name+=centr;
    eff = (TH1*)gDirectory->Get(name.Data());
    name="raw_LambdaBar_";
    name+=centr;
    raw = (TH1D*)gDirectory->Get(name.Data());
    spe = new TH1D(*raw);
    for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
        Double_t pt=spe->GetBinCenter(i);
        Double_t c=spe->GetBinContent(i);
        c -= (c*fd(pt));
        spe->SetBinContent(i,c);
    }
    spe->Divide(eff);
    name="LambdaBar_";
    name+=centr;
    spe->SetName(name.Data());
    spe->Write();  
  f->Close();
}

void Generated() {
  TList *v0listMC=0;
  TH3F *h3=0;
  //
  TFile::Open("LHC11a10b_plus/Merged.root");  
  v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");

  TH3F *k0s = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
  TH3F *k0s_nonInj = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");

  TH3F *lam = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
  TH3F *lam_nonInj = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");

  TH3F *alam = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
  TH3F *alam_nonInj = 
    (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");

  TH3F *xiMinus = 
    (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
  TH3F *xiPlus = 
    (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");


  //    
  TFile::Open("LHC11a10b_bis/Merged.root");  
  v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");

  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
  k0s->Add(h3);
  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
  k0s_nonInj->Add(h3); 

  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
  lam->Add(h3); 
  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
  lam_nonInj->Add(h3); 

  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
  alam->Add(h3); 
  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
  alam_nonInj->Add(h3); 

  h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
  xiMinus->Add(h3);
  h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
  xiPlus->Add(h3); 


  //
  TFile::Open("LHC11a10a_bis/Merged.root");
  v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");

  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
  k0s_nonInj->Add(h3); 
  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
  lam_nonInj->Add(h3); 
  h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
  alam_nonInj->Add(h3); 
  h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
  xiMinus->Add(h3);
  h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
  xiPlus->Add(h3); 

  TFile *f=TFile::Open("Generated.root","new");
  k0s->Write(); k0s_nonInj->Write();
  lam->Write(); lam_nonInj->Write();
  alam->Write(); alam_nonInj->Write();
  xiMinus->Write();
  xiPlus->Write();
  f->Close();
}

extern TStyle *gStyle;

void SetOptions(TH1 *h, const Char_t *t, Int_t c, Int_t m, Float_t s) {
   h->SetLineColor(c);
   h->SetMarkerColor(c); 
   h->SetMarkerStyle(m);
   h->SetMarkerSize(s);
   h->SetTitle(t);
   TAxis *x=h->GetXaxis();
   x->SetTitle("p_{T} (GeV/c)");
   x->SetTitleSize(0.05);   
   x->SetLabelSize(0.05);   
   TAxis *y=h->GetYaxis();
   y->SetTitle("1/N_{ev}dN/dp_{T}/dy (GeV/c)^{-1}");   
   y->SetTitleSize(0.05);   
   y->SetLabelSize(0.05);   
}
void DrawSpectra() {
  TH1 *h=0;
  TString name;

  const Char_t *title[]={
    "0-5 %",
    "0-10 %",
    "10-20 %",
    "20-40 %",
    "40-60 %",
    "60-80 %",
    "80-90 %"
  };
  const Int_t nCent=sizeof(title)/sizeof(const Char_t *);
  TString centr[nCent]={"0005","0010","1020","2040","4060","6080","8090"};
  const Int_t   colour[nCent]={2,   3,   635, 419, 4 , 6,  1  };
  const Int_t   marker[nCent]={22,  29,  34,  21,  23, 33, 20 };
  const Float_t masize[nCent]={1.6, 1.7, 1.3, 1.3, 1.6,2,  1.3};

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetLegendFillColor(0);

  TCanvas *kc=new TCanvas(); kc->SetLogy();
  TCanvas *lc=new TCanvas(); lc->SetLogy();
  TCanvas *rc=new TCanvas();

  TFile::Open("SpectraV0CutVariations.root");
  
  for (Int_t i=0; i<nCent; i++) {
    name="Lambda_";
    h=(TH1*)gDirectory->Get((name+centr[i]).Data());
    SetOptions(h,title[i],colour[i],marker[i],masize[i]);
    lc->cd();
    if (i==0) h->Draw(); else h->Draw("same");

    TH1F *r=(TH1F*)h->Clone();
    r->GetYaxis()->SetTitle("#Lambda/K^{0}_{S}");

    name="K0s_";
    h=(TH1*)gDirectory->Get((name+centr[i]).Data());
    SetOptions(h,title[i],colour[i],marker[i],masize[i]);
    kc->cd();
    if (i==0) h->Draw(); else h->Draw("same");

    r->Divide(h);
    rc->cd();
    if (i==0) r->Draw(); else r->Draw("same");
  }

  TLegend *leg=0;
  kc->cd();
  leg=kc->BuildLegend(0.659,0.320,0.881,0.875,"K^{0}_{S} spectra:");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);

  lc->cd();
  leg=lc->BuildLegend(0.659,0.320,0.881,0.875,"#Lambda spectra:");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);

  rc->cd();
  leg=rc->BuildLegend(0.659,0.320,0.881,0.875,"#Lambda/K^{0}_{S} ratios:");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
}

void Loop() {
  const Char_t *cent[]={"0005","0010","1020","2040","4060","6080","8090"};
  Int_t nCent=sizeof(cent)/sizeof(const Char_t *);
  const Int_t cmin[]={0, 0,  10, 20, 40, 60, 80};
  const Int_t cmax[]={5, 10, 20, 40, 60, 80, 90};

  for (Int_t i=0; i<nCent; i++) {
    cout<<cent[i]<<endl;
    SpectraV0CutVariations(cmin[i],cmax[i],cent[i]);
  }
}
Commit	Line	Data
2c8e2a0d	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TMath.h>
	3	#include <TROOT.h>
	4	#include <Riostream.h>
	5	#include <TCanvas.h>
d5fa2d40	6	#include <TStyle.h>
2c8e2a0d	7	#include <TString.h>
d5fa2d40	8	#include <TLegend.h>
2c8e2a0d	9	#include <TFile.h>
	10	#include <TList.h>
	11	#include <TH1F.h>
	12	#include <TH1D.h>
	13	#include <TF2.h>
	14	#include <TFitResult.h>
	15	#include <TFitResultPtr.h>
	16	#include <TH2F.h>
	17	#include <TH3F.h>
	18	#endif
	19
2c8e2a0d	20	TFile *fAss;
	21	TFile *fGen;
	22	TFile *fRaw;
	23
	24	void Generated();
ebebf970	25	void Corrections(Float_t cMin, Float_t cMax, TString centr);
2c8e2a0d	26	void RawYields(Float_t cMin, Float_t cMax, TString centr);
	27	void Spectra(TString centr);
	28
	29	void SpectraV0CutVariations(Float_t cMin, Float_t cMax, TString centr) {
	30
	31	fRaw=TFile::Open((centr+"/AliV0CutVariations.root").Data());
	32	fAss=TFile::Open((centr+"/AliV0CutVariationsMC.root").Data());
	33	fGen=TFile::Open("Generated.root");
	34	if (!fGen) {
	35	Generated();
	36	fGen=TFile::Open("Generated.root");
	37	}
	38
ebebf970	39	Corrections(cMin,cMax,centr);
2c8e2a0d	40	RawYields(cMin,cMax,centr);
	41	Spectra(centr);
	42
	43	fAss->Close();
	44	fGen->Close();
	45	fRaw->Close();
	46	}
	47
	48	TH1 GetEfficiency(Float_t, Float_t, const Char_t , const Char_t *);
ebebf970	49	TH1 *GetFeedDown(Float_t, Float_t, TString);
2c8e2a0d	50
ebebf970	51	void Corrections(Float_t cmin, Float_t cmax, TString centr) {
2c8e2a0d	52
	53	TString name;
	54
	55	TH1 *effK0s=
	56	GetEfficiency(cmin, cmax, "fK0sAs","f3dHistPrimRawPtVsYVsMultK0Short");
	57	name="eff_K0s_";
	58	name+=centr;
	59	effK0s->SetName(name.Data());
	60
	61	TH1 *effLambda=
	62	GetEfficiency(cmin, cmax, "fLambdaAs","f3dHistPrimRawPtVsYVsMultLambda");
	63	name="eff_Lambda_";
	64	name+=centr;
	65	effLambda->SetName(name.Data());
ebebf970	66	TH1 *fdLambda=GetFeedDown(cmin, cmax, centr);
	67	name="fd_Lambda_";
	68	name+=centr;
	69	fdLambda->SetName(name.Data());
	70
2c8e2a0d	71
	72	TH1 *effLambdaBar=
	73	GetEfficiency(cmin,cmax,"fLambdaBarAs","f3dHistPrimRawPtVsYVsMultAntiLambda");
	74	name="eff_LambdaBar_";
	75	name+=centr;
	76	effLambdaBar->SetName(name.Data());
	77
ea54d2e3	78	TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d	79	effK0s->Write();
ebebf970	80	effLambda->Write(); fdLambda->Write();
2c8e2a0d	81	effLambdaBar->Write();
	82	f->Close();
	83	}
	84
	85	TH1 *GetEfficiency(Float_t cMin, Float_t cMax,
	86	const Char_t chis, const Char_t znam) {
9f1ae35f	87	Double_t xbins[]={
	88	0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
	89	1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
	90	2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
	91	4.5,5.0,5.5,6.5,8.0,10.0,12.0
	92	};
	93	Int_t nb=sizeof(xbins)/sizeof(Double_t);
	94	Int_t nb1=nb-1;
	95
2c8e2a0d	96	// Numerator
	97	fAss->cd();
	98	TH2F f2d=(TH2F)gDirectory->Get(chis); f2d->Sumw2();
	99	TH1D *hAs=f2d->ProjectionX("hAs",0,-1,"e");
	100
	101	// Denominator
	102	fGen->cd();
	103	TH3F f3d = (TH3F)gDirectory->Get(znam);
	104	f3d->Sumw2();
	105
	106	TH1D *fpt = f3d->ProjectionX("fpt",
	107	f3d->GetYaxis()->FindBin(-0.5+1e-2),
	108	f3d->GetYaxis()->FindBin(+0.5-1e-2),
d58451d7	109	f3d->GetZaxis()->FindBin(cMin+1e-2),
d58451d7	110	f3d->GetZaxis()->FindBin(cMax-1e-2)
2c8e2a0d	111	);
	112	TH1 *hMc=fpt->Rebin(nb1,"hMc",xbins);
	113
	114	//Efficiency
	115	TH1 eff = (TH1)hAs->Clone();
	116	eff->Divide(hAs,hMc,1,1,"b");
	117
	118	return eff;
	119	}
	120
ebebf970	121	Bool_t GetBinContentError(const TH1 *h, Double_t x, Double_t &c, Double_t &e) {
	122	Int_t i1=h->GetXaxis()->FindFixBin(x);
	123	if (i1 <=0) return kFALSE;
	124	Int_t nb=h->GetNbinsX();
	125	if (i1 >nb) return kFALSE;
	126
	127	Double_t x1=h->GetBinCenter(i1);
	128
	129	if (TMath::Abs(x1-x) < 1e-13) {
	130	c=h->GetBinContent(i1);
	131	e=h->GetBinError(i1);
	132	return kTRUE;
	133	}
	134
	135	Int_t i2 = (x1 < x) ? i1+1 : i1-1;
	136	if (i2 <=0) i2=i1+1;
	137	if (i2 >nb) i2=i1-1;
	138	Double_t x2=h->GetBinCenter(i2);
	139
	140	Double_t c1=h->GetBinContent(i1);
	141	Double_t c2=h->GetBinContent(i2);
	142	c = c1 + (x - x1)*(c2 - c1)/(x2 - x1);
	143	e = 0.5*(h->GetBinError(i1) + h->GetBinError(i2));
	144
	145	return kTRUE;
	146	}
	147
	148	TH1 *GetFeedDown(Float_t cMin, Float_t cMax, TString cent) {
	149	// Get the FD matrix
	150	fAss->cd();
	151	TH3F f3d = (TH3F)gDirectory->Get("fLambdaFromXi");
	152	f3d->Sumw2();
	153	TH2D fdMat = (TH2D)f3d->Project3D("zxe");
	154	const Int_t nbx=fdMat->GetNbinsX();
	155	const Int_t nby=fdMat->GetNbinsY();
	156	TAxis *xiPtAxis=fdMat->GetYaxis();
	157
	158	// Re-normalise the FD matrix with MC Xi spectrum
	159	fGen->cd();
	160	f3d = (TH3F*)gDirectory->Get("f3dHistGenPtVsYVsMultXiMinus");
	161	f3d->Sumw2();
	162
	163	TH1D *hMcXi = f3d->ProjectionX("fpt",
	164	f3d->GetYaxis()->FindBin(-0.5+1e-2),
	165	f3d->GetYaxis()->FindBin(+0.5-1e-2),
d58451d7	166	f3d->GetZaxis()->FindBin(cMin+1e-2),
d58451d7	167	f3d->GetZaxis()->FindBin(cMax-1e-2)
ebebf970	168	);
	169
	170	TH2D h2McXi=new TH2D(fdMat);
	171	h2McXi->Reset();
	172	h2McXi->Sumw2();
	173
	174	for (Int_t i=1; i<=nbx; i++) {
	175	for (Int_t j=1; j<=nby; j++) {
	176	Double_t pt=xiPtAxis->GetBinCenter(j);
	177	Double_t c=0.,e=0.;
	178	if (!GetBinContentError(hMcXi,pt,c,e)) continue;
	179	h2McXi->SetBinContent(i,j,c);
	180	h2McXi->SetBinError(i,j,e);
	181	}
	182	}
	183	fdMat->Divide(h2McXi);
	184	delete h2McXi;
	185
	186
	187	// Multiply the re-normalised matrix with the real Xi spectrum
	188	TString fileName("systcorrectedpt_cent");
	189	//* "Dictionary" *
	190	if (cent=="0005") cent="0010";
	191	if (cent=="6080") cent="6090";
	192	if (cent=="8090") cent="6090";
	193	TFile *f=TFile::Open((fileName+cent+".root").Data());
	194	TH1F hReXi=(TH1F)gDirectory->Get("correctedpt_0");
	195
	196	TH2D h2ReXi=new TH2D(fdMat);
	197	h2ReXi->Reset();
	198	h2ReXi->Sumw2();
	199
	200	for (Int_t i=1; i<=nbx; i++) {
	201	for (Int_t j=1; j<=nby; j++) {
	202	Double_t pt=xiPtAxis->GetBinCenter(j);
	203	Double_t c=0.,e=0.;
	204	if (!GetBinContentError(hReXi,pt,c,e)) continue;
	205	h2ReXi->SetBinContent(i,j,c);
	206	h2ReXi->SetBinError(i,j,e);
	207	}
	208	}
	209	fdMat->Multiply(h2ReXi);
	210	delete h2ReXi;
	211	f->Close();
	212
	213	TH1 *fd=fdMat->ProjectionX("_px",0,-1,"e");
	214	fd->Scale(0.1,"width");
	215	return fd;
	216	}
	217
	218
2c8e2a0d	219	void RawYields(Float_t cMin, Float_t cMax, TString centr) {
	220	TString name;
	221	TH2F *f2d=0;
	222
	223	//+++ Number of events for normalisation
ea54d2e3	224	TFile *file=TFile::Open("LHC10h_pass2/Merged.root");
2c8e2a0d	225	TList v0list=(TList )gFile->Get("PWGLFExtractV0_PP/clistV0");
	226	TH1F fMult=(TH1F)v0list->FindObject("fHistMultiplicity");
	227	Int_t i1=fMult->GetXaxis()->FindBin(cMin+1e-2);
d58451d7	228	Int_t i2=fMult->GetXaxis()->FindBin(cMax-1e-2);
2c8e2a0d	229	Float_t nEvents=fMult->Integral(i1,i2);
	230	file->Close();
	231
	232	fRaw->cd();
	233
	234	name="raw_K0s_";
	235	name+=centr;
	236	f2d=(TH2F*)gDirectory->Get("fK0sSi"); f2d->Sumw2();
	237	TH1D *rawK0s=f2d->ProjectionX(name,0,-1,"e");
	238	rawK0s->Scale(1/nEvents,"width");
	239
	240	name="raw_Lambda_";
	241	name+=centr;
	242	f2d=(TH2F*)gDirectory->Get("fLambdaSi"); f2d->Sumw2();
	243	TH1D *rawLambda=f2d->ProjectionX(name,0,-1,"e");
	244	rawLambda->Scale(1/nEvents,"width");
	245
	246	name="raw_LambdaBar_";
	247	name+=centr;
	248	f2d=(TH2F*)gDirectory->Get("fLambdaBarSi"); f2d->Sumw2();
	249	TH1D *rawLambdaBar=f2d->ProjectionX(name,0,-1,"e");
	250	rawLambdaBar->Scale(1/nEvents,"width");
	251
ea54d2e3	252	TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d	253	rawK0s->Write();
	254	rawLambda->Write();
	255	rawLambdaBar->Write();
	256	f->Close();
	257	}
	258
8fb8d5db	259	Double_t fd(Double_t x) {
8fb8d5db	260	//Effective FD correction
b4fd0d82	261	return 0.219901 + 0.0322588x - 0.0173934xx + 0.00179039xxx - 5.67881e-05xxxx;
8fb8d5db	262	}
8fb8d5db	263
2c8e2a0d	264	void Spectra(TString centr) {
	265	TString name;
	266	TH1 *eff=0;
	267	TH1D *raw=0;
	268	TH1D *spe=0;
	269
ea54d2e3	270	TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
2c8e2a0d	271	name="eff_K0s_";
	272	name+=centr;
	273	eff = (TH1*)gDirectory->Get(name.Data());
	274	name="raw_K0s_";
	275	name+=centr;
	276	raw = (TH1D*)gDirectory->Get(name.Data());
	277	spe = new TH1D(*raw);
	278	spe->Divide(eff);
	279	name="K0s_";
	280	name+=centr;
	281	spe->SetName(name.Data());
	282	spe->Write();
	283
	284	name="eff_Lambda_";
	285	name+=centr;
	286	eff = (TH1*)gDirectory->Get(name.Data());
	287	name="raw_Lambda_";
	288	name+=centr;
	289	raw = (TH1D*)gDirectory->Get(name.Data());
	290	spe = new TH1D(*raw);
8fb8d5db	291	for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
	292	Double_t pt=spe->GetBinCenter(i);
	293	Double_t c=spe->GetBinContent(i);
	294	c -= (c*fd(pt));
	295	spe->SetBinContent(i,c);
	296	}
2c8e2a0d	297	spe->Divide(eff);
	298	name="Lambda_";
	299	name+=centr;
	300	spe->SetName(name.Data());
	301	spe->Write();
	302
	303	name="eff_LambdaBar_";
	304	name+=centr;
	305	eff = (TH1*)gDirectory->Get(name.Data());
	306	name="raw_LambdaBar_";
	307	name+=centr;
	308	raw = (TH1D*)gDirectory->Get(name.Data());
	309	spe = new TH1D(*raw);
8fb8d5db	310	for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
	311	Double_t pt=spe->GetBinCenter(i);
	312	Double_t c=spe->GetBinContent(i);
	313	c -= (c*fd(pt));
	314	spe->SetBinContent(i,c);
	315	}
2c8e2a0d	316	spe->Divide(eff);
	317	name="LambdaBar_";
	318	name+=centr;
	319	spe->SetName(name.Data());
	320	spe->Write();
	321	f->Close();
	322	}
	323
	324	void Generated() {
	325	TList *v0listMC=0;
	326	TH3F *h3=0;
	327	//
	328	TFile::Open("LHC11a10b_plus/Merged.root");
	329	v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
	330
	331	TH3F *k0s =
	332	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
	333	TH3F *k0s_nonInj =
	334	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
	335
	336	TH3F *lam =
	337	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
	338	TH3F *lam_nonInj =
	339	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
	340
	341	TH3F *alam =
	342	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
	343	TH3F *alam_nonInj =
	344	(TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
	345
ebebf970	346	TH3F *xiMinus =
	347	(TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
	348	TH3F *xiPlus =
	349	(TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
	350
2c8e2a0d	351
	352	//
	353	TFile::Open("LHC11a10b_bis/Merged.root");
	354	v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
	355
	356	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
	357	k0s->Add(h3);
	358	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
	359	k0s_nonInj->Add(h3);
	360
	361	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
	362	lam->Add(h3);
	363	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
	364	lam_nonInj->Add(h3);
	365
	366	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
	367	alam->Add(h3);
	368	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
	369	alam_nonInj->Add(h3);
	370
ebebf970	371	h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
	372	xiMinus->Add(h3);
	373	h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
	374	xiPlus->Add(h3);
	375
	376
2c8e2a0d	377	//
2c8e2a0d	378	TFile::Open("LHC11a10a_bis/Merged.root");
7e9df831	379	v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
7e9df831	380
2c8e2a0d	381	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
	382	k0s_nonInj->Add(h3);
	383	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
	384	lam_nonInj->Add(h3);
	385	h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
	386	alam_nonInj->Add(h3);
ebebf970	387	h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
	388	xiMinus->Add(h3);
	389	h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
	390	xiPlus->Add(h3);
2c8e2a0d	391
	392	TFile *f=TFile::Open("Generated.root","new");
	393	k0s->Write(); k0s_nonInj->Write();
	394	lam->Write(); lam_nonInj->Write();
	395	alam->Write(); alam_nonInj->Write();
ebebf970	396	xiMinus->Write();
ebebf970	397	xiPlus->Write();
2c8e2a0d	398	f->Close();
2c8e2a0d	399	}
ea54d2e3	400
d5fa2d40	401	extern TStyle *gStyle;
	402
	403	void SetOptions(TH1 h, const Char_t t, Int_t c, Int_t m, Float_t s) {
	404	h->SetLineColor(c);
	405	h->SetMarkerColor(c);
	406	h->SetMarkerStyle(m);
	407	h->SetMarkerSize(s);
	408	h->SetTitle(t);
	409	TAxis *x=h->GetXaxis();
	410	x->SetTitle("p_{T} (GeV/c)");
	411	x->SetTitleSize(0.05);
	412	x->SetLabelSize(0.05);
	413	TAxis *y=h->GetYaxis();
	414	y->SetTitle("1/N_{ev}dN/dp_{T}/dy (GeV/c)^{-1}");
	415	y->SetTitleSize(0.05);
	416	y->SetLabelSize(0.05);
	417	}
	418	void DrawSpectra() {
	419	TH1 *h=0;
	420	TString name;
	421
	422	const Char_t *title[]={
	423	"0-5 %",
	424	"0-10 %",
	425	"10-20 %",
	426	"20-40 %",
	427	"40-60 %",
	428	"60-80 %",
	429	"80-90 %"
	430	};
	431	const Int_t nCent=sizeof(title)/sizeof(const Char_t *);
	432	TString centr[nCent]={"0005","0010","1020","2040","4060","6080","8090"};
	433	const Int_t colour[nCent]={2, 3, 635, 419, 4 , 6, 1 };
	434	const Int_t marker[nCent]={22, 29, 34, 21, 23, 33, 20 };
	435	const Float_t masize[nCent]={1.6, 1.7, 1.3, 1.3, 1.6,2, 1.3};
	436
	437	gStyle->SetOptStat(0);
	438	gStyle->SetOptTitle(0);
	439	gStyle->SetLegendFillColor(0);
	440
	441	TCanvas *kc=new TCanvas(); kc->SetLogy();
	442	TCanvas *lc=new TCanvas(); lc->SetLogy();
c4209f2b	443	TCanvas *rc=new TCanvas();
d5fa2d40	444
	445	TFile::Open("SpectraV0CutVariations.root");
	446
	447	for (Int_t i=0; i<nCent; i++) {
c4209f2b	448	name="Lambda_";
d5fa2d40	449	h=(TH1*)gDirectory->Get((name+centr[i]).Data());
d5fa2d40	450	SetOptions(h,title[i],colour[i],marker[i],masize[i]);
c4209f2b	451	lc->cd();
d5fa2d40	452	if (i==0) h->Draw(); else h->Draw("same");
d5fa2d40	453
c4209f2b	454	TH1F r=(TH1F)h->Clone();
	455	r->GetYaxis()->SetTitle("#Lambda/K^{0}_{S}");
	456
	457	name="K0s_";
d5fa2d40	458	h=(TH1*)gDirectory->Get((name+centr[i]).Data());
d5fa2d40	459	SetOptions(h,title[i],colour[i],marker[i],masize[i]);
c4209f2b	460	kc->cd();
d5fa2d40	461	if (i==0) h->Draw(); else h->Draw("same");
c4209f2b	462
	463	r->Divide(h);
	464	rc->cd();
	465	if (i==0) r->Draw(); else r->Draw("same");
d5fa2d40	466	}
d5fa2d40	467
c4209f2b	468	TLegend *leg=0;
d5fa2d40	469	kc->cd();
c4209f2b	470	leg=kc->BuildLegend(0.659,0.320,0.881,0.875,"K^{0}_{S} spectra:");
	471	leg->SetBorderSize(0);
	472	leg->SetFillColor(0);
	473
d5fa2d40	474	lc->cd();
c4209f2b	475	leg=lc->BuildLegend(0.659,0.320,0.881,0.875,"#Lambda spectra:");
	476	leg->SetBorderSize(0);
	477	leg->SetFillColor(0);
	478
	479	rc->cd();
	480	leg=rc->BuildLegend(0.659,0.320,0.881,0.875,"#Lambda/K^{0}_{S} ratios:");
	481	leg->SetBorderSize(0);
	482	leg->SetFillColor(0);
d5fa2d40	483	}
0fc5ae6c	484
	485	void Loop() {
	486	const Char_t *cent[]={"0005","0010","1020","2040","4060","6080","8090"};
	487	Int_t nCent=sizeof(cent)/sizeof(const Char_t *);
	488	const Int_t cmin[]={0, 0, 10, 20, 40, 60, 80};
	489	const Int_t cmax[]={5, 10, 20, 40, 60, 80, 90};
	490
	491	for (Int_t i=0; i<nCent; i++) {
	492	cout<<cent[i]<<endl;
	493	SpectraV0CutVariations(cmin[i],cmax[i],cent[i]);
	494	}
	495	}