[u/mrichter/AliRoot.git] / PWG3 / dielectron / macros / PlotDataResults.C

#include <TROOT.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TObjArray.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TAxis.h>
#include <TLatex.h>
#include <TFile.h>
#include <TString.h>
#include <TLegend.h>
#include <TList.h>

#include "AliDielectronSignalExt.h"
#include "AliDielectronCFdraw.h"
#include "AliDielectron.h"

AliDielectronSignalBase* GetSignalLS(AliDielectronCFdraw &d, Int_t step);
AliDielectronSignalBase* GetSignalRot(AliDielectronCFdraw &d, Int_t step);
void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd);
void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat=0x0, Bool_t save=kFALSE);
void FitMCshape(AliDielectronSignalBase *sig);

const char *mcLineShapeFile="$ALICE_ROOT/PWG3/dielectron/macros/mcMinv_LHC10e2.root";

//_______________________________________
void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t save=kFALSE)
{
  AliDielectronCFdraw d(filenameData);
  AliDielectronCFdraw dCorr("corrCont","corrCont");
  TString nameCorr(filenameMC);
  if (!nameCorr.IsNull()) d.SetCFContainers(nameCorr.Data());
  TFile f(filenameData);
  TH1 *hStats=(TH1*)f.Get("hEventStat");
  hStats->SetDirectory(0);
  f.Close();
  
  Int_t stepFirst=0, stepAny=1, stepTOFmix=2;
  
  gStyle->SetOptStat(0);
  //Set common Ranges
//   d.SetRangeUser("Leg1_Pt",0.8,1000.);
//   d.SetRangeUser("Leg2_Pt",0.8,1000.);
  d.SetRangeUser("Leg1_NclsTPC",140.,170.);
  d.SetRangeUser("Leg2_NclsTPC",140.,170.);
  d.SetRangeUser("Leg1_Pt",1.01,100000);
  d.SetRangeUser("Leg2_Pt",1.01,100000);
//   d.SetRangeUser("Leg1_TPC_nSigma_Electrons",-3,3);
//   d.SetRangeUser("Leg2_TPC_nSigma_Electrons",-3,3);
//   d.SetRangeUser("Leg1_TPC_nSigma_Pions",3,20);
//   d.SetRangeUser("Leg2_TPC_nSigma_Pions",3,20);
//   d.SetRangeUser("Leg1_TPC_nSigma_Protons",3,20);
//   d.SetRangeUser("Leg2_TPC_nSigma_Protons",3,20);
  
  d.SetRangeUser("M",0.5,5.);
  //============================
  //SPD first
  //
  
  //--- Like sign subtraction
  AliDielectronSignalBase *sigFirst=GetSignalLS(d,stepFirst);
  SetStyle(sigFirst,"ITS First - Like Sign subtraction");
//   DrawSpectra(sigFirst,"cFirst",hStats,save);
  //--- Rotation subtraction
  AliDielectronSignalBase *sigFirstRot=GetSignalRot(d,stepFirst);
  SetStyle(sigFirstRot,"ITS First - Track rotation subtraction");
//   DrawSpectra(sigFirstRot,"cFirstRot",hStats,save);
  
  //============================
  //SPD any
  //
  AliDielectronSignalBase *sigAny=GetSignalLS(d,stepAny);
  SetStyle(sigAny,"ITS First - Like Sign subtraction");
  DrawSpectra(sigAny,"cAny",hStats,save);
  //--- Rotation subtraction
  AliDielectronSignalBase *sigAnyRot=GetSignalRot(d,stepAny);
  SetStyle(sigAnyRot,"ITS First - Track rotation subtraction");
//   DrawSpectra(sigAnyRot,"cAnyRot",hStats,save);
  
  //=============================
  //TOF up to 1.2, parametrisation in TPC ele
  //
  AliDielectronSignalBase *sigTOFmix=GetSignalLS(d,stepTOFmix);
  SetStyle(sigTOFmix,"TOF + TPC - Like Sign subtraction");
//   DrawSpectra(sigTOFmix,"cTOFTPC",hStats,save);
  //--- Rotation subtraction
  AliDielectronSignalBase *sigTOFmixRot=GetSignalRot(d,stepTOFmix);
  SetStyle(sigTOFmixRot,"TOF + TPC - Track rotation subtraction");
//   DrawSpectra(sigTOFmixRot,"cTOFTPCrot",hStats,save);
  
  if (hStats) delete hStats;
}


//_______________________________________
AliDielectronSignalBase *GetSignalLS(AliDielectronCFdraw &d, Int_t step)
{
  //
  // Get Extracted signal from likesign method
  //
  
  TObjArray *arr=new TObjArray;
  arr->SetOwner();
  
  for (Int_t iType=0;iType<3;++iType){
    d.SetRangeUser("PairType",iType,iType);
    arr->AddAt(d.Project("M",step),iType);
  }
  
  AliDielectronSignalExt *sig=new AliDielectronSignalExt;
  sig->SetScaleRawToBackground(3.2,5.);
  sig->SetIntegralRange(2.92,3.15);
  sig->SetMethod(AliDielectronSignalBase::kLikeSign);
  sig->Process(arr);
  
  delete arr;
  return sig;
}

