[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronSignalBase.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 SignalBase                                  //
//                                                                       //
//                                                                       //
/*
Base class for signal extraction from a histogram or an array of histograms
The histogram is assumed to be an inv. mass spectrum,
the array of histograms is assumed to be an array with inv. mass histograms
resulting from single and mixed events, as defined in AliDielectron.cxx

*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////


#include <TVectorT.h>
#include <TPaveText.h>
#include <TF1.h>
#include <TH1.h>
#include <TDatabasePDG.h>

#include "AliDielectronSignalFunc.h"
#include "AliDielectronSignalBase.h"

ClassImp(AliDielectronSignalBase)

TH1F* AliDielectronSignalBase::fgHistSimPM=0x0;
TObject* AliDielectronSignalBase::fgPeakShape=0x0;
const char* AliDielectronSignalBase::fgkValueNames[6] = {
  "Signal","Background","Significance","Signal/Background","Mass","MassWidth"};

AliDielectronSignalBase::AliDielectronSignalBase() :
  TNamed(),
  fHistSignal(0),
  fHistBackground(0),
  fHistDataPM(0),
  fHistDataPP(0),
  fHistDataMM(0),
  fHistDataME(0),
  fHistRfactor(0),
  fValues(6),
  fErrors(6),
  fIntMin(0),
  fIntMax(0),
  fFitMin(0),
  fFitMax(0),
  fRebin(1),
  fMethod(kLikeSign),
  fScaleMin(0.),
  fScaleMax(0.),
  fScaleMin2(0.),
  fScaleMax2(0.),
  fScaleFactor(1.),
  fMixingCorr(kFALSE),
  fPeakMethod(kBinCounting),
  fProcessed(kFALSE),
  fPOIpdg(443)
{
  //
  // Default Constructor
  //
}

//______________________________________________
AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) :
  TNamed(name, title),
  fHistSignal(0),
  fHistBackground(0),
  fHistDataPM(0),
  fHistDataPP(0),
  fHistDataMM(0),
  fHistDataME(0),
  fHistRfactor(0),
  fValues(6),
  fErrors(6),
  fIntMin(0),
  fIntMax(0),
  fFitMin(0),
  fFitMax(0),
  fRebin(1),
  fMethod(kLikeSign),
  fScaleMin(0.),
  fScaleMax(0.),
  fScaleMin2(0.),
  fScaleMax2(0.),
  fScaleFactor(1.),
  fMixingCorr(kFALSE),
  fPeakMethod(kBinCounting),
  fProcessed(kFALSE),
  fPOIpdg(443)
{
  //
  // Named Constructor
  //
}

//______________________________________________
AliDielectronSignalBase::~AliDielectronSignalBase()
{
  //
  // Default Destructor
  //
  if(fHistSignal) delete fHistSignal;
  if(fHistBackground) delete fHistBackground;
  if (fHistDataPP) delete fHistDataPP;
  if (fHistDataPM) delete fHistDataPM;
  if (fHistDataMM) delete fHistDataMM;
  if (fHistDataME) delete fHistDataME;
  if (fHistRfactor)delete fHistRfactor;  
}

//______________________________________________
TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/*=0.*/, Double_t y1/*=0.*/, Double_t x2/*=0.*/, Double_t y2/*=0.*/)
{
  //
  // Draw extracted values in a TPaveText
  // with the corners x1,y2,x2,y2
  //
  if (TMath::Abs(x1)<1e-20&&TMath::Abs(x2)<1e-20){
    x1=.6;
    x2=.9;
    y1=.7;
    y2=.9;
  }
  TPaveText *t=new TPaveText(x1,y1,x2,y2,"brNDC");
  t->SetFillColor(kWhite);
  t->SetBorderSize(1);
  t->SetTextAlign(12);
  t->AddText(Form("Range  : %.2f - %.2f GeV/c^{2}", fIntMin, fIntMax));
  t->AddText(Form("Signal : %.1f #pm %.1f", fValues(0), fErrors(0)));
  t->AddText(Form("Backgnd: %.1f #pm %.1f", fValues(1), fErrors(1)));
  t->AddText(Form("Signif.: %.2f #pm %.2f", fValues(2), fErrors(2)));
  t->AddText(Form("S/B    : %.2f #pm %.2f", fValues(3), fErrors(3)));
  if(fValues(4)>0)
    t->AddText(Form("Mass: %.2f #pm %.2f GeV/c^{2}", fValues(4), fErrors(4)));
  if(fValues(5)>0)
    t->AddText(Form("Mass res.: %.1f #pm %.1f MeV/c^{2}", 1000*fValues(5), 1000*fErrors(5)));
  t->Draw();

  return t;
}

