Reference Manual

Go to the documentation of this file.
00001 #include "RooGaussian.h"
00002 #include "RooChebychev.h"
00003 #include "RooVoigtian.h"
00004 #include "RooExponential.h"
00005 #include "RooRealVar.h"
00006 #include "RooFormulaVar.h"
00007 #include "RooDataHist.h"
00008 #include "RooAddPdf.h"
00009 #include "RooGlobalFunc.h"
00010 #include "RooCategory.h"
00011 #include "RooSimultaneous.h"
00012 #include "RooChi2Var.h"
00013 #include "TFile.h"
00014 #include "TH1F.h"
00015 #include "TH2F.h"
00016 #include "TH3F.h"
00017 #include "TProfile2D.h"
00018 #include "TCanvas.h"
00019 #include "RooPlot.h"
00020 
00021 namespace dqmTnP{
00022 
00023 class AbstractFitter{
00024   protected:
00025   RooRealVar mass;
00026   RooRealVar mean;
00027   double expectedMean;
00028   RooRealVar sigma;
00029   double expectedSigma;
00030   RooRealVar efficiency;
00031   RooRealVar nSignalAll;
00032   RooFormulaVar nSignalPass;
00033   RooFormulaVar nSignalFail;
00034   RooRealVar nBackgroundFail;
00035   RooRealVar nBackgroundPass;
00036   RooCategory category;
00037   RooSimultaneous simPdf;
00038   RooDataHist *data;
00039   double chi2;
00040   bool verbose;
00041 
00042   public:
00043   AbstractFitter(bool verbose_ = false):
00044     mass("mass","mass",0.,100.,"GeV"),
00045     mean("mean","mean",0.,100.,"GeV"),
00046     sigma("sigma","sigma",0.,100.,"GeV"),
00047     efficiency("efficiency","efficiency",0.5,0.0,1.0),
00048     nSignalAll("nSignalAll","nSignalAll",0.,1e10),
00049     nSignalPass("nSignalPass","nSignalAll*efficiency",RooArgList(nSignalAll,efficiency)),
00050     nSignalFail("nSignalFail","nSignalAll*(1-efficiency)",RooArgList(nSignalAll,efficiency)),
00051     nBackgroundFail("nBackgroundFail","nBackgroundFail",0.,1e10),
00052     nBackgroundPass("nBackgroundPass","nBackgroundPass",0.,1e10),
00053     category("category","category"),
00054     simPdf("simPdf","simPdf",category),
00055     data(0),
00056     verbose(verbose_)
00057   {
00058     //turn on/off default messaging of roofit
00059     RooMsgService::instance().setSilentMode( !verbose?kTRUE:kFALSE );
00060     for(int i=0; i<RooMsgService::instance().numStreams(); i++){
00061       RooMsgService::instance().setStreamStatus( i, verbose?kTRUE:kFALSE );
00062     }
00063     category.defineType("pass");
00064     category.defineType("fail");
00065   };
00066   virtual ~AbstractFitter(){};
00067   void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_){
00068     expectedMean = expectedMean_;
00069     expectedSigma = expectedSigma_;
00070     mass.setRange(massLow,massHigh);
00071     mean.setRange(massLow,massHigh);
00072   }
00073   virtual void fit(TH1* num, TH1* den) = 0;
00074   double getEfficiency(){ return efficiency.getVal(); }
00075   double getEfficiencyError(){ return efficiency.getError(); }
00076   double getChi2(){ return chi2; }
00077   void savePlot(TString name){
00078     using namespace RooFit;
00079     RooPlot* frame = mass.frame(Name(name), Title("Failing and Passing Probe Distributions"));
00080     data->plotOn(frame,Cut("category==category::pass"),LineColor(kGreen),MarkerColor(kGreen));
00081     data->plotOn(frame,Cut("category==category::fail"),LineColor(kRed),MarkerColor(kRed));
00082     simPdf.plotOn(frame,Slice(category,"pass"),ProjWData(category,*data),LineColor(kGreen));
00083     simPdf.plotOn(frame,Slice(category,"fail"),ProjWData(category,*data),LineColor(kRed));
00084     simPdf.paramOn(frame,Layout(0.58,0.99,0.99));
00085     data->statOn(frame,Layout(0.70,0.99,0.5));
00086     frame->Write();
00087     delete frame;
00088   }
00089 
00090   TString calculateEfficiency(TH3 *pass, TH3 *all, int massDimension, TProfile2D* &eff, TProfile2D* &effChi2, TString plotName=""){
00091     //sort out the TAxis
00092     TAxis *par1Axis, *par2Axis, *massAxis;
00093     int par1C, par2C, massC;
00094     if(massDimension==1){
00095       massAxis = all->GetXaxis();
00096       massC = 1;
00097       par1Axis = all->GetYaxis();
00098       par1C = all->GetXaxis()->GetNbins()+2;
00099       par2Axis = all->GetZaxis();
00100       par2C = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
00101     }else if(massDimension==2){
00102       par1Axis = all->GetXaxis();
00103       par1C = 1;
00104       massAxis = all->GetYaxis();
00105       massC = all->GetXaxis()->GetNbins()+2;
00106       par2Axis = all->GetZaxis();
00107       par2C = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
00108     }else if(massDimension==3){
00109       par1Axis = all->GetXaxis();
00110       par1C = 1;
00111       par2Axis = all->GetYaxis();
00112       par2C = all->GetXaxis()->GetNbins()+2;
00113       massAxis = all->GetZaxis();
00114       massC = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
00115     }else{
00116       return "massDimension > 3 !, skipping...";
00117     }
00118 
00119     //create eff and effChi2 TProfiles
00120     if(!par1Axis || !par2Axis) return "No par1Axis or par2Axis!";
00121     if( par1Axis->GetXbins()->GetSize()==0 && par2Axis->GetXbins()->GetSize()==0 ){
00122       eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax(),par2Axis->GetNbins(),par2Axis->GetXmin(),par2Axis->GetXmax());
00123     }else if( par1Axis->GetXbins()->GetSize()==0 ){
00124       eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax(),par2Axis->GetNbins(),par2Axis->GetXbins()->GetArray());
00125     }else if( par2Axis->GetXbins()->GetSize()==0 ){
00126       eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray(),par2Axis->GetNbins(),par2Axis->GetXmin(),par2Axis->GetXmax());
00127     }else{
00128       eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray(),par2Axis->GetNbins(),par2Axis->GetXbins()->GetArray());
00129     }
00130     eff->SetTitle("");
00131     eff->SetXTitle( par1Axis->GetTitle() );
00132     eff->SetYTitle( par2Axis->GetTitle() );
00133     eff->SetStats(kFALSE);
00134     effChi2 = (TProfile2D*)eff->Clone("efficiencyChi2");
00135     eff->SetZTitle("Efficiency");
00136     eff->SetOption("colztexte");
00137     eff->GetZaxis()->SetRangeUser(-0.001,1.001);
00138     effChi2->SetZTitle("Chi^2/NDF");
00139     effChi2->SetOption("colztext");
00140 
00141     //create the 1D mass distribution container histograms
00142     TH1D* all1D = (massAxis->GetXbins()->GetSize()==0) ?
00143       new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXmin(),massAxis->GetXmax()):
00144       new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXbins()->GetArray()); 
00145     TH1D* pass1D = (TH1D *)all1D->Clone("pass1D");
00146 
00147     //for each parameter bin fit the mass distributions
00148     for(int par1=1; par1<=par1Axis->GetNbins(); par1++){
00149       for(int par2=1; par2<=par2Axis->GetNbins(); par2++){
00150         for(int mass=1; mass<=massAxis->GetNbins(); mass++){
00151           int index = par1*par1C + par2*par2C + mass*massC;
00152           all1D->SetBinContent(mass,all->GetBinContent(index));
00153           pass1D->SetBinContent(mass,pass->GetBinContent(index));
00154         }
00155         fit( pass1D, all1D );
00156         int index = par1 + par2*(par1Axis->GetNbins()+2);
00157         eff->SetBinContent( index, getEfficiency() );
00158         eff->SetBinEntries( index, 1 );
00159         eff->SetBinError( index, sqrt( getEfficiency()*getEfficiency() + getEfficiencyError()*getEfficiencyError() ) );
00160         effChi2->SetBinContent( index, getChi2() );
00161         effChi2->SetBinEntries( index, 1 );
00162         if(plotName!=""){
00163           savePlot( TString::Format("%s_%d_%d",plotName.Data(),par1,par2) );
00164         }
00165       }
00166     }
00167     delete all1D;
00168     delete pass1D;
00169     return "";//OK
00170   }
00171 
00172   TString calculateEfficiency(TH2 *pass, TH2 *all, int massDimension, TProfile* &eff, TProfile* &effChi2, TString plotName=""){
00173     //sort out the TAxis
00174     TAxis *par1Axis, *massAxis;
00175     int par1C, massC;
00176     if(massDimension==1){
00177       massAxis = all->GetXaxis();
00178       massC = 1;
00179       par1Axis = all->GetYaxis();
00180       par1C = all->GetXaxis()->GetNbins()+2;
00181     }else if(massDimension==2){
00182       par1Axis = all->GetXaxis();
00183       par1C = 1;
00184       massAxis = all->GetYaxis();
00185       massC = all->GetXaxis()->GetNbins()+2;
00186     }else{
00187       return "massDimension > 2 !, skipping...";
00188     }
00189 
00190     //create eff and effChi2 TProfiles
00191     if(!par1Axis) return "No par1Axis!";
00192     eff = (par1Axis->GetXbins()->GetSize()==0)?
00193       new TProfile("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax()):
00194       new TProfile("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray());
00195     eff->SetTitle("");
00196     eff->SetXTitle( par1Axis->GetTitle() );
00197     eff->SetLineColor(2);
00198     eff->SetLineWidth(2);
00199     eff->SetMarkerStyle(20);
00200     eff->SetMarkerSize(0.8);
00201     eff->SetStats(kFALSE);
00202     effChi2 = (TProfile*)eff->Clone("efficiencyChi2");
00203     eff->SetYTitle("Efficiency");
00204     eff->SetOption("PE");
00205     eff->GetYaxis()->SetRangeUser(-0.001,1.001);
00206     effChi2->SetYTitle("Chi^2/NDF");
00207     effChi2->SetOption("HIST");
00208 
00209     //create the 1D mass distribution container histograms
00210     TH1D * all1D = (massAxis->GetXbins()->GetSize()==0) ?
00211       new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXmin(),massAxis->GetXmax()):
00212       new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXbins()->GetArray()); 
00213     TH1D * pass1D = (TH1D *)all1D->Clone("pass1D");
00214 
00215     //for each parameter bin fit the mass distributions
00216     for(int par1=1; par1<=par1Axis->GetNbins(); par1++){
00217       for(int mass=1; mass<=massAxis->GetNbins(); mass++){
00218         int index = par1*par1C + mass*massC;
00219         all1D->SetBinContent(mass,all->GetBinContent(index));
00220         pass1D->SetBinContent(mass,pass->GetBinContent(index));
00221       }
00222       fit( pass1D, all1D );
00223       int index = par1;
00224       eff->SetBinContent( index, getEfficiency() );
00225       eff->SetBinEntries( index, 1 );
00226       eff->SetBinError( index, sqrt( getEfficiency()*getEfficiency() + getEfficiencyError()*getEfficiencyError() ) );
00227       effChi2->SetBinContent( index, getChi2() );
00228       effChi2->SetBinEntries( index, 1 );
00229       if(plotName!=""){
00230         savePlot( TString::Format("%s_%d",plotName.Data(),par1) );
00231       }
00232     }
00233     delete all1D;
00234     delete pass1D;
00235     return "";//OK
00236   }
00237 };
00238 
00239 //concrete fitter: Gaussian signal plus linear background
00240 class GaussianPlusLinearFitter: public AbstractFitter{
00241   protected:
00242   RooGaussian gaussian;
00243   RooRealVar slopeFail;
00244   RooChebychev linearFail;
00245   RooRealVar slopePass;
00246   RooChebychev linearPass;
00247   RooAddPdf pdfFail;
00248   RooAddPdf pdfPass;
00249   public:
00250   GaussianPlusLinearFitter(bool verbose = false):
00251     AbstractFitter(verbose),
00252     gaussian("gaussian","gaussian",mass,mean,sigma),
00253     slopeFail("slopeFail","slopeFail",0.,-1.,1.),
00254     linearFail("linearFail","linearFail",mass,slopeFail),
00255     slopePass("slopePass","slopePass",0.,-1.,1.),
00256     linearPass("linearPass","linearPass",mass,slopePass),
00257     pdfFail("pdfFail","pdfFail", RooArgList(gaussian,linearFail), RooArgList(nSignalFail,nBackgroundFail)),
00258     pdfPass("pdfPass","pdfPass", RooArgList(gaussian,linearPass), RooArgList(nSignalPass,nBackgroundPass))
00259   {
00260     simPdf.addPdf(pdfFail,"fail");
00261     simPdf.addPdf(pdfPass,"pass");
00262   };
00263   ~GaussianPlusLinearFitter(){};
00264   void fit(TH1* pass, TH1* all){
00265     using namespace RooFit;
00266     all->Add(pass,-1);
00267     TH1* &fail = all;
00268     if(!data) delete data;
00269     data = new RooDataHist("data", "data", mass, Index(category), Import("fail",*fail), Import("pass",*pass) );
00270     if(pass->Integral()+fail->Integral() < 5){
00271       efficiency.setVal(0.5);
00272       efficiency.setError(0.5);
00273       chi2 = 0;
00274       return;
00275     }
00276     mean.setVal(expectedMean);
00277     sigma.setVal(expectedSigma);
00278     efficiency.setVal(pass->Integral()/(pass->Integral()+fail->Integral()));
00279     nSignalAll.setVal(0.5*(fail->Integral()+pass->Integral()));
00280     nBackgroundFail.setVal(0.5*fail->Integral());
00281     nBackgroundPass.setVal(0.5*pass->Integral());
00282     slopeFail.setVal(0.);
00283     slopePass.setVal(0.);
00284     if(verbose){
00285       simPdf.fitTo( *data );
00286     }else{
00287       simPdf.fitTo( *data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1) );
00288     }
00289     RooDataHist dataFail("fail", "fail", mass, fail );
00290     RooDataHist dataPass("pass", "pass", mass, pass );
00291     chi2 = ( RooChi2Var("chi2Fail","chi2Fail",pdfFail,dataFail,DataError(RooAbsData::Poisson)).getVal()
00292       +RooChi2Var("chi2Pass","chi2Pass",pdfPass,dataPass,DataError(RooAbsData::Poisson)).getVal() )/(2*pass->GetNbinsX()-8);
00293     if(chi2>3){
00294       efficiency.setVal(0.5);
00295       efficiency.setError(0.5);
00296     }
00297   }
00298 };
00299 
00300 //concrete fitter: voigtian signal plus exponential background
00301 class VoigtianPlusExponentialFitter: public AbstractFitter{
00302   protected:
00303   RooRealVar width;
00304   RooVoigtian voigtian;
00305   RooRealVar slopeFail;
00306   RooExponential exponentialFail;
00307   RooRealVar slopePass;
00308   RooExponential exponentialPass;
00309   RooAddPdf pdfFail;
00310   RooAddPdf pdfPass;
00311   public:
00312   VoigtianPlusExponentialFitter(bool verbose = false):
00313     AbstractFitter(verbose),
00314     width("width","width",2.5,"GeV"),
00315     voigtian("voigtian","voigtian",mass,mean,width,sigma),
00316     slopeFail("slopeFail","slopeFail",0.,-1.,0.),
00317     exponentialFail("linearFail","linearFail",mass,slopeFail),
00318     slopePass("slopePass","slopePass",0.,-1.,0.),
00319     exponentialPass("linearPass","linearPass",mass,slopePass),
00320     pdfFail("pdfFail","pdfFail", RooArgList(voigtian,exponentialFail), RooArgList(nSignalFail,nBackgroundFail)),
00321     pdfPass("pdfPass","pdfPass", RooArgList(voigtian,exponentialPass), RooArgList(nSignalPass,nBackgroundPass))
00322   {
00323     width.setConstant(kTRUE);
00324     simPdf.addPdf(pdfFail,"fail");
00325     simPdf.addPdf(pdfPass,"pass");
00326   };
00327   ~VoigtianPlusExponentialFitter(){};
00328   void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_, double width_){
00329     expectedMean = expectedMean_;
00330     expectedSigma = expectedSigma_;
00331     mass.setRange(massLow,massHigh);
00332     mean.setRange(massLow,massHigh);
00333     width.setVal(width_);
00334   }
00335   void fit(TH1* pass, TH1* all){
00336     using namespace RooFit;
00337     all->Add(pass,-1);
00338     TH1* &fail = all;
00339     if(!data) delete data;
00340     data = new RooDataHist("data", "data", mass, Index(category), Import("fail",*fail), Import("pass",*pass) );
00341     if(pass->Integral()+fail->Integral() < 5){
00342       efficiency.setVal(0.5);
00343       efficiency.setError(0.5);
00344       chi2 = 0;
00345       return;
00346     }
00347     mean.setVal(expectedMean);
00348     sigma.setVal(expectedSigma);
00349     efficiency.setVal(pass->Integral()/(pass->Integral()+fail->Integral()));
00350     nSignalAll.setVal(0.5*(fail->Integral()+pass->Integral()));
00351     nBackgroundFail.setVal(0.5*fail->Integral());
00352     nBackgroundPass.setVal(0.5*pass->Integral());
00353     slopeFail.setVal(0.);
00354     slopePass.setVal(0.);
00355     if(verbose){
00356       simPdf.fitTo( *data );
00357     }else{
00358       simPdf.fitTo( *data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1) );
00359     }
00360     RooDataHist dataFail("fail", "fail", mass, fail );
00361     RooDataHist dataPass("pass", "pass", mass, pass );
00362     chi2 = ( RooChi2Var("chi2Fail","chi2Fail",pdfFail,dataFail,DataError(RooAbsData::Poisson)).getVal()
00363       +RooChi2Var("chi2Pass","chi2Pass",pdfPass,dataPass,DataError(RooAbsData::Poisson)).getVal() )/(2*all->GetNbinsX()-8);
00364     if(chi2>3){
00365       efficiency.setVal(0.5);
00366       efficiency.setError(0.5);
00367     }
00368   }
00369 };
00370 
00371 }//namespace dqmTnP
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/DQMOffline/Trigger/plugins/GenericTnPFitter.h