//_______________________________________
AliDielectronSignalBase *GetSignalRot(AliDielectronCFdraw &d, Int_t step)
{
  //
  // Get Extracted signal from likesign method
  //
  
  TObjArray *arr=new TObjArray;
  arr->SetOwner();
  
  Int_t iType=AliDielectron::kEv1PM;
  d.SetRangeUser("PairType",iType,iType);
  arr->AddAt(d.Project("M",step),iType);
  
  iType=AliDielectron::kEv1PMRot;
  d.SetRangeUser("PairType",iType,iType);
  arr->AddAt(d.Project("M",step),iType);
  
  AliDielectronSignalExt *sig=new AliDielectronSignalExt;
//   sig->SetScaleRawToBackground(3.2,4.);
  sig->SetIntegralRange(2.93,3.15);
  sig->SetMethod(AliDielectronSignalBase::kRotation);
  sig->Process(arr);
  
  
  delete arr;
  return sig;
}

//_______________________________________
void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd)
{
  //
  //
  //
  TH1 *hUS=sig->GetUnlikeSignHistogram();
  hUS->SetMarkerStyle(20);
  hUS->SetMarkerSize(0.7);
  hUS->SetMarkerColor(kRed);
  hUS->SetLineColor(kRed);
  hUS->SetStats(0);
  hUS->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
                     nameAdd,1000*hUS->GetXaxis()->GetBinWidth(1)));
  
  TH1* hBackground=sig->GetBackgroundHistogram();
  hBackground->SetMarkerStyle(24);
  hBackground->SetMarkerSize(0.7);
  hBackground->SetStats(0);
  hBackground->SetMarkerColor(kBlue);
  hBackground->SetLineColor(kBlue);
  hBackground->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
                             nameAdd,1000*hBackground->GetXaxis()->GetBinWidth(1)));
  
  TH1* hSignal=sig->GetSignalHistogram();
  hSignal->SetMarkerStyle(20);
  hSignal->SetMarkerSize(0.7);
  hSignal->SetMarkerColor(kRed);
  hSignal->SetLineColor(kRed);
  hSignal->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
                         nameAdd,1000*hSignal->GetXaxis()->GetBinWidth(1)));
  
  
}

//_______________________________________
void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat, Bool_t save)
{
  //
  //
  //
  
  
  FitMCshape(sig);
  
  gStyle->SetOptTitle(0);
  TCanvas *c=(TCanvas*)gROOT->FindObject(cname);
  if (!c) c=new TCanvas(cname,cname,400,500);
  c->SetTopMargin(0.04);  c->SetRightMargin(0.04);
  c->SetLeftMargin(0.13); c->SetBottomMargin(0.14);
  c->Clear();
  c->Divide(1,2,0,0);
  c->cd(1);
  gPad->SetTopMargin(0.01);
  gPad->SetRightMargin(0.01);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.001);
  /*
  gPad->SetGridx();
  gPad->SetGridy();
*/
  gPad->SetTickx();
  gPad->SetTicky();
  
  TH1 *hUS=sig->GetUnlikeSignHistogram();
  
  TH1* hBackground=sig->GetBackgroundHistogram();
  
  hUS->GetXaxis()->CenterTitle();
  hUS->GetXaxis()->SetTitleSize(0.07); hUS->GetXaxis()->SetLabelSize(0.059);
