[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronSignalExt.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                                                                       //
//                      Dielectron SignalExt                             //
//                                                                       //
/*
Dielectron signal extraction class using functions as input.

A function to describe the signal as well as one to describe the background
has to be deployed by the user. Alternatively on of the default implementaions
can be used.

*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////
#include <TF1.h>
#include <TH1.h>
#include <TH2F.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TString.h>
#include <TLine.h>

#include <AliLog.h>

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

ClassImp(AliDielectronSignalExt)

AliDielectronSignalExt::AliDielectronSignalExt() :
  AliDielectronSignalBase()
{
  //
  // Default Constructor
  //
}

//______________________________________________
AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
  AliDielectronSignalBase(name, title)
{
  //
  // Named Constructor
  //
}

//______________________________________________
AliDielectronSignalExt::~AliDielectronSignalExt()
{
  //
  // Default Destructor
  //
}

//______________________________________________
void AliDielectronSignalExt::Process(TObjArray* const arrhist)
{
  // 
  // signal subtraction. support like-sign subtraction and event mixing method
  //
  switch ( fMethod ){
    case kLikeSign :
    case kLikeSignArithm :
    case kLikeSignRcorr:
    case kLikeSignArithmRcorr:
      ProcessLS(arrhist);    // process like-sign subtraction method
      break;

    case kEventMixing : 
      ProcessEM(arrhist);    // process event mixing method
      break;

  case kRotation:
      ProcessRotation(arrhist);
      break;

    default :
      AliWarning("Subtraction method not supported. Please check SetMethod() function.");
  }
}

//______________________________________________
void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
{
  //
  // signal subtraction 
  //
  fHistDataPP = (TH1*)(arrhist->At(AliDielectron::kEv1PP))->Clone("histPP");  // ++    SE
  fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
  fHistDataMM = (TH1*)(arrhist->At(AliDielectron::kEv1MM))->Clone("histMM");  // --    SE
  fHistDataPP->Sumw2();
  fHistDataPM->Sumw2();
  fHistDataMM->Sumw2();
  fHistDataPP->SetDirectory(0);
  fHistDataPM->SetDirectory(0);
  fHistDataMM->SetDirectory(0);
  
  // rebin the histograms
  if (fRebin>1) { 
    fHistDataPP->Rebin(fRebin);
    fHistDataPM->Rebin(fRebin);
    fHistDataMM->Rebin(fRebin);
  }       

  fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
                          fHistDataPM->GetXaxis()->GetNbins(),
                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  fHistRfactor->Sumw2();
  fHistRfactor->SetDirectory(0);
  
  fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
			 fHistDataPM->GetXaxis()->GetNbins(),
			 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  fHistSignal->SetDirectory(0);
  fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
			     fHistDataPM->GetXaxis()->GetNbins(),
			     fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  fHistBackground->SetDirectory(0);
  
  // fill out background histogram
  for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
    Float_t pp = fHistDataPP->GetBinContent(ibin);
    Float_t mm = fHistDataMM->GetBinContent(ibin);

    Float_t background = 2*TMath::Sqrt(pp*mm);
    Float_t ebackground = TMath::Sqrt(mm+pp);
    if (fMethod==kLikeSignArithm || fMethod==kLikeSignArithmRcorr ){
      //Arithmetic mean instead of geometric
      background=(pp+mm);
      ebackground=TMath::Sqrt(pp+mm);
      if (TMath::Abs(ebackground)<1e-30) ebackground=1;
    }

    fHistBackground->SetBinContent(ibin, background);
    fHistBackground->SetBinError(ibin, ebackground);
  }

  //correct LS spectrum bin-by-bin with R factor obtained in mixed events
  if(fMixingCorr || fMethod==kLikeSignRcorr || fMethod==kLikeSignArithmRcorr) {
    
    TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP");  // ++    ME
    TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM");  // --    ME

    TH1* histMixPM = 0x0;
    if (arrhist->At(AliDielectron::kEv1MEv2P)){
      histMixPM   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM");  // -+    ME
    }

    if (arrhist->At(AliDielectron::kEv1PEv2M)){
      TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
      if (!histMixPM) fHistDataME   = (TH1*)htmp->Clone("mixPM");                   // +-    ME
      else histMixPM->Add(htmp);
    }

    histMixPM->Sumw2();
    
    // rebin the histograms
    if (fRebin>1) { 
      histMixPP->Rebin(fRebin);
      histMixMM->Rebin(fRebin);
      histMixPM->Rebin(fRebin);
    }       
    
    // fill out rfactor histogram
    for(Int_t ibin=1; ibin<=histMixPM->GetXaxis()->GetNbins(); ibin++) {
      Float_t pp  = histMixPP->GetBinContent(ibin);
      Float_t ppe = histMixPP->GetBinError(ibin);
      Float_t mm  = histMixMM->GetBinContent(ibin);
      Float_t mme = histMixMM->GetBinError(ibin);
      Float_t pm  = histMixPM->GetBinContent(ibin);
      Float_t pme = histMixPM->GetBinError(ibin);
      
      Float_t background = 2*TMath::Sqrt(pp*mm);
      // do not use it since ME could be weighted err!=sqrt(entries)
      //      Float_t ebackground = TMath::Sqrt(mm+pp); 
      Float_t ebackground = TMath::Sqrt(mm/pp*ppe*ppe + pp/mm*mme*mme);
      if (fMethod==kLikeSignArithm){
        //Arithmetic mean instead of geometric
        background=(pp+mm);
        ebackground=TMath::Sqrt(ppe*ppe+mme*mme);
        if (TMath::Abs(ebackground)<1e-30) ebackground=1;
      }
      
      Float_t rcon = 1.0;
      Float_t rerr = 0.0;
      if(background>0.0) {
        rcon = pm/background;
        rerr = TMath::Sqrt((1./background)*(1./background) * pme*pme +
                           (pm/(background*background))*(pm/(background*background)) * ebackground*ebackground);
      }
      fHistRfactor->SetBinContent(ibin, rcon);
      fHistRfactor->SetBinError(ibin, rerr);
    }
    
    fHistBackground->Multiply(fHistRfactor);
    
    if (histMixPP) delete histMixPP;
    if (histMixMM) delete histMixMM;
    if (histMixPM) delete histMixPM;
  }
  
  //scale histograms to match integral between fScaleMin and fScaleMax
  // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
  if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
  else if (fScaleMin>0.){
    fScaleFactor=fScaleMin;
    fHistBackground->Scale(fScaleFactor);
  }

  //subract background
  fHistSignal->Add(fHistDataPM);
  fHistSignal->Add(fHistBackground,-1);
  
  // signal
  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
	  			             fHistSignal->FindBin(fIntMax), fErrors(0));
  // background
  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
						 fHistBackground->FindBin(fIntMax), 
						 fErrors(1));
  //printf("%f  %f\n",fValues(0),fValues(1));
  // S/B and significance
  SetSignificanceAndSOB();

  fProcessed = kTRUE;
}

//______________________________________________
void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
{
  //
  // event mixing of +- and -+
  //

  if (!arrhist->At(AliDielectron::kEv1PM) || !(arrhist->At(AliDielectron::kEv1MEv2P) || arrhist->At(AliDielectron::kEv1PEv2M)) ){
    AliError("Either OS or mixed histogram missing");
    return;
  }

  delete fHistDataPM; fHistDataPM=0x0;
  delete fHistDataME; fHistDataME=0x0;
  delete fHistBackground; fHistBackground=0x0;
  
  fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
  fHistDataPM->Sumw2();
  fHistDataPM->SetDirectory(0x0);

  if (arrhist->At(AliDielectron::kEv1MEv2P)){
    fHistDataME   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME");  // -+    ME
  }
  
  if (arrhist->At(AliDielectron::kEv1PEv2M)){
    TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
    if (!fHistDataME) fHistDataME   = (TH1*)htmp->Clone("histMPME");  // -+    ME
    else fHistDataME->Add(htmp);
  }

  fHistBackground = (TH1*)fHistDataME->Clone("ME_Background");
  fHistBackground->SetDirectory(0x0);
  fHistBackground->Sumw2();

  // rebin the histograms
  if (fRebin>1) {
    fHistDataPM->Rebin(fRebin);
    fHistDataME->Rebin(fRebin);
    fHistBackground->Rebin(fRebin);
  }

  //scale histograms to match integral between fScaleMin and fScaleMax
  // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
  if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
  else if (fScaleMin>0.){
    fScaleFactor=fScaleMin;
    fHistBackground->Scale(fScaleFactor);
  }

  fHistSignal=(TH1*)fHistDataPM->Clone("Signal");
  fHistSignal->Add(fHistBackground,-1.);

    // signal
  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
                                             fHistSignal->FindBin(fIntMax), fErrors(0));
  // background
  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
                                                 fHistBackground->FindBin(fIntMax),
                                                 fErrors(1));
  // S/B and significance
  SetSignificanceAndSOB();

  fProcessed = kTRUE;
}

//______________________________________________
void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
{
  //
  // signal subtraction
  //
  fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
  if (!fHistDataPM){
    AliError("Unlike sign histogram not available. Cannot extract the signal.");
    return;
  }
  fHistDataPM->Sumw2();

  fHistBackground = (TH1*)(arrhist->At(AliDielectron::kEv1PMRot))->Clone("histRotation");
  if (!fHistBackground){
    AliError("Histgram from rotation not available. Cannot extract the signal.");
    delete fHistDataPM;
    fHistDataPM=0x0;
    return;
  }
  fHistBackground->Sumw2();

  // rebin the histograms
  if (fRebin>1) {
    fHistDataPM->Rebin(fRebin);
    fHistBackground->Rebin(fRebin);
  }

  //scale histograms to match integral between fScaleMin and fScaleMax
  // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
  if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
  else if (fScaleMin>0.){
    fScaleFactor=fScaleMin;
    fHistBackground->Scale(fScaleFactor);
  }

  fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
  fHistSignal->Add(fHistBackground,-1.);

    // signal
  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
                                             fHistSignal->FindBin(fIntMax), fErrors(0));
  // background
  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
                                                 fHistBackground->FindBin(fIntMax),
                                                 fErrors(1));
  // S/B and significance
  SetSignificanceAndSOB();
  
  fProcessed = kTRUE;
  
}

//______________________________________________
void AliDielectronSignalExt::Draw(const Option_t* option)
{
  //
  // Draw the fitted function
  //
  TString drawOpt(option); 
  drawOpt.ToLower();   

  Float_t minY = 0.001;
  Float_t maxY = 1.2*fHistDataPM->GetMaximum();
  Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
  Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
  Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
  Float_t minMinY = fHistSignal->GetMinimum();

  TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
  cSub->SetLeftMargin(0.15);
  cSub->SetRightMargin(0.0);
  cSub->SetTopMargin(0.002);
  cSub->SetBottomMargin(0.0);
  cSub->Divide(2,2,0.,0.);
  cSub->Draw();

  TVirtualPad* pad = cSub->cd(1);
  pad->SetLeftMargin(0.15);
  pad->SetRightMargin(0.0);
  pad->SetTopMargin(0.005);
  pad->SetBottomMargin(0.0);
  TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
  range1->SetStats(kFALSE);
  range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
  range1->GetYaxis()->CenterTitle();
  range1->GetYaxis()->SetLabelSize(0.05);
  range1->GetYaxis()->SetTitleSize(0.06);
  range1->GetYaxis()->SetTitleOffset(0.8);
  range1->Draw();
  fHistDataPM->SetLineColor(1);
  fHistDataPM->SetLineWidth(2);
  //  fHistDataPM->SetMarkerStyle(21);
  fHistDataPM->Draw("Psame");
  TLatex *latex = new TLatex();
  latex->SetNDC();
  latex->SetTextSize(0.05);
  latex->DrawLatex(0.2, 0.95, "Background un-substracted");
  TLine line;
  line.SetLineWidth(1);
  line.SetLineStyle(2);
  line.DrawLine(fIntMin, minY, fIntMin, maxY);
  line.DrawLine(fIntMax, minY, fIntMax, maxY);

  pad = cSub->cd(2);
  pad->SetLeftMargin(0.);
  pad->SetRightMargin(0.005);
  pad->SetTopMargin(0.005);
  pad->SetBottomMargin(0.0);
  TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
  range2->SetStats(kFALSE);
  range2->Draw();
  fHistBackground->SetLineColor(4);
  fHistBackground->SetLineWidth(2);
  //  fHistBackground->SetMarkerColor(4);
  //  fHistBackground->SetMarkerStyle(6);
  fHistBackground->Draw("Psame");
  latex->DrawLatex(0.05, 0.95, "Like-sign background");
  line.DrawLine(fIntMin, minY, fIntMin, maxY);
  line.DrawLine(fIntMax, minY, fIntMax, maxY);
  TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
  legend->SetFillColor(0);
  legend->SetMargin(0.15);
  legend->AddEntry(fHistDataPM, "N_{+-}", "l");
  legend->AddEntry(fHistDataPP, "N_{++}", "l");
  legend->AddEntry(fHistDataMM, "N_{--}", "l");
  legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
  legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
  legend->Draw();

  
  pad = cSub->cd(3);
  pad->SetLeftMargin(0.15);
  pad->SetRightMargin(0.0);
  pad->SetTopMargin(0.0);
  pad->SetBottomMargin(0.15);
  TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
  range3->SetStats(kFALSE);
  range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
  range3->GetYaxis()->CenterTitle();
  range3->GetYaxis()->SetLabelSize(0.05);
  range3->GetYaxis()->SetTitleSize(0.06);
  range3->GetYaxis()->SetTitleOffset(0.8);
  range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
  range3->GetXaxis()->CenterTitle();
  range3->GetXaxis()->SetLabelSize(0.05);
  range3->GetXaxis()->SetTitleSize(0.06);
  range3->GetXaxis()->SetTitleOffset(1.0);
  range3->Draw();
  fHistDataPM->Draw("Psame");
  fHistDataPP->SetLineWidth(2);
  fHistDataPP->SetLineColor(6);
  fHistDataMM->SetLineWidth(2);
  fHistDataMM->SetLineColor(8);
  fHistDataPP->Draw("Psame");
  fHistDataMM->Draw("Psame");
  line.DrawLine(minX, 0.,maxX, 0.);
  line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
  line.DrawLine(fIntMax, minMinY, fIntMax, maxY);

  pad = cSub->cd(4);
  pad->SetLeftMargin(0.0);
  pad->SetRightMargin(0.005);
  pad->SetTopMargin(0.0);
  pad->SetBottomMargin(0.15);
  TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
  range4->SetStats(kFALSE);
  range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
  range4->GetXaxis()->CenterTitle();
  range4->GetXaxis()->SetLabelSize(0.05);
  range4->GetXaxis()->SetTitleSize(0.06);
  range4->GetXaxis()->SetTitleOffset(1.0);
  range4->Draw();
  fHistSignal->SetLineWidth(2);
  fHistSignal->SetLineColor(2);
  fHistSignal->Draw("Psame");
  latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
  if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
  line.DrawLine(minX, 0.,maxX, 0.);
  line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
  line.DrawLine(fIntMax, minMinY, fIntMax, maxY);

  cSub->SaveAs(Form("%s_summary.png", fName.Data()));
}
Commit	Line	Data
572b0139	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	///////////////////////////////////////////////////////////////////////////
	17	// //
	18	// Dielectron SignalExt //
	19	// //
	20	/*
	21	Dielectron signal extraction class using functions as input.
	22
	23	A function to describe the signal as well as one to describe the background
	24	has to be deployed by the user. Alternatively on of the default implementaions
	25	can be used.
	26
	27	*/
	28	// //
	29	///////////////////////////////////////////////////////////////////////////
572b0139	30	#include <TF1.h>
572b0139	31	#include <TH1.h>
bc75eeb5	32	#include <TH2F.h>
	33	#include <TLatex.h>
	34	#include <TLegend.h>
572b0139	35	#include <TCanvas.h>
	36	#include <TMath.h>
	37	#include <TString.h>
bc75eeb5	38	#include <TLine.h>
572b0139	39
	40	#include <AliLog.h>
	41
	42	#include "AliDielectronSignalExt.h"
ac390e40	43	#include "AliDielectron.h"
572b0139	44
	45	ClassImp(AliDielectronSignalExt)
	46
	47	AliDielectronSignalExt::AliDielectronSignalExt() :
bc75eeb5	48	AliDielectronSignalBase()
572b0139	49	{
	50	//
	51	// Default Constructor
	52	//
572b0139	53	}
	54
	55	//______________________________________________
	56	AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
bc75eeb5	57	AliDielectronSignalBase(name, title)
572b0139	58	{
	59	//
	60	// Named Constructor
	61	//
	62	}
	63
	64	//______________________________________________
	65	AliDielectronSignalExt::~AliDielectronSignalExt()
	66	{
	67	//
	68	// Default Destructor
	69	//
572b0139	70	}
	71
	72	//______________________________________________
	73	void AliDielectronSignalExt::Process(TObjArray* const arrhist)
	74	{
	75	//
	76	// signal subtraction. support like-sign subtraction and event mixing method
	77	//
	78	switch ( fMethod ){
bc75eeb5	79	case kLikeSign :
45b2b1b8	80	case kLikeSignArithm :
b9d223bb	81	case kLikeSignRcorr:
b9d223bb	82	case kLikeSignArithmRcorr:
572b0139	83	ProcessLS(arrhist); // process like-sign subtraction method
	84	break;
	85
bc75eeb5	86	case kEventMixing :
572b0139	87	ProcessEM(arrhist); // process event mixing method
	88	break;
	89
2a14a7b1	90	case kRotation:
	91	ProcessRotation(arrhist);
	92	break;
	93
572b0139	94	default :
	95	AliWarning("Subtraction method not supported. Please check SetMethod() function.");
	96	}
	97	}
	98
	99	//______________________________________________
	100	void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
	101	{
	102	//
	103	// signal subtraction
	104	//
ac390e40	105	fHistDataPP = (TH1*)(arrhist->At(AliDielectron::kEv1PP))->Clone("histPP"); // ++ SE
	106	fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
	107	fHistDataMM = (TH1*)(arrhist->At(AliDielectron::kEv1MM))->Clone("histMM"); // -- SE
bc75eeb5	108	fHistDataPP->Sumw2();
	109	fHistDataPM->Sumw2();
	110	fHistDataMM->Sumw2();
554e40f8	111	fHistDataPP->SetDirectory(0);
	112	fHistDataPM->SetDirectory(0);
	113	fHistDataMM->SetDirectory(0);
572b0139	114
bc75eeb5	115	// rebin the histograms
	116	if (fRebin>1) {
	117	fHistDataPP->Rebin(fRebin);
	118	fHistDataPM->Rebin(fRebin);
	119	fHistDataMM->Rebin(fRebin);
	120	}
	121
ac390e40	122	fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
	123	fHistDataPM->GetXaxis()->GetNbins(),
	124	fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
	125	fHistRfactor->Sumw2();
	126	fHistRfactor->SetDirectory(0);
	127
2a14a7b1	128	fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
bc75eeb5	129	fHistDataPM->GetXaxis()->GetNbins(),
bc75eeb5	130	fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
554e40f8	131	fHistSignal->SetDirectory(0);
2a14a7b1	132	fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
bc75eeb5	133	fHistDataPM->GetXaxis()->GetNbins(),
bc75eeb5	134	fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
554e40f8	135	fHistBackground->SetDirectory(0);
554e40f8	136
ac390e40	137	// fill out background histogram
bc75eeb5	138	for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
bc75eeb5	139	Float_t pp = fHistDataPP->GetBinContent(ibin);
bc75eeb5	140	Float_t mm = fHistDataMM->GetBinContent(ibin);
bc75eeb5	141
	142	Float_t background = 2TMath::Sqrt(ppmm);
	143	Float_t ebackground = TMath::Sqrt(mm+pp);
b9d223bb	144	if (fMethod==kLikeSignArithm \|\| fMethod==kLikeSignArithmRcorr ){
45b2b1b8	145	//Arithmetic mean instead of geometric
	146	background=(pp+mm);
	147	ebackground=TMath::Sqrt(pp+mm);
	148	if (TMath::Abs(ebackground)<1e-30) ebackground=1;
	149	}
bc75eeb5	150
bc75eeb5	151	fHistBackground->SetBinContent(ibin, background);
bc75eeb5	152	fHistBackground->SetBinError(ibin, ebackground);
572b0139	153	}
ac390e40	154
ac390e40	155	//correct LS spectrum bin-by-bin with R factor obtained in mixed events
b9d223bb	156	if(fMixingCorr \|\| fMethod==kLikeSignRcorr \|\| fMethod==kLikeSignArithmRcorr) {
ac390e40	157
	158	TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP"); // ++ ME
	159	TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM"); // -- ME
b9d223bb	160
	161	TH1* histMixPM = 0x0;
	162	if (arrhist->At(AliDielectron::kEv1MEv2P)){
	163	histMixPM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM"); // -+ ME
	164	}
	165
	166	if (arrhist->At(AliDielectron::kEv1PEv2M)){
	167	TH1 htmp=(TH1)(arrhist->At(AliDielectron::kEv1PEv2M));
	168	if (!histMixPM) fHistDataME = (TH1*)htmp->Clone("mixPM"); // +- ME
	169	else histMixPM->Add(htmp);
	170	}
	171
ac390e40	172	histMixPM->Sumw2();
ac390e40	173
	174	// rebin the histograms
	175	if (fRebin>1) {
	176	histMixPP->Rebin(fRebin);
	177	histMixMM->Rebin(fRebin);
	178	histMixPM->Rebin(fRebin);
	179	}
	180
	181	// fill out rfactor histogram
	182	for(Int_t ibin=1; ibin<=histMixPM->GetXaxis()->GetNbins(); ibin++) {
	183	Float_t pp = histMixPP->GetBinContent(ibin);
	184	Float_t ppe = histMixPP->GetBinError(ibin);
	185	Float_t mm = histMixMM->GetBinContent(ibin);
	186	Float_t mme = histMixMM->GetBinError(ibin);
	187	Float_t pm = histMixPM->GetBinContent(ibin);
	188	Float_t pme = histMixPM->GetBinError(ibin);
	189
	190	Float_t background = 2TMath::Sqrt(ppmm);
	191	// do not use it since ME could be weighted err!=sqrt(entries)
	192	// Float_t ebackground = TMath::Sqrt(mm+pp);
	193	Float_t ebackground = TMath::Sqrt(mm/ppppeppe + pp/mmmmemme);
	194	if (fMethod==kLikeSignArithm){
	195	//Arithmetic mean instead of geometric
	196	background=(pp+mm);
	197	ebackground=TMath::Sqrt(ppeppe+mmemme);
	198	if (TMath::Abs(ebackground)<1e-30) ebackground=1;
	199	}
	200
	201	Float_t rcon = 1.0;
	202	Float_t rerr = 0.0;
	203	if(background>0.0) {
	204	rcon = pm/background;
	205	rerr = TMath::Sqrt((1./background)(1./background) pme*pme +
	206	(pm/(backgroundbackground))(pm/(backgroundbackground)) ebackground*ebackground);
	207	}
	208	fHistRfactor->SetBinContent(ibin, rcon);
	209	fHistRfactor->SetBinError(ibin, rerr);
	210	}
	211
	212	fHistBackground->Multiply(fHistRfactor);
	213
	214	if (histMixPP) delete histMixPP;
	215	if (histMixMM) delete histMixMM;
	216	if (histMixPM) delete histMixPM;
	217	}
	218
fb7d2d99	219	//scale histograms to match integral between fScaleMin and fScaleMax
5720c765	220	// or if fScaleMax < fScaleMin use fScaleMin as scale factor
554e40f8	221	if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
5720c765	222	else if (fScaleMin>0.){
	223	fScaleFactor=fScaleMin;
	224	fHistBackground->Scale(fScaleFactor);
	225	}
554e40f8	226
554e40f8	227	//subract background
ac390e40	228	fHistSignal->Add(fHistDataPM);
554e40f8	229	fHistSignal->Add(fHistBackground,-1);
554e40f8	230
bc75eeb5	231	// signal
	232	fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
	233	fHistSignal->FindBin(fIntMax), fErrors(0));
	234	// background
	235	fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
