#include <SideBandSubtraction.h>

Public Member Functions
void	addSideBandRegion (Double_t min, Double_t max)
void	addSignalRegion (Double_t min, Double_t max)
void	doFastSubtraction (TH1F &Total, TH1F &Result, SbsRegion &leftRegion, SbsRegion &rightRegion)
int	doGlobalFit ()
int	doSubtraction (RooRealVar *variable, Double_t stsratio, Int_t index)
std::vector< TH1F * >	getBaseHistos ()
RooFitResult *	getFitResult ()
std::vector< TH1F >	getRawHistos ()
std::vector< TH1F >	getSBSHistos ()
Double_t	getSTSRatio ()
void	printResults (std::string prefix="")
void	resetSBSProducts ()
void	saveResults (std::string outname)
void	setDataSet (RooDataSet *newData)
	SideBandSubtract (RooAbsPdf model_shape, RooAbsPdf bkg_shape, RooDataSet data, RooRealVar sep_var, std::vector< TH1F * > base, bool verb)
	SideBandSubtract ()
	~SideBandSubtract ()
Private Member Functions
Double_t	getYield (std::vector< SbsRegion > Regions, RooAbsPdf *PDF)
void	print_plot (RooRealVar *printVar, std::string outname)
Private Attributes
RooAbsPdf *	BackgroundPDF
TH1F *	base_histo
std::vector< TH1F * >	BaseHistos
RooDataSet *	Data
RooFitResult *	fit_result
RooAbsPdf *	ModelPDF
std::vector< TH1F >	RawHistos
std::vector< TH1F >	SBSHistos
RooRealVar *	SeparationVariable
std::vector< TH1F >	SideBandHistos
std::vector< SbsRegion >	SideBandRegions
std::vector< SbsRegion >	SignalRegions
Double_t	SignalSidebandRatio
bool	verbose

Detailed Description

Definition at line 20 of file SideBandSubtraction.h.

Constructor & Destructor Documentation

SideBandSubtract::SideBandSubtract ( )

Definition at line 315 of file SideBandSubtraction.cc.

  : BackgroundPDF(0), 
    ModelPDF(0), 
    Data(0),
    SeparationVariable(0),
    verbose(0),
    SignalRegions(),
    SideBandRegions(),
    SideBandHistos(0),
    RawHistos(0),
    SBSHistos(0),
    BaseHistos(0),
    fit_result(0),
    SignalSidebandRatio(0)
{
  // Default constructor so we can do fast subtraction
}

SideBandSubtract::SideBandSubtract	(	RooAbsPdf *	model_shape,
		RooAbsPdf *	bkg_shape,
		RooDataSet *	data,
		RooRealVar *	sep_var,
		std::vector< TH1F * >	base,
		bool	verb
	)

SideBandSubtract::~SideBandSubtract ( )

Definition at line 400 of file SideBandSubtraction.cc.

References resetSBSProducts().

{
  //**WARNING** 

  // We don't delete objects that we don't own (duh) so, all of our
  // pointers just hang out and get handled by other people :)
 
  // We DO own the BaseHistos IFF the user didn't provide us with 
  // them and we constructed them from the dataset... in this case 
  // base_histo will be non-null and we want to delete the vector 
  // of histograms...
  resetSBSProducts();
  /*
  if(base_histo!=NULL)
    {
      BaseHistos.clear();
    }
  */
}

Member Function Documentation

void SideBandSubtract::addSideBandRegion	(	Double_t	min,
		Double_t	max
	)

Definition at line 430 of file SideBandSubtraction.cc.

References max(), SbsRegion::max, SbsRegion::min, min, NULL, SbsRegion::RegionName, SeparationVariable, SideBandRegions, and stringify().

{
  SbsRegion sideband;
  sideband.min=min;
  sideband.max=max;
  sideband.RegionName="SideBand" + stringify(SideBandRegions.size());
  if(SeparationVariable!=NULL)
    SeparationVariable->setRange(sideband.RegionName.c_str(),sideband.min,sideband.max);
  SideBandRegions.push_back(sideband);
  return;
}