//   hUS->GetXaxis()->SetLabelOffset(1.);
  hUS->GetXaxis()->SetTitleOffset(1.1);
  hUS->GetXaxis()->SetNdivisions(510);
  hUS->GetYaxis()->CenterTitle();
  hUS->GetYaxis()->SetTitleSize(0.07); hUS->GetYaxis()->SetLabelSize(0.059);
  hUS->GetYaxis()->SetTitleOffset(1.);
  hUS->GetXaxis()->SetLabelFont(42); hUS->GetYaxis()->SetLabelFont(42);
  hUS->GetXaxis()->SetTitleFont(42); hUS->GetYaxis()->SetTitleFont(42);
  hUS->GetYaxis()->SetNdivisions(510);
  //hUS->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
  //hUS->SetLineWidth(.7);  
  hUS->SetMarkerSize(.8);
  hUS->SetMarkerColor(2); hUS->SetMarkerStyle(20); hUS->SetLineColor(2);
  hUS->Draw();
  hUS->SetMaximum(hUS->GetMaximum()*1.3);

  //hBackground->SetLineWidth(.4);  
  hBackground->SetMarkerSize(.7);
  hBackground->SetMarkerColor(4); hBackground->SetMarkerStyle(27); hBackground->SetLineColor(4);

  hBackground->Draw("same");
  
  TLatex *lat=new TLatex;
  lat->SetNDC(kTRUE);
  //  lat->DrawLatex(0.68, 0.67, "ALICE Performance");
  lat->SetTextColor(1);lat->SetTextFont(42);lat->SetTextSize(.045);
  
  Double_t sigN=sig->GetSignal();
  Double_t sigEr=sig->GetSignalError();
  Double_t sigS2B=sig->GetSB();
  Double_t sigS2Ber=sig->GetSBError();
  Double_t sigSignif= sig->GetSignificance();
  Double_t sigSignifEr= sig->GetSignificanceError();
  lat->DrawLatex(0.18, 0.92, Form("S: %3d#pm%4.1f, S/B: %3.1f#pm %4.2f, Signif.: %4.1f#pm%4.2f (%4.2f-%4.2f GeV) ",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,
                       hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMin())),
                       hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMax())+1)));
  
  TLegend *leg=new TLegend(0.17,0.72,0.42,0.88);
  leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetTextFont(42);
  leg->SetFillStyle(0); leg->SetMargin(0.25); //separation symbol-text
  leg->SetEntrySeparation(0.15);
  leg->AddEntry(hUS,"OS", "p");
  leg->AddEntry(hBackground, Form("LS*%4.2f",sig->GetScaleFactor()), "p");
  leg->Draw();

  c->cd(2);
  gPad->SetRightMargin(0.01);
  gPad->SetTopMargin(0);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(.15);
  /*
  gPad->SetGridx();
  gPad->SetGridy();
  */
  gPad->SetTickx();
  gPad->SetTicky();
  
  TH1* hSignal=sig->GetSignalHistogram();
  hSignal->GetXaxis()->SetTitle("M_{ee} (GeV/c^{2})");
  hSignal->GetXaxis()->CenterTitle();
  hSignal->GetXaxis()->SetTitleSize(0.06); hSignal->GetXaxis()->SetLabelSize(0.05);
  hSignal->GetXaxis()->SetTitleOffset(1.1);
  hSignal->GetXaxis()->SetNdivisions(510);
  hSignal->GetYaxis()->CenterTitle();
  hSignal->GetYaxis()->SetTitleSize(0.06); hSignal->GetYaxis()->SetLabelSize(0.05);
  hSignal->GetYaxis()->SetTitleOffset(1.1);
  hSignal->GetXaxis()->SetLabelFont(42); hSignal->GetYaxis()->SetLabelFont(42);
  hSignal->GetXaxis()->SetTitleFont(42); hSignal->GetYaxis()->SetTitleFont(42);
  hSignal->GetYaxis()->SetNdivisions(510);
  //hSignal->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
  //hSignal->SetLineWidth(.7);  
  hSignal->SetMarkerSize(.8);
  hSignal->SetMarkerColor(2); hSignal->SetMarkerStyle(20); hSignal->SetLineColor(2);

  hSignal->Draw();
  
  TLegend *leg2=new TLegend(0.16,0.8,0.57,0.93);
  leg2->SetBorderSize(0); leg2->SetFillColor(0); leg2->SetTextFont(42);
  leg2->SetFillStyle(0); leg2->SetMargin(0.25); //separation symbol-text
  leg2->SetEntrySeparation(0.15);
  leg2->AddEntry(hSignal, Form("OS-%4.2f*LS",sig->GetScaleFactor()), "p");
  TObject *hMmc=hSignal->GetListOfFunctions()->FindObject("mcMinv");
  if (hMmc) leg2->AddEntry(hMmc, Form("MC (#chi^{2}/dof=%3.1f)",TString(hMmc->GetTitle()).Atof()), "l");
  leg2->Draw();
  
  Double_t beforePhys=0;
  Double_t afterPhys=0;
  Double_t afterV0and=0;
  Double_t afterEventSel=0;
  Double_t afterPileupRej=0;

  if (hEventStat){
    Int_t beforePhysBin=hEventStat->GetXaxis()->FindBin("Before Phys. Sel.");
    Int_t afterPhysBin=hEventStat->GetXaxis()->FindBin("After Phys. Sel.");
    Int_t afterV0andBin=hEventStat->GetXaxis()->FindBin("V0and triggers");
    Int_t afterEventSelBin=hEventStat->GetXaxis()->FindBin("After Event Filter");
    Int_t afterPileupRejBin=hEventStat->GetXaxis()->FindBin("After Pileup rejection");
    
    beforePhys=hEventStat->GetBinContent(beforePhysBin);
    afterPhys=hEventStat->GetBinContent(afterPhysBin);
    afterV0and=hEventStat->GetBinContent(afterV0andBin);
    afterEventSel=hEventStat->GetBinContent(afterEventSelBin);
    afterPileupRej=hEventStat->GetBinContent(afterPileupRejBin);
    
    printf("Mevents: all: %4.1f, PhysSel: %4.1f, V0AND: %4.1f, EventSel: %4.1f, PileupRej: %4.1f\n",(Float_t)beforePhys/1.e+6,(Float_t)afterPhys/1.e+6,(Float_t)afterV0and/1.e+6,(Float_t)afterEventSel/1.e+6,(Float_t)afterPileupRej/1.e+6);
    
    lat->DrawLatex(0.18, 0.2, Form("%4.1f Mevents",(Float_t)afterPhys/1.e+6));
  }
  
  if (save){
//     c->SaveAs(Form("%s.eps",cname));
    c->SaveAs(Form("%s.png",cname));
/*
    FILE *out_file;
    if ( (out_file = fopen(Form("sig_%s.txt",cname), "w")) == NULL )
    {   fprintf(stderr, "Cannot open file %s\n", Form("sig_%s.txt",cname)); }
    fprintf(stdout, "Signal file: %s\n", Form("sig_%s.txt",cname));
    fprintf(out_file,"%3d %4.1f  %3.1f %4.2f  %4.1f %4.2f %d\n",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,(Int_t)afterPhys);
    fclose(out_file);

    TFile outMinv(Form("Minv_%s.root",cname), "RECREATE");
    hUS->Write();
    hBackground->Write();
    hSignal->Write();
    if (hMmc) hMmc->Write();
    outMinv.Close();*/

  }
  
  return;
}

