SideBandSubtract Class Reference

#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 ()
	SideBandSubtract (RooAbsPdf *model_shape, RooAbsPdf *bkg_shape, RooDataSet *data, RooRealVar *sep_var, const std::vector< TH1F * > &base, bool verb)
	~SideBandSubtract ()

Private Member Functions
Double_t	getYield (const 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 317 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,
		const std::vector< TH1F * > &	base,
		bool	verb
	)

Definition at line 334 of file SideBandSubtraction.cc.

  : BackgroundPDF(bkg_shape), 
    ModelPDF(model_shape), 
    Data(data),
    SeparationVariable(sep_var),
    verbose(verb),
    SignalRegions(),
    SideBandRegions(),
    SideBandHistos(),
    RawHistos(),
    SBSHistos(),
    BaseHistos(base),
    base_histo(0),
    fit_result(0),
    SignalSidebandRatio(0)
{ }

SideBandSubtract::~SideBandSubtract

(

)

Definition at line 402 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 432 of file SideBandSubtraction.cc.

References SbsRegion::max, bookConverter::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 421 of file SideBandSubtraction.cc.

References SbsRegion::max, bookConverter::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 490 of file SideBandSubtraction.cc.

References 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 443 of file SideBandSubtraction.cc.

References BackgroundPDF, BaseHistos, ecal_dqm_sourceclient-live_cfg::cerr, gather_cfg::cout, Data, doSubtraction(), fit_result, getYield(), i, ModelPDF, NULL, resetSBSProducts(), SeparationVariable, SideBandRegions, SignalRegions, SignalSidebandRatio, and AlCaHLTBitMon_QueryRunRegistry::string.

{
  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 81 of file SideBandSubtraction.cc.

References ecal_dqm_sourceclient-live_cfg::cerr, i, j, bookConverter::max, min(), NULL, setHistOptions(), AlCaHLTBitMon_QueryRunRegistry::string, relativeConstraints::value, and MetTreeProducer::var().

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 532 of file SideBandSubtraction.cc.

References BaseHistos.

{
  return BaseHistos;
}

RooFitResult * SideBandSubtract::getFitResult

(

)

Definition at line 528 of file SideBandSubtraction.cc.

References fit_result.

{
  return fit_result;
}

vector< TH1F > SideBandSubtract::getRawHistos

(

)

Definition at line 536 of file SideBandSubtraction.cc.

References RawHistos.

{
  return RawHistos;
}

vector< TH1F > SideBandSubtract::getSBSHistos

(

)

Definition at line 540 of file SideBandSubtraction.cc.

References SBSHistos.

{
  return SBSHistos;
}

Double_t SideBandSubtract::getSTSRatio

(

)

Definition at line 544 of file SideBandSubtraction.cc.

References SignalSidebandRatio.

{
  return SignalSidebandRatio;
}

Double_t SideBandSubtract::getYield ( const std::vector< SbsRegion > &  Regions, RooAbsPdf *  PDF  )

private

Definition at line 52 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 299 of file SideBandSubtraction.cc.

References ecal_dqm_sourceclient-live_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 151 of file SideBandSubtraction.cc.

References ecal_dqm_sourceclient-live_cfg::cerr, gather_cfg::cout, lut2db_cfg::filename, i, NULL, GenerateHcalLaserBadRunList::out, convertSQLitetoXML_cfg::output, print_histo(), AlCaHLTBitMon_QueryRunRegistry::string, 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";
  std::ofstream output(result_outname.c_str(),std::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 548 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

)

Definition at line 214 of file SideBandSubtraction.cc.

References ecal_dqm_sourceclient-live_cfg::cerr, HILowLumiHLTOfflineSource_cfi::dirname, i, relval_steps::key, NULL, getGTfromDQMFile::obj, convertSQLitetoXML_cfg::output, and stringify().

{
  //saves the ephemeral stuff to a root file for future
  //use/modification (ie everything printed by printResults())

  TFile output(outname.c_str(),"UPDATE"); //open the output file,
                      //create it if it doesn't
                      //exist
  //Since keys are only available from files on disk, we need to write
  //out a new file.  If the file already existed, then we opened to
  //update, and are writing nothing new.  
  output.Write();
  TString dirname;
  TIter nextkey(output.GetListOfKeys());
  TKey *key;
  TDirectory* curDir=NULL;
  while((key=(TKey*)nextkey.Next()))
    {

      if(key==NULL)
    break;
      TObject *obj = key->ReadObj();
      if(obj->IsA()->InheritsFrom("TDirectory"))
      dirname=obj->GetName();

    }

  if(dirname=="")
    {
      //we didn't find any directories so, we'll make a new one
      curDir = output.mkdir("run0","Run 0");
      output.cd("run0");
    }
  else
    {
      //manipulate last dir string, make a new one, and get ready to fill
      dirname.Remove(0,3);
      Int_t run_num = dirname.Atoi();
      run_num++;
      dirname = "run" + stringify(run_num);
      curDir = output.mkdir(dirname.Data(),("Run "+stringify(run_num)).c_str());
      output.cd(dirname.Data());
    }
  if(curDir==NULL)
    curDir = output.GetDirectory("",kTRUE,"");

  //these should all be the same size, but to be pedantic we'll loop
  //over each one individually, also, we need to associate them with
  //the directory because by default they "float" in memory to avoid
  //conflicts with other root files the user has open. If they want to
  //write to those files, they need to close their file, pass the name
  //here, and then let us work.
  for(unsigned int i=0; i < SideBandHistos.size(); ++i)
    {
      SideBandHistos[i].SetDirectory(curDir);
      SideBandHistos[i].Write();
    }
  for(unsigned int i=0; i < RawHistos.size(); ++i)
    {
      RawHistos[i].SetDirectory(curDir);
      RawHistos[i].Write();
    }
  for(unsigned int i=0; i < SBSHistos.size(); ++i)
    {
      SBSHistos[i].SetDirectory(curDir);
      SBSHistos[i].Write();
    }
  if(Data!=NULL && ModelPDF!=NULL && BackgroundPDF!=NULL && SeparationVariable!=NULL)
    {
      RooPlot *sep_varFrame = SeparationVariable->frame();
      Data->plotOn(sep_varFrame);
      ModelPDF->plotOn(sep_varFrame);
      BackgroundPDF->plotOn(sep_varFrame);
      sep_varFrame->Write();
    }
  else
    cerr <<"ERROR: saveResults, did not save RooPlot of data and fit\n";
  output.Write();

}

void SideBandSubtract::setDataSet

(

RooDataSet *

newData

)

Definition at line 294 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.

Referenced by doGlobalFit().

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.