//______________________________________________
void AliDielectronSignalBase::Print(Option_t */*option*/) const
{
  //
  // Print the statistics
  //
  printf("Signal : %.5g #pm %.5g\n",fValues(0), fErrors(0));
  printf("Backgnd: %.5g #pm %.5g\n",fValues(1), fErrors(1));
  printf("Signif.: %.5g #pm %.5g\n",fValues(2), fErrors(2));
  printf("SoB    : %.5g #pm %.5g\n",fValues(3), fErrors(3));
  if(fValues(4)>0)
    printf("Mass: %.5g #pm %.5g\n", fValues(4), fErrors(4));
  if(fValues(5)>0)
    printf("Mass res.: %.5g #pm %.5g\n", fValues(5), fErrors(5));
}

//______________________________________________
Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax)
{
  //
  // scale histBackground to match the integral of histRaw in the interval intMin, intMax
  //

  //protect using over and underflow bins in normalisation calculation
  if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
  if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
  
  if (intMax>histRaw->GetXaxis()->GetXmax())
    intMax=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
  if (intMax>histBackground->GetXaxis()->GetXmax())
    intMax=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
  
  Double_t intRaw  = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
  Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
  Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
  if (intBack>0){
    histBackground->Sumw2();
    histBackground->Scale(scaleFactor);
  }

  return scaleFactor;
}
//______________________________________________
Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax, Double_t intMin2, Double_t intMax2)
{
  //
  // scale histBackground to match the integral of histRaw in the interval intMin, intMax and intMin2, intMax2
  //

  if(intMin2==intMax2) return (ScaleHistograms(histRaw, histBackground, intMin, intMax));

  //protect using over and underflow bins in normalisation calculation
  if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
  if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
  
  if (intMax2>histRaw->GetXaxis()->GetXmax())
    intMax2=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
  if (intMax2>histBackground->GetXaxis()->GetXmax())
    intMax2=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
  
  Double_t intRaw  = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
  Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
  intRaw  += histRaw->Integral(histRaw->FindBin(intMin2),histRaw->FindBin(intMax2));
  intBack += histBackground->Integral(histBackground->FindBin(intMin2),histBackground->FindBin(intMax2));

  Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
  if (intBack>0){
    histBackground->Sumw2();
    histBackground->Scale(scaleFactor);
  }

  return scaleFactor;
}