48609e3d	236	fHistBackground->FindBin(fIntMax),
48609e3d	237	fErrors(1));
ac390e40	238	//printf("%f %f\n",fValues(0),fValues(1));
bc75eeb5	239	// S/B and significance
572b0139	240	SetSignificanceAndSOB();
572b0139	241
bc75eeb5	242	fProcessed = kTRUE;
572b0139	243	}
	244
	245	//______________________________________________
	246	void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
	247	{
	248	//
ac390e40	249	// event mixing of +- and -+
572b0139	250	//
5720c765	251
b9d223bb	252	if (!arrhist->At(AliDielectron::kEv1PM) \|\| !(arrhist->At(AliDielectron::kEv1MEv2P) \|\| arrhist->At(AliDielectron::kEv1PEv2M)) ){
	253	AliError("Either OS or mixed histogram missing");
	254	return;
5720c765	255	}
5720c765	256
b9d223bb	257	delete fHistDataPM; fHistDataPM=0x0;
	258	delete fHistDataME; fHistDataME=0x0;
	259	delete fHistBackground; fHistBackground=0x0;
	260
	261	fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
	262	fHistDataPM->Sumw2();
	263	fHistDataPM->SetDirectory(0x0);
5720c765	264
b9d223bb	265	if (arrhist->At(AliDielectron::kEv1MEv2P)){
	266	fHistDataME = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME"); // -+ ME
	267	}
	268
	269	if (arrhist->At(AliDielectron::kEv1PEv2M)){
	270	TH1 htmp=(TH1)(arrhist->At(AliDielectron::kEv1PEv2M));
	271	if (!fHistDataME) fHistDataME = (TH1*)htmp->Clone("histMPME"); // -+ ME
	272	else fHistDataME->Add(htmp);
	273	}
5720c765	274
b9d223bb	275	fHistBackground = (TH1*)fHistDataME->Clone("ME_Background");
	276	fHistBackground->SetDirectory(0x0);
	277	fHistBackground->Sumw2();
	278
	279	// rebin the histograms
	280	if (fRebin>1) {
	281	fHistDataPM->Rebin(fRebin);
	282	fHistDataME->Rebin(fRebin);
	283	fHistBackground->Rebin(fRebin);
5720c765	284	}
5720c765	285
b9d223bb	286	//scale histograms to match integral between fScaleMin and fScaleMax
b9d223bb	287	// or if fScaleMax < fScaleMin use fScaleMin as scale factor
5720c765	288	if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
b9d223bb	289	else if (fScaleMin>0.){
5720c765	290	fScaleFactor=fScaleMin;
	291	fHistBackground->Scale(fScaleFactor);
	292	}
	293
b9d223bb	294	fHistSignal=(TH1*)fHistDataPM->Clone("Signal");
b9d223bb	295	fHistSignal->Add(fHistBackground,-1.);
5720c765	296
b9d223bb	297	// signal
5720c765	298	fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
b9d223bb	299	fHistSignal->FindBin(fIntMax), fErrors(0));
5720c765	300	// background
5720c765	301	fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
b9d223bb	302	fHistBackground->FindBin(fIntMax),
b9d223bb	303	fErrors(1));
5720c765	304	// S/B and significance
	305	SetSignificanceAndSOB();
	306
	307	fProcessed = kTRUE;