void SideBandSubtract::addSignalRegion	(	Double_t	min,
		Double_t	max
	)

Definition at line 419 of file SideBandSubtraction.cc.

References max(), SbsRegion::max, SbsRegion::min, min, NULL, SbsRegion::RegionName, SeparationVariable, SignalRegions, and stringify().

{
  SbsRegion signal;
  signal.min=min;
  signal.max=max;
  signal.RegionName="Signal" + stringify(SignalRegions.size());
  if(SeparationVariable!=NULL)
    SeparationVariable->setRange(signal.RegionName.c_str(),signal.min,signal.max);
  SignalRegions.push_back(signal);
  return;
}

void SideBandSubtract::doFastSubtraction	(	TH1F &	Total,
		TH1F &	Result,
		SbsRegion &	leftRegion,
		SbsRegion &	rightRegion
	)

Definition at line 488 of file SideBandSubtraction.cc.

References root::function(), SbsRegion::max, and SbsRegion::min.

{
  Int_t binMin = Total.FindBin(leftRegion.max,0.0,0.0);
  Int_t binMax = Total.FindBin(leftRegion.min,0.0,0.0);
  double numLeft = Total.Integral( binMin, binMax );

  binMin = Total.FindBin(rightRegion.max,0.0,0.0);
  binMax = Total.FindBin(rightRegion.min,0.0,0.0);
  double numRight = Total.Integral( binMin, binMax );
  
  const unsigned int nbinsx = Total.GetNbinsX();
  const double x1 = (leftRegion.max + leftRegion.min)/2.0;
  const double x2 = (rightRegion.max + rightRegion.min)/2.0;

  const double y1 = numLeft/(leftRegion.max - leftRegion.min);
  const double y2 = numRight/(rightRegion.max - rightRegion.min);
    
  const double Slope = (y2 - y1)/(x2 - x1);
  const double Intercept = y1 - Slope*x1;
  // Evantually we want to handle more complicated functions, but for
  // now... just use a linear one
  TF1 function("sbs_function_line", "[0]*x + [1]",Total.GetMinimum(), Total.GetMaximum());
  for ( unsigned int binx = 1;  binx <= nbinsx; ++binx ) 
    {
      double binWidth = Total.GetBinWidth(binx);
      function.SetParameter(0,Slope*binWidth);
      function.SetParameter(1,Intercept*binWidth);

      double xx = Total.GetBinCenter(binx);
      double cu = Total.GetBinContent(binx) - function.Eval(xx);    
      // TODO: Propogate the  error on the parameters in function.
      double error1 = Total.GetBinError(binx);
    
      Result.SetBinContent(binx, cu);
      Result.SetBinError(binx, error1);
    }
  Result.SetEntries(Result.Integral() );
}

int SideBandSubtract::doGlobalFit ( )

Definition at line 441 of file SideBandSubtraction.cc.

References BackgroundPDF, BaseHistos, benchmark_cfg::cerr, gather_cfg::cout, doSubtraction(), fit_result, getYield(), i, ModelPDF, NULL, resetSBSProducts(), SeparationVariable, SideBandRegions, SignalRegions, and SignalSidebandRatio.