//_______________________________________
void FitMCshape(AliDielectronSignalBase *sig)
{
  TFile mcFile(mcLineShapeFile);
  if (!mcFile.IsOpen()) {
    printf("mcMinv_LHC10e2 not found!!!\n");
    return;
  }
  TH1D *hMmc = (TH1D*)mcFile.Get("mcMinv");
  if (!hMmc) {
    printf("mcMinv not found!!\n");
    return;
  }
  hMmc->SetDirectory(0);
  hMmc->SetName("mcMinv");
  mcFile.Close();
  
  TH1* hMsub=sig->GetSignalHistogram();
  Double_t mb1=sig->GetIntegralMin();
  Double_t mb2=sig->GetIntegralMax();
  
  Double_t effInt = 0.;
  for(Int_t iBin=hMmc->FindBin(mb1); iBin<hMmc->FindBin(mb2); iBin++) {
    effInt += hMmc->GetBinContent(iBin);
  }
  effInt/=hMmc->Integral();
  printf("MC signal fraction in range %4.2f-%4.2f GeV: %5.3f \n",hMmc->GetBinLowEdge(hMmc->FindBin(mb1)),hMmc->GetBinLowEdge(hMmc->FindBin(mb2)+1),effInt);
 
  Float_t mcScale1=(hMsub->GetXaxis()->GetBinWidth(1)/hMmc->GetXaxis()->GetBinWidth(1))*
    hMsub->Integral(hMsub->FindBin(mb1),hMsub->FindBin(mb2))/
    hMmc->Integral(hMmc->FindBin(mb1),hMmc->FindBin(mb2));
  
  printf("1st guess of MC scale factor: %6.3f\n",mcScale1);
  
  Float_t mcScale=0.;
  Float_t chi2_min=100000.;
  Int_t iMin=0;
  Int_t ndf=0;
  
  for(Int_t i=0; i<20; i++){
    Float_t chi2=0.;
    Float_t scale=(0.4+0.05*(Float_t)i)*mcScale1;
    ndf=0;
    for(Int_t ib=1; ib<=hMsub->GetXaxis()->GetNbins(); ib++){
      Float_t data=(Float_t)hMsub->GetBinContent(ib);
      Float_t err=(Float_t)hMsub->GetBinError(ib);
      Float_t mc=scale*((Float_t)hMmc->GetBinContent(hMmc->FindBin(hMsub->GetBinCenter(ib))));
      if (err>0) {
        chi2 += ((data-mc)*(data-mc))/(err*err);
        ndf++;
      } else {
        printf("bin %d Err: %6.3f, chi2: %6.1f\n",ib,err,chi2);
      }
    }
      //printf("%d scale factor: %6.3f, chi2: %6.1f\n",i,scale,chi2);
    if(chi2 < chi2_min){
      chi2_min = chi2;
      mcScale = scale;
      iMin=i;
    }
  }
  //Float_t chi2dof=chi2_min/(Float_t)(hMinv->GetXaxis()->GetNbins()-1);
  Float_t chi2dof=chi2_min/((Float_t)(ndf-1));
  printf("MC fit (i=%d): chi2/dof: %6.3f/%d, Scale: %7.4f \n",iMin,chi2_min,(ndf-1),mcScale);
  hMmc->SetTitle(Form("%f",chi2dof));
  
  //mcScale=IntData/IntMC;printf("Int Data, MC: %10.1f %10.1f, MC scale: %6.3f\n",IntData,IntMC,mcScale);
  
  hMmc->Scale(mcScale);
  hMmc->SetOption("sameHISTC");
  hMsub->GetListOfFunctions()->Add(hMmc);
}
Commit	Line	Data
2a14a7b1	1	#include <TROOT.h>
	2	#include <TStyle.h>
	3	#include <TH1D.h>
	4	#include <TObjArray.h>
	5	#include <TCanvas.h>
	6	#include <TMath.h>
	7	#include <TAxis.h>
554e40f8	8	#include <TLatex.h>
	9	#include <TFile.h>
	10	#include <TString.h>
	11	#include <TLegend.h>
	12	#include <TList.h>
2a14a7b1	13
	14	#include "AliDielectronSignalExt.h"
	15	#include "AliDielectronCFdraw.h"
	16	#include "AliDielectron.h"
	17
554e40f8	18	AliDielectronSignalBase* GetSignalLS(AliDielectronCFdraw &d, Int_t step);
	19	AliDielectronSignalBase* GetSignalRot(AliDielectronCFdraw &d, Int_t step);
	20	void SetStyle(AliDielectronSignalBase sig, const char nameAdd);
	21	void DrawSpectra(AliDielectronSignalBase sig, const char cname, TH1 *hEventStat=0x0, Bool_t save=kFALSE);
	22	void FitMCshape(AliDielectronSignalBase *sig);