572b0139	308	}
572b0139	309
2a14a7b1	310	//______________________________________________
	311	void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
	312	{
	313	//
	314	// signal subtraction
	315	//
ac390e40	316	fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
2a14a7b1	317	if (!fHistDataPM){
	318	AliError("Unlike sign histogram not available. Cannot extract the signal.");
	319	return;
	320	}
	321	fHistDataPM->Sumw2();
	322
ac390e40	323	fHistBackground = (TH1*)(arrhist->At(AliDielectron::kEv1PMRot))->Clone("histRotation");
2a14a7b1	324	if (!fHistBackground){
	325	AliError("Histgram from rotation not available. Cannot extract the signal.");
	326	delete fHistDataPM;
	327	fHistDataPM=0x0;
	328	return;
	329	}
ba15fdfb	330	fHistBackground->Sumw2();
	331
	332	// rebin the histograms
	333	if (fRebin>1) {
	334	fHistDataPM->Rebin(fRebin);
	335	fHistBackground->Rebin(fRebin);
	336	}
2a14a7b1	337
554e40f8	338	//scale histograms to match integral between fScaleMin and fScaleMax
5720c765	339	// or if fScaleMax < fScaleMin use fScaleMin as scale factor
554e40f8	340	if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
5720c765	341	else if (fScaleMin>0.){
	342	fScaleFactor=fScaleMin;
	343	fHistBackground->Scale(fScaleFactor);
	344	}
	345
2a14a7b1	346	fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
	347	fHistSignal->Add(fHistBackground,-1.);
	348
	349	// signal
	350	fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
	351	fHistSignal->FindBin(fIntMax), fErrors(0));
	352	// background
	353	fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
	354	fHistBackground->FindBin(fIntMax),
	355	fErrors(1));
	356	// S/B and significance
	357	SetSignificanceAndSOB();
	358
	359	fProcessed = kTRUE;
	360
	361	}
	362