TObject* AliDielectronSignalBase::DescribePeakShape(ESignalExtractionMethod method, Bool_t replaceValErr,  TH1F *mcShape) {
  //
  // Describe the extracted peak by the selected method and overwrite signal etc if needed
  //

  fPeakMethod=method;
  Double_t data=0.;
  Double_t mc=0.;
  Double_t massPOI=TDatabasePDG::Instance()->GetParticle(fPOIpdg)->Mass();
  Double_t nPOI     = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
  Double_t binWidth = fHistSignal->GetBinWidth(  fHistSignal->FindBin(massPOI));
  TF1 *fit=0x0;
  Int_t parMass =-1;
  Int_t parSigma=-1;

  // do the scaling/fitting
  switch(fPeakMethod) {
  case kBinCounting: /*nothing needs to be done*/ 
    break;

  case kMCScaledMax:
    if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
    data = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
    mc   = mcShape->GetBinContent(fHistSignal->FindBin(massPOI));
    mcShape->Scale(data / mc );
    break;

  case kMCScaledInt:
    if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
    if(mcShape->GetBinWidth(1)!=fHistSignal->GetBinWidth(1)) 
      printf(" WARNING: MC and signal histogram have different bin widths. \n");
    data = fHistSignal->Integral(fHistSignal->FindBin(fIntMin),fHistSignal->FindBin(fIntMax));
    mc   = mcShape->Integral(mcShape->FindBin(fIntMin),mcShape->FindBin(fIntMax));
    mcShape->Scale(data / mc );
    break;

  case kMCFitted:
    if(!mcShape && !fgHistSimPM) { printf(" ERROR: No MC histogram passed or set. Returning. \n"); return 0x0; }
    if(!fgHistSimPM) fgHistSimPM=mcShape;
    fit = new TF1("fitMC",AliDielectronSignalFunc::PeakFunMC,fFitMin,fFitMax,1);
    fit->SetParNames("N");
    fHistSignal->Fit(fit,"RNI0");
    break;

  case kCrystalBall:
    fit = new TF1("fitCB",AliDielectronSignalFunc::PeakFunCB,fFitMin,fFitMax,5);
    fit->SetParNames("alpha","n","meanx","sigma","N");
    //  fit->SetParameters(-.2,5.,gMjpsi,.06,20);
    //  fit->SetParameters(1.,3.6,gMjpsi,.08,700);
    fit->SetParameters(0.4, 4.0, massPOI, 0.025, 1.3*nPOI);
    fit->SetParLimits(0, 0.0,           1.           );
    fit->SetParLimits(1, 0.01,          10.          );
    fit->SetParLimits(2, massPOI-0.02,  massPOI+0.02 );
    fit->SetParLimits(3, 0.001,          0.2         );
    fit->SetParLimits(4, 0.2*nPOI,      2.0*nPOI     );
    parMass=2;
    parSigma=3;
    fHistSignal->Fit(fit,"RNI0");
    break;

  case kGaus:
    fit = new TF1("fitGaus",AliDielectronSignalFunc::PeakFunGaus,fFitMin,fFitMax,3);
    //fit = new TF1("fitGaus","gaus",fFitMin,fFitMax);
    fit->SetParNames("N","meanx","sigma");
    fit->SetParameters(1.3*nPOI, massPOI, 0.025);
    fit->SetParLimits(0, 0.2*nPOI,      2.0*nPOI     );
    fit->SetParLimits(1, massPOI-0.02, massPOI+0.02);
    fit->SetParLimits(2, 0.001,         1.           );
    parMass=1;
    parSigma=2;
    fHistSignal->Fit(fit,"RNI0");
    break;

  }

  // overwrite values and errors if requested
  if(replaceValErr) {
    switch(fPeakMethod) {
    case kBinCounting:
      fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), fHistSignal->FindBin(fIntMax), fErrors(0));
      SetSignificanceAndSOB();
      break;
    case kMCScaledMax:
    case kMCScaledInt:
      fValues(0) = mcShape->IntegralAndError(mcShape->FindBin(fIntMin), mcShape->FindBin(fIntMax), fErrors(0));
      SetSignificanceAndSOB();
      break;

    case kMCFitted:
    case kCrystalBall:
    case kGaus:
      fValues(0) = fit->Integral(fIntMin, fIntMax)/binWidth;
      fErrors(0) = fit->IntegralError(fIntMin, fIntMax)/binWidth;
      SetSignificanceAndSOB();
      break;
    }
    // set mass position and width
    if(parMass>=0) {
      fValues(4) = fit->GetParameter(parMass);
      fErrors(4) = fit->GetParError(parMass);
    }
    if(parSigma>=0) {
      fValues(5) = fit->GetParameter(parSigma);
      fErrors(5) = fit->GetParError(parSigma);
    }
    else {
      // calculate FWHM
      SetFWHM();
    }
  }

  // set the peak method obj
  switch(fPeakMethod) {
  case kBinCounting:
    fgPeakShape=(TH1F*)fHistSignal->Clone("BinCount");
    break;
  case kMCScaledMax:
  case kMCScaledInt:
    fgPeakShape=mcShape;
    break;
  case kMCFitted:
  case kCrystalBall:
  case kGaus:
    fgPeakShape=fit;
    break;
  }

  return fgPeakShape;

}
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	// Dielectron SignalBase //
	18	// //
	19	// //
	20	/*
	21	Base class for signal extraction from a histogram or an array of histograms
	22	The histogram is assumed to be an inv. mass spectrum,
	23	the array of histograms is assumed to be an array with inv. mass histograms
	24	resulting from single and mixed events, as defined in AliDielectron.cxx
	25
	26	*/
	27	// //
	28	///////////////////////////////////////////////////////////////////////////
	29
	30