2a14a7b1	23
554e40f8	24	const char *mcLineShapeFile="$ALICE_ROOT/PWG3/dielectron/macros/mcMinv_LHC10e2.root";
2a14a7b1	25
2a14a7b1	26	//_______________________________________
554e40f8	27	void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t save=kFALSE)
2a14a7b1	28	{
554e40f8	29	AliDielectronCFdraw d(filenameData);
	30	AliDielectronCFdraw dCorr("corrCont","corrCont");
	31	TString nameCorr(filenameMC);
	32	if (!nameCorr.IsNull()) d.SetCFContainers(nameCorr.Data());
	33	TFile f(filenameData);
	34	TH1 hStats=(TH1)f.Get("hEventStat");
	35	hStats->SetDirectory(0);
	36	f.Close();
2a14a7b1	37
554e40f8	38	Int_t stepFirst=0, stepAny=1, stepTOFmix=2;
2a14a7b1	39
	40	gStyle->SetOptStat(0);
	41	//Set common Ranges
554e40f8	42	// d.SetRangeUser("Leg1_Pt",0.8,1000.);
554e40f8	43	// d.SetRangeUser("Leg2_Pt",0.8,1000.);
ffbede40	44	d.SetRangeUser("Leg1_NclsTPC",140.,170.);
	45	d.SetRangeUser("Leg2_NclsTPC",140.,170.);
	46	d.SetRangeUser("Leg1_Pt",1.01,100000);
	47	d.SetRangeUser("Leg2_Pt",1.01,100000);
	48	// d.SetRangeUser("Leg1_TPC_nSigma_Electrons",-3,3);
	49	// d.SetRangeUser("Leg2_TPC_nSigma_Electrons",-3,3);
	50	// d.SetRangeUser("Leg1_TPC_nSigma_Pions",3,20);
	51	// d.SetRangeUser("Leg2_TPC_nSigma_Pions",3,20);
	52	// d.SetRangeUser("Leg1_TPC_nSigma_Protons",3,20);
	53	// d.SetRangeUser("Leg2_TPC_nSigma_Protons",3,20);
	54
	55	d.SetRangeUser("M",0.5,5.);
2a14a7b1	56	//============================
	57	//SPD first
	58	//
554e40f8	59
2a14a7b1	60	//--- Like sign subtraction
554e40f8	61	AliDielectronSignalBase *sigFirst=GetSignalLS(d,stepFirst);
554e40f8	62	SetStyle(sigFirst,"ITS First - Like Sign subtraction");
ffbede40	63	// DrawSpectra(sigFirst,"cFirst",hStats,save);
2a14a7b1	64	//--- Rotation subtraction
554e40f8	65	AliDielectronSignalBase *sigFirstRot=GetSignalRot(d,stepFirst);
	66	SetStyle(sigFirstRot,"ITS First - Track rotation subtraction");
	67	// DrawSpectra(sigFirstRot,"cFirstRot",hStats,save);
2a14a7b1	68
	69	//============================
	70	//SPD any
	71	//
554e40f8	72	AliDielectronSignalBase *sigAny=GetSignalLS(d,stepAny);
	73	SetStyle(sigAny,"ITS First - Like Sign subtraction");
	74	DrawSpectra(sigAny,"cAny",hStats,save);
	75	//--- Rotation subtraction
	76	AliDielectronSignalBase *sigAnyRot=GetSignalRot(d,stepAny);
	77	SetStyle(sigAnyRot,"ITS First - Track rotation subtraction");
ffbede40	78	// DrawSpectra(sigAnyRot,"cAnyRot",hStats,save);
554e40f8	79
	80	//=============================
	81	//TOF up to 1.2, parametrisation in TPC ele
	82	//
	83	AliDielectronSignalBase *sigTOFmix=GetSignalLS(d,stepTOFmix);
	84	SetStyle(sigTOFmix,"TOF + TPC - Like Sign subtraction");
	85	// DrawSpectra(sigTOFmix,"cTOFTPC",hStats,save);
2a14a7b1	86	//--- Rotation subtraction
554e40f8	87	AliDielectronSignalBase *sigTOFmixRot=GetSignalRot(d,stepTOFmix);
554e40f8	88	SetStyle(sigTOFmixRot,"TOF + TPC - Track rotation subtraction");
ffbede40	89	// DrawSpectra(sigTOFmixRot,"cTOFTPCrot",hStats,save);
2a14a7b1	90
554e40f8	91	if (hStats) delete hStats;
2a14a7b1	92	}
	93
	94
	95	//_______________________________________
554e40f8	96	AliDielectronSignalBase *GetSignalLS(AliDielectronCFdraw &d, Int_t step)
2a14a7b1	97	{
	98	//
	99	// Get Extracted signal from likesign method
	100	//
	101
	102	TObjArray *arr=new TObjArray;
	103	arr->SetOwner();
554e40f8	104
2a14a7b1	105	for (Int_t iType=0;iType<3;++iType){
	106	d.SetRangeUser("PairType",iType,iType);
	107	arr->AddAt(d.Project("M",step),iType);
	108	}
554e40f8	109
2a14a7b1	110	AliDielectronSignalExt *sig=new AliDielectronSignalExt;
554e40f8	111	sig->SetScaleRawToBackground(3.2,5.);
554e40f8	112	sig->SetIntegralRange(2.92,3.15);
2a14a7b1	113	sig->SetMethod(AliDielectronSignalBase::kLikeSign);
	114	sig->Process(arr);
	115
	116	delete arr;
	117	return sig;
	118	}
	119
	120	//_______________________________________
554e40f8	121	AliDielectronSignalBase *GetSignalRot(AliDielectronCFdraw &d, Int_t step)
2a14a7b1	122	{
	123	//
	124	// Get Extracted signal from likesign method
	125	//
	126
	127	TObjArray *arr=new TObjArray;
	128	arr->SetOwner();
554e40f8	129
2a14a7b1	130	Int_t iType=AliDielectron::kEv1PM;
	131	d.SetRangeUser("PairType",iType,iType);
	132	arr->AddAt(d.Project("M",step),iType);
	133
	134	iType=AliDielectron::kEv1PMRot;
	135	d.SetRangeUser("PairType",iType,iType);
	136	arr->AddAt(d.Project("M",step),iType);
	137
	138	AliDielectronSignalExt *sig=new AliDielectronSignalExt;
ffbede40	139	// sig->SetScaleRawToBackground(3.2,4.);
554e40f8	140	sig->SetIntegralRange(2.93,3.15);
2a14a7b1	141	sig->SetMethod(AliDielectronSignalBase::kRotation);
	142	sig->Process(arr);
	143
554e40f8	144
2a14a7b1	145	delete arr;
	146	return sig;
	147	}
	148
	149	//_______________________________________