572b0139	363	//______________________________________________
	364	void AliDielectronSignalExt::Draw(const Option_t* option)
	365	{
	366	//
	367	// Draw the fitted function
	368	//
	369	TString drawOpt(option);
	370	drawOpt.ToLower();
	371
bc75eeb5	372	Float_t minY = 0.001;
	373	Float_t maxY = 1.2*fHistDataPM->GetMaximum();
	374	Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
	375	Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
	376	Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1)); // in MeV
	377	Float_t minMinY = fHistSignal->GetMinimum();
	378
48609e3d	379	TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
bc75eeb5	380	cSub->SetLeftMargin(0.15);
	381	cSub->SetRightMargin(0.0);
	382	cSub->SetTopMargin(0.002);
	383	cSub->SetBottomMargin(0.0);
	384	cSub->Divide(2,2,0.,0.);
572b0139	385	cSub->Draw();
572b0139	386
bc75eeb5	387	TVirtualPad* pad = cSub->cd(1);
	388	pad->SetLeftMargin(0.15);
	389	pad->SetRightMargin(0.0);
	390	pad->SetTopMargin(0.005);
	391	pad->SetBottomMargin(0.0);
	392	TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
	393	range1->SetStats(kFALSE);
	394	range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
	395	range1->GetYaxis()->CenterTitle();
	396	range1->GetYaxis()->SetLabelSize(0.05);
	397	range1->GetYaxis()->SetTitleSize(0.06);
	398	range1->GetYaxis()->SetTitleOffset(0.8);
	399	range1->Draw();
	400	fHistDataPM->SetLineColor(1);
	401	fHistDataPM->SetLineWidth(2);
	402	// fHistDataPM->SetMarkerStyle(21);
	403	fHistDataPM->Draw("Psame");
	404	TLatex *latex = new TLatex();
	405	latex->SetNDC();
	406	latex->SetTextSize(0.05);
	407	latex->DrawLatex(0.2, 0.95, "Background un-substracted");
	408	TLine line;
	409	line.SetLineWidth(1);
	410	line.SetLineStyle(2);
	411	line.DrawLine(fIntMin, minY, fIntMin, maxY);
	412	line.DrawLine(fIntMax, minY, fIntMax, maxY);
	413
	414	pad = cSub->cd(2);
	415	pad->SetLeftMargin(0.);
	416	pad->SetRightMargin(0.005);
	417	pad->SetTopMargin(0.005);
	418	pad->SetBottomMargin(0.0);
	419	TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
	420	range2->SetStats(kFALSE);
	421	range2->Draw();
	422	fHistBackground->SetLineColor(4);
	423	fHistBackground->SetLineWidth(2);
	424	// fHistBackground->SetMarkerColor(4);
	425	// fHistBackground->SetMarkerStyle(6);
	426	fHistBackground->Draw("Psame");
	427	latex->DrawLatex(0.05, 0.95, "Like-sign background");
	428	line.DrawLine(fIntMin, minY, fIntMin, maxY);
	429	line.DrawLine(fIntMax, minY, fIntMax, maxY);
	430	TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
	431	legend->SetFillColor(0);
	432	legend->SetMargin(0.15);
	433	legend->AddEntry(fHistDataPM, "N_{+-}", "l");
	434	legend->AddEntry(fHistDataPP, "N_{++}", "l");
	435	legend->AddEntry(fHistDataMM, "N_{--}", "l");
	436	legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
	437	legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
	438	legend->Draw();
	439