{
  if(verbose)
    cout <<"Beginning SideBand Subtraction\n";

  if(ModelPDF!=NULL && Data!=NULL && SeparationVariable!=NULL)
    {
      fit_result = ModelPDF->fitTo(*Data,"r");
    }
  else
    {
      cerr <<"ERROR: doGobalFit, no ModelPDF, SeparationVariable or Data specified\n";
      return -1;
    }

  Double_t SideBandYield=getYield(SideBandRegions,BackgroundPDF);
  Double_t  BackgroundInSignal=getYield(SignalRegions,BackgroundPDF);

  SignalSidebandRatio = BackgroundInSignal/SideBandYield;
  if(verbose)
    {
      cout <<"Finished fitting background!\n";
      cout <<"Attained a Signal to Sideband ratio of: " << SignalSidebandRatio<<endl;
    }
  //I can't see a way around a double loop for each variable.  Maybe I
  //can get a C++/RooFit guru to save me the trouble here?


  //need to grab sbs objects after each global fit, because they get reset
  resetSBSProducts();
  TIterator* iter = (TIterator*) Data->get()->createIterator();
  RooAbsArg *variable;
  while((variable = (RooAbsArg*)iter->Next()))
    {
      for(unsigned int i=0; i < BaseHistos.size(); i++)
        {
          if((string)variable->GetName()!=(string)SeparationVariable->GetName() 
             && (string)variable->GetName()==(string)BaseHistos[i]->GetName())
            doSubtraction((RooRealVar*)variable,SignalSidebandRatio,i);
        }
    }

  //  clean up our memory...
  if(variable)      delete variable;
  if(iter)          delete iter;
  return 0;
}

int SideBandSubtract::doSubtraction	(	RooRealVar *	variable,
		Double_t	stsratio,
		Int_t	index
	)

Definition at line 79 of file SideBandSubtraction.cc.

References benchmark_cfg::cerr, i, j, max(), min, NULL, setHistOptions(), and relativeConstraints::value.

Referenced by doGlobalFit().

{
  if(Data==NULL || SeparationVariable==NULL)
    {
      cerr << "ERROR: Data or SeparationVariable is NULL returning now!\n";
      return -1;
    }
  TH1F* SideBandHist = (TH1F*)BaseHistos[index]->Clone();
  setHistOptions(SideBandHist,(string)variable->GetName()+"Sideband",(string)SideBandHist->GetTitle() + " Sideband",(string)variable->getUnit());

  TH1F* SignalHist = (TH1F*)BaseHistos[index]->Clone();
  setHistOptions(SignalHist,(string)variable->GetName()+"SignalHist",(string)SignalHist->GetTitle() + " Raw Signal",(string)variable->getUnit());

  //Begin a loop over the data to fill our histograms. I should figure
  //out how to do this in one shot to avoid a loop
  //O(N_vars*N_events)...

  TIterator* iter = (TIterator*) Data->get()->createIterator();
  RooAbsArg *var=NULL;
  RooRealVar *sep_var=NULL;
  while((var = (RooAbsArg*)iter->Next()))
    {
      if((string)var->GetName()==(string)SeparationVariable->GetName())
        {
          sep_var = (RooRealVar*)var;
          break;
        }
    }
  for(int i=0; i < Data->numEntries(); i++)
    {
      Data->get(i);
      Double_t value = variable->getVal();
      Double_t cutval = sep_var->getVal();

      for(unsigned int j=0; j < SideBandRegions.size(); j++) //UGLY!  Find better solution!
        {
          if(cutval > SideBandRegions[j].min && cutval < SideBandRegions[j].max)
            SideBandHist->Fill(value);
        }
      for(unsigned int j=0; j < SignalRegions.size(); j++)
        {
          if(cutval > SignalRegions[j].min && cutval < SignalRegions[j].max)
            SignalHist->Fill(value);
        }
    }
  //Save pre-subtracted histo
  SignalHist->Sumw2(); SideBandHist->Sumw2(); 
  //SignalHist->SetDirectory(0); SideBandHist->SetDirectory(0);
  RawHistos.push_back(*SignalHist);

  SignalHist->Add(SideBandHist, -stsratio);

  SignalHist->SetTitle(((string)BaseHistos[index]->GetTitle() + " SBS Signal").c_str());
  //Save subtracted histo
  SBSHistos.push_back(*SignalHist);
  //Save Sideband histo
  SideBandHistos.push_back(*SideBandHist);

  if(SideBandHist) delete SideBandHist;
  return 0;
}

vector< TH1F * > SideBandSubtract::getBaseHistos ( )

Definition at line 530 of file SideBandSubtraction.cc.

References BaseHistos.

{
  return BaseHistos;
}

RooFitResult * SideBandSubtract::getFitResult ( )

Definition at line 526 of file SideBandSubtraction.cc.