48609e3d	31	#include <TVectorT.h>
48609e3d	32	#include <TPaveText.h>
37a6f270	33	#include <TF1.h>
	34	#include <TH1.h>
	35	#include <TDatabasePDG.h>
	36
	37	#include "AliDielectronSignalFunc.h"
48609e3d	38	#include "AliDielectronSignalBase.h"
572b0139	39
	40	ClassImp(AliDielectronSignalBase)
	41
37a6f270	42	TH1F* AliDielectronSignalBase::fgHistSimPM=0x0;
	43	TObject* AliDielectronSignalBase::fgPeakShape=0x0;
	44	const char* AliDielectronSignalBase::fgkValueNames[6] = {
	45	"Signal","Background","Significance","Signal/Background","Mass","MassWidth"};
	46
572b0139	47	AliDielectronSignalBase::AliDielectronSignalBase() :
572b0139	48	TNamed(),
bc75eeb5	49	fHistSignal(0),
	50	fHistBackground(0),
	51	fHistDataPM(0),
	52	fHistDataPP(0),
	53	fHistDataMM(0),
5720c765	54	fHistDataME(0),
ac390e40	55	fHistRfactor(0),
8df8e382	56	fValues(6),
8df8e382	57	fErrors(6),
bc75eeb5	58	fIntMin(0),
	59	fIntMax(0),
	60	fFitMin(0),
	61	fFitMax(0),
	62	fRebin(1),
	63	fMethod(kLikeSign),
554e40f8	64	fScaleMin(0.),
554e40f8	65	fScaleMax(0.),
fd6ebd85	66	fScaleMin2(0.),
fd6ebd85	67	fScaleMax2(0.),
ffbede40	68	fScaleFactor(1.),
ac390e40	69	fMixingCorr(kFALSE),
37a6f270	70	fPeakMethod(kBinCounting),
	71	fProcessed(kFALSE),
	72	fPOIpdg(443)
572b0139	73	{
	74	//
	75	// Default Constructor
	76	//
572b0139	77	}
	78
	79	//______________________________________________
	80	AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) :
	81	TNamed(name, title),
bc75eeb5	82	fHistSignal(0),
	83	fHistBackground(0),
	84	fHistDataPM(0),
	85	fHistDataPP(0),
	86	fHistDataMM(0),
5720c765	87	fHistDataME(0),
ac390e40	88	fHistRfactor(0),
8df8e382	89	fValues(6),
8df8e382	90	fErrors(6),
bc75eeb5	91	fIntMin(0),
	92	fIntMax(0),
	93	fFitMin(0),
	94	fFitMax(0),
	95	fRebin(1),
	96	fMethod(kLikeSign),