554e40f8	150	void SetStyle(AliDielectronSignalBase sig, const char nameAdd)
	151	{
	152	//
	153	//
	154	//
	155	TH1 *hUS=sig->GetUnlikeSignHistogram();
	156	hUS->SetMarkerStyle(20);
	157	hUS->SetMarkerSize(0.7);
	158	hUS->SetMarkerColor(kRed);
	159	hUS->SetLineColor(kRed);
	160	hUS->SetStats(0);
	161	hUS->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
	162	nameAdd,1000*hUS->GetXaxis()->GetBinWidth(1)));
	163
	164	TH1* hBackground=sig->GetBackgroundHistogram();
	165	hBackground->SetMarkerStyle(24);
	166	hBackground->SetMarkerSize(0.7);
	167	hBackground->SetStats(0);
	168	hBackground->SetMarkerColor(kBlue);
	169	hBackground->SetLineColor(kBlue);
	170	hBackground->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
	171	nameAdd,1000*hBackground->GetXaxis()->GetBinWidth(1)));
	172
	173	TH1* hSignal=sig->GetSignalHistogram();
	174	hSignal->SetMarkerStyle(20);
	175	hSignal->SetMarkerSize(0.7);
	176	hSignal->SetMarkerColor(kRed);
	177	hSignal->SetLineColor(kRed);
	178	hSignal->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
	179	nameAdd,1000*hSignal->GetXaxis()->GetBinWidth(1)));
	180
	181
	182	}
	183
	184	//_______________________________________
	185	void DrawSpectra(AliDielectronSignalBase sig, const char cname, TH1 *hEventStat, Bool_t save)
2a14a7b1	186	{
	187	//
	188	//
	189	//
554e40f8	190
	191
	192	FitMCshape(sig);
	193
2a14a7b1	194	gStyle->SetOptTitle(0);
2a14a7b1	195	TCanvas c=(TCanvas)gROOT->FindObject(cname);
554e40f8	196	if (!c) c=new TCanvas(cname,cname,400,500);
	197	c->SetTopMargin(0.04); c->SetRightMargin(0.04);
	198	c->SetLeftMargin(0.13); c->SetBottomMargin(0.14);
2a14a7b1	199	c->Clear();
	200	c->Divide(1,2,0,0);
	201	c->cd(1);
	202	gPad->SetTopMargin(0.01);
	203	gPad->SetRightMargin(0.01);
554e40f8	204	gPad->SetLeftMargin(0.15);
	205	gPad->SetBottomMargin(0.001);
	206	/*
2a14a7b1	207	gPad->SetGridx();
2a14a7b1	208	gPad->SetGridy();
554e40f8	209	*/
2a14a7b1	210	gPad->SetTickx();
2a14a7b1	211	gPad->SetTicky();
554e40f8	212
2a14a7b1	213	TH1 *hUS=sig->GetUnlikeSignHistogram();
2a14a7b1	214
2a14a7b1	215	TH1* hBackground=sig->GetBackgroundHistogram();
2a14a7b1	216
554e40f8	217	hUS->GetXaxis()->CenterTitle();
	218	hUS->GetXaxis()->SetTitleSize(0.07); hUS->GetXaxis()->SetLabelSize(0.059);
	219	// hUS->GetXaxis()->SetLabelOffset(1.);
	220	hUS->GetXaxis()->SetTitleOffset(1.1);
	221	hUS->GetXaxis()->SetNdivisions(510);
	222	hUS->GetYaxis()->CenterTitle();
	223	hUS->GetYaxis()->SetTitleSize(0.07); hUS->GetYaxis()->SetLabelSize(0.059);
	224	hUS->GetYaxis()->SetTitleOffset(1.);
	225	hUS->GetXaxis()->SetLabelFont(42); hUS->GetYaxis()->SetLabelFont(42);
	226	hUS->GetXaxis()->SetTitleFont(42); hUS->GetYaxis()->SetTitleFont(42);
	227	hUS->GetYaxis()->SetNdivisions(510);
	228	//hUS->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
	229	//hUS->SetLineWidth(.7);
	230	hUS->SetMarkerSize(.8);
	231	hUS->SetMarkerColor(2); hUS->SetMarkerStyle(20); hUS->SetLineColor(2);
2a14a7b1	232	hUS->Draw();
554e40f8	233	hUS->SetMaximum(hUS->GetMaximum()*1.3);
	234
	235	//hBackground->SetLineWidth(.4);
	236	hBackground->SetMarkerSize(.7);
	237	hBackground->SetMarkerColor(4); hBackground->SetMarkerStyle(27); hBackground->SetLineColor(4);
	238
2a14a7b1	239	hBackground->Draw("same");
2a14a7b1	240
554e40f8	241	TLatex *lat=new TLatex;
	242	lat->SetNDC(kTRUE);
	243	// lat->DrawLatex(0.68, 0.67, "ALICE Performance");
	244	lat->SetTextColor(1);lat->SetTextFont(42);lat->SetTextSize(.045);
	245
	246	Double_t sigN=sig->GetSignal();
	247	Double_t sigEr=sig->GetSignalError();
	248	Double_t sigS2B=sig->GetSB();
	249	Double_t sigS2Ber=sig->GetSBError();
	250	Double_t sigSignif= sig->GetSignificance();
	251	Double_t sigSignifEr= sig->GetSignificanceError();
	252	lat->DrawLatex(0.18, 0.92, Form("S: %3d#pm%4.1f, S/B: %3.1f#pm %4.2f, Signif.: %4.1f#pm%4.2f (%4.2f-%4.2f GeV) ",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,
	253	hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMin())),
	254	hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMax())+1)));
	255
	256	TLegend *leg=new TLegend(0.17,0.72,0.42,0.88);
	257	leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetTextFont(42);
	258	leg->SetFillStyle(0); leg->SetMargin(0.25); //separation symbol-text
	259	leg->SetEntrySeparation(0.15);
	260	leg->AddEntry(hUS,"OS", "p");
	261	leg->AddEntry(hBackground, Form("LS*%4.2f",sig->GetScaleFactor()), "p");
	262	leg->Draw();
	263
2a14a7b1	264	c->cd(2);
2a14a7b1	265	gPad->SetRightMargin(0.01);
554e40f8	266	gPad->SetTopMargin(0);
	267	gPad->SetLeftMargin(0.15);
	268	gPad->SetBottomMargin(.15);
	269	/*
2a14a7b1	270	gPad->SetGridx();
2a14a7b1	271	gPad->SetGridy();
554e40f8	272	*/
2a14a7b1	273	gPad->SetTickx();
2a14a7b1	274	gPad->SetTicky();
554e40f8	275
2a14a7b1	276	TH1* hSignal=sig->GetSignalHistogram();
554e40f8	277	hSignal->GetXaxis()->SetTitle("M_{ee} (GeV/c^{2})");
	278	hSignal->GetXaxis()->CenterTitle();
	279	hSignal->GetXaxis()->SetTitleSize(0.06); hSignal->GetXaxis()->SetLabelSize(0.05);
	280	hSignal->GetXaxis()->SetTitleOffset(1.1);
	281	hSignal->GetXaxis()->SetNdivisions(510);
	282	hSignal->GetYaxis()->CenterTitle();
	283	hSignal->GetYaxis()->SetTitleSize(0.06); hSignal->GetYaxis()->SetLabelSize(0.05);
	284	hSignal->GetYaxis()->SetTitleOffset(1.1);
	285	hSignal->GetXaxis()->SetLabelFont(42); hSignal->GetYaxis()->SetLabelFont(42);
	286	hSignal->GetXaxis()->SetTitleFont(42); hSignal->GetYaxis()->SetTitleFont(42);
	287	hSignal->GetYaxis()->SetNdivisions(510);
	288	//hSignal->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
	289	//hSignal->SetLineWidth(.7);
	290	hSignal->SetMarkerSize(.8);
	291	hSignal->SetMarkerColor(2); hSignal->SetMarkerStyle(20); hSignal->SetLineColor(2);
	292
2a14a7b1	293	hSignal->Draw();
554e40f8	294
	295	TLegend *leg2=new TLegend(0.16,0.8,0.57,0.93);
	296	leg2->SetBorderSize(0); leg2->SetFillColor(0); leg2->SetTextFont(42);
	297	leg2->SetFillStyle(0); leg2->SetMargin(0.25); //separation symbol-text
	298	leg2->SetEntrySeparation(0.15);
	299	leg2->AddEntry(hSignal, Form("OS-%4.2f*LS",sig->GetScaleFactor()), "p");
	300	TObject *hMmc=hSignal->GetListOfFunctions()->FindObject("mcMinv");
	301	if (hMmc) leg2->AddEntry(hMmc, Form("MC (#chi^{2}/dof=%3.1f)",TString(hMmc->GetTitle()).Atof()), "l");
	302	leg2->Draw();
	303
	304	Double_t beforePhys=0;
	305	Double_t afterPhys=0;
	306	Double_t afterV0and=0;
	307	Double_t afterEventSel=0;
	308	Double_t afterPileupRej=0;
	309
	310	if (hEventStat){
	311	Int_t beforePhysBin=hEventStat->GetXaxis()->FindBin("Before Phys. Sel.");
	312	Int_t afterPhysBin=hEventStat->GetXaxis()->FindBin("After Phys. Sel.");
	313	Int_t afterV0andBin=hEventStat->GetXaxis()->FindBin("V0and triggers");
	314	Int_t afterEventSelBin=hEventStat->GetXaxis()->FindBin("After Event Filter");
	315	Int_t afterPileupRejBin=hEventStat->GetXaxis()->FindBin("After Pileup rejection");
	316
	317	beforePhys=hEventStat->GetBinContent(beforePhysBin);
	318	afterPhys=hEventStat->GetBinContent(afterPhysBin);
	319	afterV0and=hEventStat->GetBinContent(afterV0andBin);
	320	afterEventSel=hEventStat->GetBinContent(afterEventSelBin);
	321	afterPileupRej=hEventStat->GetBinContent(afterPileupRejBin);
	322
	323	printf("Mevents: all: %4.1f, PhysSel: %4.1f, V0AND: %4.1f, EventSel: %4.1f, PileupRej: %4.1f\n",(Float_t)beforePhys/1.e+6,(Float_t)afterPhys/1.e+6,(Float_t)afterV0and/1.e+6,(Float_t)afterEventSel/1.e+6,(Float_t)afterPileupRej/1.e+6);
	324
	325	lat->DrawLatex(0.18, 0.2, Form("%4.1f Mevents",(Float_t)afterPhys/1.e+6));
	326	}
	327
	328	if (save){
ffbede40	329	// c->SaveAs(Form("%s.eps",cname));
554e40f8	330	c->SaveAs(Form("%s.png",cname));
ffbede40	331	/*
554e40f8	332	FILE *out_file;
	333	if ( (out_file = fopen(Form("sig_%s.txt",cname), "w")) == NULL )
	334	{ fprintf(stderr, "Cannot open file %s\n", Form("sig_%s.txt",cname)); }
	335	fprintf(stdout, "Signal file: %s\n", Form("sig_%s.txt",cname));
	336	fprintf(out_file,"%3d %4.1f %3.1f %4.2f %4.1f %4.2f %d\n",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,(Int_t)afterPhys);
	337	fclose(out_file);
ffbede40	338
554e40f8	339	TFile outMinv(Form("Minv_%s.root",cname), "RECREATE");
	340	hUS->Write();
	341	hBackground->Write();
	342	hSignal->Write();
ffbede40	343	if (hMmc) hMmc->Write();
ffbede40	344	outMinv.Close();*/
554e40f8	345
	346	}
	347
	348	return;