References fit_result.

{
  return fit_result;
}

vector< TH1F > SideBandSubtract::getRawHistos ( )

Definition at line 534 of file SideBandSubtraction.cc.

References RawHistos.

{
  return RawHistos;
}

vector< TH1F > SideBandSubtract::getSBSHistos ( )

Definition at line 538 of file SideBandSubtraction.cc.

References SBSHistos.

{
  return SBSHistos;
}

Double_t SideBandSubtract::getSTSRatio ( )

Definition at line 542 of file SideBandSubtraction.cc.

References SignalSidebandRatio.

{
  return SignalSidebandRatio;
}

Double_t SideBandSubtract::getYield	(	std::vector< SbsRegion >	Regions,
		RooAbsPdf *	PDF
	)		`[private]`

Definition at line 50 of file SideBandSubtraction.cc.

References gather_cfg::cout, i, and NULL.

Referenced by doGlobalFit().

{
  if(PDF==NULL || SeparationVariable==NULL)
    return 0.0;

  Double_t yield=0;
  RooAbsReal* intPDF;
  for(unsigned int i=0; i < Regions.size(); i++)
    {
      if(verbose)
        cout<<"Integrating over Range: "<<Regions[i].RegionName<<" from "<<Regions[i].min<<" to "<<Regions[i].max<<endl;
      intPDF = PDF->createIntegral(*SeparationVariable, Range(Regions[i].RegionName.c_str()));
      yield += intPDF->getVal();
      if(verbose)
        cout<<"Current yield: "<<yield<<endl;
    }
  return yield;
}

void SideBandSubtract::print_plot	(	RooRealVar *	printVar,
		std::string	outname
	)		`[private]`

Definition at line 297 of file SideBandSubtraction.cc.

References benchmark_cfg::cerr, and NULL.

{
  if(Data==NULL || ModelPDF==NULL)
    {
      cerr << "ERROR: print_plot, Data or ModelPDF are NULL\n";
      return;
    }

  RooPlot *genericFrame = printVar->frame();
  Data->plotOn(genericFrame);
  ModelPDF->plotOn(genericFrame);
  TCanvas genericCanvas;
  genericFrame->Draw();
  outname = outname + ".eps";
  genericCanvas.SaveAs(outname.c_str());
  outname.replace(outname.size()-3,3,"gif");
  genericCanvas.SaveAs(outname.c_str());
}

void SideBandSubtract::printResults ( std::string prefix = "" )

Definition at line 149 of file SideBandSubtraction.cc.

References benchmark_cfg::cerr, gather_cfg::cout, lut2db_cfg::filename, i, NULL, dbtoconf::out, convertSQLitetoXML_cfg::output, print_histo(), and EcalTangentSkim_cfg::Verbose.

{//handles *all* printing
  //spool over vectors of histograms and print them, then print
  //separation variable plots and the results text file.
  if(SeparationVariable==NULL)
    {
      cerr << "ERROR: printResults, SeparationVariable is NULL!\n";
      return;
    }
  string filename; //output file name
  for(unsigned int i=0; i < RawHistos.size(); ++i)
    {
      filename=prefix + "Raw_" + (string)RawHistos[i].GetName();
      print_histo(&RawHistos[i], filename);
    }
  for(unsigned int i=0; i < SBSHistos.size(); ++i)
    {
      filename=prefix + "SBS_" + (string)RawHistos[i].GetName();
      print_histo(&SBSHistos[i], filename);
    }
  if(verbose)
    {
      for(unsigned int i=0; i < SideBandHistos.size(); ++i)
        {
          filename=prefix + "SideBand_" + (string)RawHistos[i].GetName();
          print_histo(&SideBandHistos[i], filename);
        }
    }

  string outname = prefix + (string)SeparationVariable->GetName() + "_fitted.eps";
  if(Data!=NULL && ModelPDF!=NULL)
    {
      RooPlot *SepVarFrame = SeparationVariable->frame();
      Data->plotOn(SepVarFrame);
      ModelPDF->plotOn(SepVarFrame);
      TCanvas Canvas;
      SepVarFrame->Draw();
      Canvas.SaveAs(outname.c_str());
      outname.replace(outname.size()-3,3,"gif");
      Canvas.SaveAs(outname.c_str());
    }
  else
    cerr <<"ERROR: printResults, Data or ModelPDF is NULL!\n";

  string result_outname = prefix + "_fit_results.txt";
  ofstream output(result_outname.c_str(),ios::out);
  if(!output)
    {
      cout <<"ERROR: Could not open file for writing!\n";
      return;
    }

  if(fit_result!=NULL)
    {
#if ROOT_VERSION_CODE <= ROOT_VERSION(5,19,0)
      fit_result->printToStream(output,Verbose);
#else
      fit_result->printMultiline(output,kTRUE);
#endif
    }
}