554e40f8	97	fScaleMin(0.),
554e40f8	98	fScaleMax(0.),
fd6ebd85	99	fScaleMin2(0.),
fd6ebd85	100	fScaleMax2(0.),
ffbede40	101	fScaleFactor(1.),
ac390e40	102	fMixingCorr(kFALSE),
37a6f270	103	fPeakMethod(kBinCounting),
	104	fProcessed(kFALSE),
	105	fPOIpdg(443)
572b0139	106	{
	107	//
	108	// Named Constructor
	109	//
	110	}
	111
	112	//______________________________________________
	113	AliDielectronSignalBase::~AliDielectronSignalBase()
	114	{
	115	//
	116	// Default Destructor
	117	//
bc75eeb5	118	if(fHistSignal) delete fHistSignal;
bc75eeb5	119	if(fHistBackground) delete fHistBackground;
554e40f8	120	if (fHistDataPP) delete fHistDataPP;
	121	if (fHistDataPM) delete fHistDataPM;
	122	if (fHistDataMM) delete fHistDataMM;
5720c765	123	if (fHistDataME) delete fHistDataME;
ac390e40	124	if (fHistRfactor)delete fHistRfactor;
572b0139	125	}
	126
	127	//______________________________________________
	128	TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/=0./, Double_t y1/=0./, Double_t x2/=0./, Double_t y2/=0./)
	129	{
	130	//
	131	// Draw extracted values in a TPaveText
	132	// with the corners x1,y2,x2,y2
	133	//
	134	if (TMath::Abs(x1)<1e-20&&TMath::Abs(x2)<1e-20){
	135	x1=.6;
	136	x2=.9;
	137	y1=.7;
	138	y2=.9;
	139	}
	140	TPaveText *t=new TPaveText(x1,y1,x2,y2,"brNDC");
	141	t->SetFillColor(kWhite);
	142	t->SetBorderSize(1);
	143	t->SetTextAlign(12);
bc75eeb5	144	t->AddText(Form("Range : %.2f - %.2f GeV/c^{2}", fIntMin, fIntMax));
	145	t->AddText(Form("Signal : %.1f #pm %.1f", fValues(0), fErrors(0)));
	146	t->AddText(Form("Backgnd: %.1f #pm %.1f", fValues(1), fErrors(1)));
	147	t->AddText(Form("Signif.: %.2f #pm %.2f", fValues(2), fErrors(2)));
	148	t->AddText(Form("S/B : %.2f #pm %.2f", fValues(3), fErrors(3)));
37a6f270	149	if(fValues(4)>0)
bc75eeb5	150	t->AddText(Form("Mass: %.2f #pm %.2f GeV/c^{2}", fValues(4), fErrors(4)));
37a6f270	151	if(fValues(5)>0)
bc75eeb5	152	t->AddText(Form("Mass res.: %.1f #pm %.1f MeV/c^{2}", 1000fValues(5), 1000fErrors(5)));
572b0139	153	t->Draw();
	154
	155	return t;
	156	}
	157
9143d69f	158	//______________________________________________
	159	void AliDielectronSignalBase::Print(Option_t /option*/) const
	160	{
	161	//
	162	// Print the statistics
	163	//
bc75eeb5	164	printf("Signal : %.5g #pm %.5g\n",fValues(0), fErrors(0));
	165	printf("Backgnd: %.5g #pm %.5g\n",fValues(1), fErrors(1));
	166	printf("Signif.: %.5g #pm %.5g\n",fValues(2), fErrors(2));
	167	printf("SoB : %.5g #pm %.5g\n",fValues(3), fErrors(3));
37a6f270	168	if(fValues(4)>0)
bc75eeb5	169	printf("Mass: %.5g #pm %.5g\n", fValues(4), fErrors(4));
37a6f270	170	if(fValues(5)>0)
bc75eeb5	171	printf("Mass res.: %.5g #pm %.5g\n", fValues(5), fErrors(5));
9143d69f	172	}
2a14a7b1	173
2a14a7b1	174	//______________________________________________
554e40f8	175	Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax)
2a14a7b1	176	{
	177	//
	178	// scale histBackground to match the integral of histRaw in the interval intMin, intMax
	179	//
927480a1	180
ffbede40	181	//protect using over and underflow bins in normalisation calculation
	182	if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
	183	if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
	184
	185	if (intMax>histRaw->GetXaxis()->GetXmax())
	186	intMax=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
	187	if (intMax>histBackground->GetXaxis()->GetXmax())
	188	intMax=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
2a14a7b1	189
	190	Double_t intRaw = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
	191	Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
554e40f8	192	Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
2a14a7b1	193	if (intBack>0){
2a14a7b1	194	histBackground->Sumw2();
554e40f8	195	histBackground->Scale(scaleFactor);
2a14a7b1	196	}
554e40f8	197
554e40f8	198	return scaleFactor;
2a14a7b1	199	}
fd6ebd85	200	//______________________________________________
	201	Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax, Double_t intMin2, Double_t intMax2)
	202	{
	203	//
	204	// scale histBackground to match the integral of histRaw in the interval intMin, intMax and intMin2, intMax2
	205	//
927480a1	206
fd6ebd85	207	if(intMin2==intMax2) return (ScaleHistograms(histRaw, histBackground, intMin, intMax));
	208
	209	//protect using over and underflow bins in normalisation calculation
	210	if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
	211	if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
	212
	213	if (intMax2>histRaw->GetXaxis()->GetXmax())
	214	intMax2=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
	215	if (intMax2>histBackground->GetXaxis()->GetXmax())
	216	intMax2=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
	217
	218	Double_t intRaw = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
	219	Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
	220	intRaw += histRaw->Integral(histRaw->FindBin(intMin2),histRaw->FindBin(intMax2));
	221	intBack += histBackground->Integral(histBackground->FindBin(intMin2),histBackground->FindBin(intMax2));
	222
	223	Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
	224	if (intBack>0){
	225	histBackground->Sumw2();
	226	histBackground->Scale(scaleFactor);
	227	}
	228
	229	return scaleFactor;
	230	}
37a6f270	231
	232	TObject* AliDielectronSignalBase::DescribePeakShape(ESignalExtractionMethod method, Bool_t replaceValErr, TH1F *mcShape) {
	233	//
	234	// Describe the extracted peak by the selected method and overwrite signal etc if needed
	235	//
	236
	237	fPeakMethod=method;
	238	Double_t data=0.;
	239	Double_t mc=0.;
	240	Double_t massPOI=TDatabasePDG::Instance()->GetParticle(fPOIpdg)->Mass();
	241	Double_t nPOI = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
	242	Double_t binWidth = fHistSignal->GetBinWidth( fHistSignal->FindBin(massPOI));
	243	TF1 *fit=0x0;
	244	Int_t parMass =-1;
	245	Int_t parSigma=-1;
	246
	247	// do the scaling/fitting
	248	switch(fPeakMethod) {
	249	case kBinCounting: /nothing needs to be done/
	250	break;
	251
	252	case kMCScaledMax:
	253	if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
	254	data = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
	255	mc = mcShape->GetBinContent(fHistSignal->FindBin(massPOI));
	256	mcShape->Scale(data / mc );
	257	break;
	258
	259	case kMCScaledInt:
	260	if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
	261	if(mcShape->GetBinWidth(1)!=fHistSignal->GetBinWidth(1))
	262	printf(" WARNING: MC and signal histogram have different bin widths. \n");
	263	data = fHistSignal->Integral(fHistSignal->FindBin(fIntMin),fHistSignal->FindBin(fIntMax));
	264	mc = mcShape->Integral(mcShape->FindBin(fIntMin),mcShape->FindBin(fIntMax));
	265	mcShape->Scale(data / mc );
	266	break;
	267
	268	case kMCFitted:
	269	if(!mcShape && !fgHistSimPM) { printf(" ERROR: No MC histogram passed or set. Returning. \n"); return 0x0; }
	270	if(!fgHistSimPM) fgHistSimPM=mcShape;
	271	fit = new TF1("fitMC",AliDielectronSignalFunc::PeakFunMC,fFitMin,fFitMax,1);
	272	fit->SetParNames("N");
	273	fHistSignal->Fit(fit,"RNI0");
	274	break;
	275
	276	case kCrystalBall:
	277	fit = new TF1("fitCB",AliDielectronSignalFunc::PeakFunCB,fFitMin,fFitMax,5);
	278	fit->SetParNames("alpha","n","meanx","sigma","N");
	279	// fit->SetParameters(-.2,5.,gMjpsi,.06,20);
	280	// fit->SetParameters(1.,3.6,gMjpsi,.08,700);
	281	fit->SetParameters(0.4, 4.0, massPOI, 0.025, 1.3*nPOI);
	282	fit->SetParLimits(0, 0.0, 1. );
	283	fit->SetParLimits(1, 0.01, 10. );
	284	fit->SetParLimits(2, massPOI-0.02, massPOI+0.02 );
	285	fit->SetParLimits(3, 0.001, 0.2 );
	286	fit->SetParLimits(4, 0.2nPOI, 2.0nPOI );
	287	parMass=2;
	288	parSigma=3;
	289	fHistSignal->Fit(fit,"RNI0");
	290	break;
	291
	292	case kGaus:
	293	fit = new TF1("fitGaus",AliDielectronSignalFunc::PeakFunGaus,fFitMin,fFitMax,3);
	294	//fit = new TF1("fitGaus","gaus",fFitMin,fFitMax);
295	fit->SetParNames("N","meanx","sigma");
296	fit->SetParameters(1.3*nPOI, massPOI, 0.025);
297	fit->SetParLimits(0, 0.2nPOI, 2.0nPOI );
298	fit->SetParLimits(1, massPOI-0.02, massPOI+0.02);
299	fit->SetParLimits(2, 0.001, 1. );
300	parMass=1;
301	parSigma=2;
302	fHistSignal->Fit(fit,"RNI0");
303	break;
304
305	}
306
307	// overwrite values and errors if requested
308	if(replaceValErr) {
309	switch(fPeakMethod) {
310	case kBinCounting:
311	fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), fHistSignal->FindBin(fIntMax), fErrors(0));
312	SetSignificanceAndSOB();
313	break;
314	case kMCScaledMax:
315	case kMCScaledInt:
316	fValues(0) = mcShape->IntegralAndError(mcShape->FindBin(fIntMin), mcShape->FindBin(fIntMax), fErrors(0));
317	SetSignificanceAndSOB();
318	break;
319
320	case kMCFitted:
321	case kCrystalBall:
322	case kGaus:
323	fValues(0) = fit->Integral(fIntMin, fIntMax)/binWidth;
324	fErrors(0) = fit->IntegralError(fIntMin, fIntMax)/binWidth;
325	SetSignificanceAndSOB();
326	break;
327	}
328	// set mass position and width
329	if(parMass>=0) {
330	fValues(4) = fit->GetParameter(parMass);
331	fErrors(4) = fit->GetParError(parMass);
332	}
333	if(parSigma>=0) {
334	fValues(5) = fit->GetParameter(parSigma);
335	fErrors(5) = fit->GetParError(parSigma);
336	}
337	else {
338	// calculate FWHM
339	SetFWHM();
340	}
341	}
342
343	// set the peak method obj
344	switch(fPeakMethod) {
345	case kBinCounting:
346	fgPeakShape=(TH1F*)fHistSignal->Clone("BinCount");
347	break;
348	case kMCScaledMax:
349	case kMCScaledInt:
350	fgPeakShape=mcShape;
351	break;
352	case kMCFitted:
353	case kCrystalBall:
354	case kGaus:
355	fgPeakShape=fit;
356	break;
357	}
358
359	return fgPeakShape;
360
361	}