572b0139	440
bc75eeb5	441	pad = cSub->cd(3);
	442	pad->SetLeftMargin(0.15);
	443	pad->SetRightMargin(0.0);
	444	pad->SetTopMargin(0.0);
	445	pad->SetBottomMargin(0.15);
	446	TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
	447	range3->SetStats(kFALSE);
	448	range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
	449	range3->GetYaxis()->CenterTitle();
	450	range3->GetYaxis()->SetLabelSize(0.05);
	451	range3->GetYaxis()->SetTitleSize(0.06);
	452	range3->GetYaxis()->SetTitleOffset(0.8);
	453	range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
	454	range3->GetXaxis()->CenterTitle();
	455	range3->GetXaxis()->SetLabelSize(0.05);
	456	range3->GetXaxis()->SetTitleSize(0.06);
	457	range3->GetXaxis()->SetTitleOffset(1.0);
	458	range3->Draw();
	459	fHistDataPM->Draw("Psame");
	460	fHistDataPP->SetLineWidth(2);
	461	fHistDataPP->SetLineColor(6);
	462	fHistDataMM->SetLineWidth(2);
	463	fHistDataMM->SetLineColor(8);
	464	fHistDataPP->Draw("Psame");
	465	fHistDataMM->Draw("Psame");
	466	line.DrawLine(minX, 0.,maxX, 0.);
	467	line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
	468	line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
	469
	470	pad = cSub->cd(4);
	471	pad->SetLeftMargin(0.0);
	472	pad->SetRightMargin(0.005);
	473	pad->SetTopMargin(0.0);
	474	pad->SetBottomMargin(0.15);
	475	TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
	476	range4->SetStats(kFALSE);
	477	range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
	478	range4->GetXaxis()->CenterTitle();
	479	range4->GetXaxis()->SetLabelSize(0.05);
	480	range4->GetXaxis()->SetTitleSize(0.06);
	481	range4->GetXaxis()->SetTitleOffset(1.0);
	482	range4->Draw();
	483	fHistSignal->SetLineWidth(2);
	484	fHistSignal->SetLineColor(2);
	485	fHistSignal->Draw("Psame");
	486	latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
	487	if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
	488	line.DrawLine(minX, 0.,maxX, 0.);
	489	line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
	490	line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
48609e3d	491
48609e3d	492	cSub->SaveAs(Form("%s_summary.png", fName.Data()));
572b0139	493	}
572b0139	494