void SideBandSubtract::resetSBSProducts ( )

Definition at line 546 of file SideBandSubtraction.cc.

References RawHistos, SBSHistos, and SideBandHistos.

Referenced by doGlobalFit(), and ~SideBandSubtract().

{
  //cout <<"Starting to reset products \n";

  if(!SideBandHistos.empty())
    SideBandHistos.clear();

  if(!RawHistos.empty())
    RawHistos.clear();
  if(!SBSHistos.empty())
    SBSHistos.clear();
}

void SideBandSubtract::saveResults ( std::string outname )

void SideBandSubtract::setDataSet ( RooDataSet * newData )

Definition at line 292 of file SideBandSubtraction.cc.

References NULL.

{
  if(newData!=NULL)
    Data=newData;
}

Member Data Documentation

RooAbsPdf* SideBandSubtract::BackgroundPDF [private]

Definition at line 25 of file SideBandSubtraction.h.

Referenced by doGlobalFit().

TH1F* SideBandSubtract::base_histo [private]

Definition at line 36 of file SideBandSubtraction.h.

std::vector<TH1F*> SideBandSubtract::BaseHistos [private]

Definition at line 35 of file SideBandSubtraction.h.

Referenced by doGlobalFit(), and getBaseHistos().

RooDataSet* SideBandSubtract::Data [private]

Definition at line 27 of file SideBandSubtraction.h.

RooFitResult* SideBandSubtract::fit_result [private]

Definition at line 37 of file SideBandSubtraction.h.

Referenced by doGlobalFit(), and getFitResult().

RooAbsPdf* SideBandSubtract::ModelPDF [private]

Definition at line 26 of file SideBandSubtraction.h.

Referenced by doGlobalFit().

std::vector<TH1F> SideBandSubtract::RawHistos [private]

Definition at line 33 of file SideBandSubtraction.h.

Referenced by getRawHistos(), and resetSBSProducts().

std::vector<TH1F> SideBandSubtract::SBSHistos [private]

Definition at line 34 of file SideBandSubtraction.h.

Referenced by getSBSHistos(), and resetSBSProducts().

RooRealVar* SideBandSubtract::SeparationVariable [private]

Definition at line 28 of file SideBandSubtraction.h.

Referenced by addSideBandRegion(), addSignalRegion(), and doGlobalFit().

std::vector<TH1F> SideBandSubtract::SideBandHistos [private]

Definition at line 32 of file SideBandSubtraction.h.

Referenced by resetSBSProducts().

std::vector<SbsRegion> SideBandSubtract::SideBandRegions [private]

Definition at line 31 of file SideBandSubtraction.h.

Referenced by addSideBandRegion(), and doGlobalFit().

std::vector<SbsRegion> SideBandSubtract::SignalRegions [private]

Definition at line 30 of file SideBandSubtraction.h.

Referenced by addSignalRegion(), and doGlobalFit().

Double_t SideBandSubtract::SignalSidebandRatio [private]

Definition at line 38 of file SideBandSubtraction.h.

Referenced by doGlobalFit(), and getSTSRatio().

bool SideBandSubtract::verbose [private]

Definition at line 29 of file SideBandSubtraction.h.

SideBandSubtract Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation