1 // TestFractionFitter
   2 #define TestFractionFitter_cxx
   3 #include "TestFractionFitter.h"
   4 
   5 using namespace std;
   6 
   7 ClassImp( TestFractionFitter );
   8 
   9 TestFractionFitter::~TestFractionFitter()
  10 {
  11   
  12 }
  13 
  14 TestFractionFitter::TestFractionFitter()
  15 {
  16   int mcfills = 100000;
  17   int datafills = 2000;
  18   double weight = 1.0;
  19   double relweight = 0.2;
  20   TH1F *h = new TH1F("h","test",20,-10,10);
  21   TH1F *hgaus = new TH1F("h1","hgaus",20,-10,10);
  22   TH1F *hlandau = new TH1F("h2","hlandau",20,-10,10);
  23   TF1 *landau = new TF1 ("mylandau","TMath::Landau(x,[0],[1],0)",-10,10);
  24   landau->SetParameters(0.0,2.0);
  25   hlandau->FillRandom("mylandau",mcfills);
  26   TF1 *gaus = new TF1 ("mygaus","TMath::Gaus(x,[0],[1],0)",-10,10);
  27   gaus->SetParameters(4.0,1.5);
  28   
  29 //   hgaus->FillRandom("mygaus",mcfills);
  30   
  31   for (int i=0;i<(int)(mcfills*relweight);++i) {
  32     double rand = gaus->GetRandom();
  33     hgaus->Fill(rand,weight); // i*weight?
  34   }
  35   
  36   h->Add(hlandau);
  37   h->Add(hgaus);
  38 //   h->FillRandom("mylandau",datafills);
  39 //   h->FillRandom("mygaus",(int)(datafills*relweight));
  40 //   h->Draw("");
  41   vector<TH1*> hmc;
  42   hmc.push_back(hlandau);
  43   hmc.push_back(hgaus);
  44   DoFit(h, hmc);
  45 
  46 }
  47 
  48 void TestFractionFitter::DoFit(TH1* hdata, vector<TH1*> hmc)
  49 {
  50   
  51   TObjArray *mc = new TObjArray(hmc.size());        // MC histograms are put in this array
  52   for (unsigned int i=0; i<hmc.size(); ++i) {
  53     mc->Add(hmc[i]);
  54   }
  55   TFractionFitter* fit = new TFractionFitter(hdata, mc); // initialise
  56 //   fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  57 //   fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
  58   Int_t status = fit->Fit();               // perform the fit
  59   cout << "fit status: " << status << endl;
  60   if (status == 0) {                       // check on fit status
  61     TH1F* result = (TH1F*) fit->GetPlot();
  62     hdata->Draw("Ep");
  63     result->Draw("same");
  64   }
  65   
  66 }