2a14a7b1	349	}
2a14a7b1	350
554e40f8	351	//_______________________________________
	352	void FitMCshape(AliDielectronSignalBase *sig)
	353	{
	354	TFile mcFile(mcLineShapeFile);
	355	if (!mcFile.IsOpen()) {
	356	printf("mcMinv_LHC10e2 not found!!!\n");
	357	return;
	358	}
	359	TH1D hMmc = (TH1D)mcFile.Get("mcMinv");
	360	if (!hMmc) {
	361	printf("mcMinv not found!!\n");
	362	return;
	363	}
	364	hMmc->SetDirectory(0);
	365	hMmc->SetName("mcMinv");
	366	mcFile.Close();
	367
	368	TH1* hMsub=sig->GetSignalHistogram();
	369	Double_t mb1=sig->GetIntegralMin();
	370	Double_t mb2=sig->GetIntegralMax();
	371
	372	Double_t effInt = 0.;
	373	for(Int_t iBin=hMmc->FindBin(mb1); iBin<hMmc->FindBin(mb2); iBin++) {
	374	effInt += hMmc->GetBinContent(iBin);
	375	}
	376	effInt/=hMmc->Integral();
	377	printf("MC signal fraction in range %4.2f-%4.2f GeV: %5.3f \n",hMmc->GetBinLowEdge(hMmc->FindBin(mb1)),hMmc->GetBinLowEdge(hMmc->FindBin(mb2)+1),effInt);
	378
	379	Float_t mcScale1=(hMsub->GetXaxis()->GetBinWidth(1)/hMmc->GetXaxis()->GetBinWidth(1))*
	380	hMsub->Integral(hMsub->FindBin(mb1),hMsub->FindBin(mb2))/
	381	hMmc->Integral(hMmc->FindBin(mb1),hMmc->FindBin(mb2));
	382
	383	printf("1st guess of MC scale factor: %6.3f\n",mcScale1);
	384
	385	Float_t mcScale=0.;
	386	Float_t chi2_min=100000.;
	387	Int_t iMin=0;
	388	Int_t ndf=0;
	389
	390	for(Int_t i=0; i<20; i++){
	391	Float_t chi2=0.;
	392	Float_t scale=(0.4+0.05(Float_t)i)mcScale1;
	393	ndf=0;
	394	for(Int_t ib=1; ib<=hMsub->GetXaxis()->GetNbins(); ib++){
	395	Float_t data=(Float_t)hMsub->GetBinContent(ib);
	396	Float_t err=(Float_t)hMsub->GetBinError(ib);
	397	Float_t mc=scale*((Float_t)hMmc->GetBinContent(hMmc->FindBin(hMsub->GetBinCenter(ib))));
	398	if (err>0) {
	399	chi2 += ((data-mc)(data-mc))/(errerr);
	400	ndf++;
	401	} else {
	402	printf("bin %d Err: %6.3f, chi2: %6.1f\n",ib,err,chi2);
	403	}
	404	}
	405	//printf("%d scale factor: %6.3f, chi2: %6.1f\n",i,scale,chi2);
	406	if(chi2 < chi2_min){
	407	chi2_min = chi2;
	408	mcScale = scale;
	409	iMin=i;
	410	}
	411	}
	412	//Float_t chi2dof=chi2_min/(Float_t)(hMinv->GetXaxis()->GetNbins()-1);
	413	Float_t chi2dof=chi2_min/((Float_t)(ndf-1));
	414	printf("MC fit (i=%d): chi2/dof: %6.3f/%d, Scale: %7.4f \n",iMin,chi2_min,(ndf-1),mcScale);
415	hMmc->SetTitle(Form("%f",chi2dof));
416
417	//mcScale=IntData/IntMC;printf("Int Data, MC: %10.1f %10.1f, MC scale: %6.3f\n",IntData,IntMC,mcScale);
418
419	hMmc->Scale(mcScale);
420	hMmc->SetOption("sameHISTC");
421	hMsub->GetListOfFunctions()->Add(hMmc);
422	}