PamUnfold/src/PamUnfold.cpp

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TROOT.h>
#include <TList.h>
#include <TString.h>
//#include <TObjectString.h>
#include <TGraphAsymmErrors.h>
#include <TGraphErrors.h>
#include <TChain.h>
#include <TCutG.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TObjString.h>
#include <TMath.h>

#include <PamUnfold.h>

using namespace std;

ClassImp(PamUnfold);


PamUnfold::PamUnfold(TString name, TString title) : TNamed(name, title){

  cout << "WARNING:: entering PamUnfold::PamUnfold(TString name, TString title)" << endl;
  cout << "           empty constructor be sure to initialize measured and smearing" << endl;

  _measured = NULL;
  _smearing = NULL;
  _prior = NULL;

}

PamUnfold::~PamUnfold(){
}

PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing) : TNamed(name, title){

  _measured = measured;
  _smearing = smearing;
  _prior = NULL;

  if( !IsBinningOK() )
    cerr << " -- ERROR in PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing)" << endl;

  if( !IsSmNormalized() ){
    cerr << " -- WARNING in PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing)" << endl;
    cerr << " ---- Remember to provide the normalization histogram for the smearing matrix" << endl;
  }

  Init();

  _is_improved = kFALSE;
  _smooth = kFALSE;
  _smooth_opt = "ROOT";

  _nsamples = 500;
  _max_steps = 50;
  _min_chi2 = 1.;

}

void PamUnfold::AddExcludedBin(Int_t bin){
  _excluded_bins.push_back(bin);
}

TH1D* PamUnfold::GetMeasured(){
  return _measured;
}

TH2D* PamUnfold::GetSmearing(){
  return _smearing;
}

TH1D* PamUnfold::GetUnfolded(){
  TH1D* __unfolded = (TH1D*) _unfolded->Clone( Form("%s_%s_unf", _measured->GetName(), this->GetName()) );
  return __unfolded;
}

TList* PamUnfold::GetBinHistList(){
  return _bin_hist_list;
}

void PamUnfold::SetMeasured(TH1D* measured){
  _measured = measured;
}

void PamUnfold::SetSmearing(TH2D* smearing){

  _smearing = smearing;

  if( !IsBinningOK() )
    cerr << " -- ERROR in PamUnfold::SetSmearing(TH2D* smearing)" << endl;

  Init();

}

void PamUnfold::SetPrior(TH1D* prior){

  _prior = (TH1D*) prior->Clone( Form("_prior_%s", this->GetName()) );

  if( !IsBinningOK() )
    cerr << " -- ERROR in PamUnfold::SetPrior(TH1D* prior)" << endl;

}

void PamUnfold::SetNormalization(TH1D* norm){
  _norm = norm;

  //If the matrix is not normalized it is normalized now.
  if( !IsSmNormalized() ){
    cout << "      Normalizing smearing matrix" << endl;

    try{ if(!_norm) throw 1; }
    catch(int e){
      if(e==1){
        cerr << " -- ERROR in PamUnfold::Init()" << endl;
        cerr << " ---- Normalization histogram is NULL" << endl;
        cerr << "      execution is likely to crash." << endl;
      }
    }

    NormalizeMatrix();

  }

}

void PamUnfold::SetImproved(Bool_t improv){
  _is_improved = improv;
}

void PamUnfold::SetSmoothing(Bool_t smooth, TString smooth_opt){

  _smooth = smooth;
  _smooth_opt = smooth_opt;

  if(!_smooth_opt.CompareTo(""))
    _smooth_opt = "ROOT";

}

void PamUnfold::SetNsamples(UInt_t nsamples){
  _nsamples = nsamples;
}

void PamUnfold::SetMaxSteps(UInt_t max_steps){
  _max_steps = max_steps;
}

void PamUnfold::SetMinChi2(Double_t min_chi2){
  _min_chi2 = min_chi2;
}


Bool_t PamUnfold::IsSmNormalized(){

  Double_t sum = 0;
  for(UInt_t i=0; i<_smearing->GetNbinsX(); i++){
    sum = 0;
    for(UInt_t j=0; j<_smearing->GetNbinsY(); j++)
      sum += _smearing->GetBinContent(_smearing->GetBin(i+1,j+1));
    if(sum>1.0000001){
      cout << "Bin: " << i+1 << " sum=" << sum << endl;
      return kFALSE;
    }
  }

  return kTRUE;

}

void PamUnfold::WMovAvSmooth(TH1D* input, vector<Int_t>&excl){

  Double_t* xarr = new Double_t [input->GetNbinsX()];

  for(UInt_t ib=0; ib<input->GetNbinsX(); ib++)
    xarr[ib] = (Double_t) input->GetBinContent(ib+1);

  Bool_t excluding;
  for(UInt_t ib=1; ib<input->GetNbinsX()-1; ib++){

    excluding = kFALSE;
    for(UInt_t iex=0; iex<excl.size(); iex++)
      if(
         ib   == excl[iex] ||
         ib+1 == excl[iex] ||
         ib+2 == excl[iex]
         ) excluding = kTRUE;
    if(excluding) continue;

    Double_t xc = xarr[ib];
    Double_t xl = xarr[ib-1];
    Double_t xh = xarr[ib+1];
    Double_t x = 0.25*(xh+xl) + 0.5*xc;
    input->SetBinContent(ib+1, (Float_t) x);
  }

  delete xarr;

  return;

}

Double_t PamUnfold::GetChi2H( TH1D* h1, TH1D* h2 ){

  if(h1->GetNbinsX() != h2->GetNbinsX()){
    cout << " -- Warning in PamUnfold::GetChi2H(TH1D*, TH1D*) : Histograms have different number of bins" << endl;
    return -1;
  }

  const Double_t* x_1 = h1->GetXaxis()->GetXbins()->GetArray();
  const Double_t* x_2 = h2->GetXaxis()->GetXbins()->GetArray();

  for(Int_t ib=0; ib < h1->GetNbinsX()+1; ib++){
    if(x_1[ib] != x_2[ib]){
      cout << " -- Warning in PamUnfold::GetChi2H(TH1D*, TH1D*) : Histograms have different number of bins" << endl;
      return -1;
    }
    else
      continue;
  }


  Double_t chi2 = 0;
  for(UInt_t i=0; i<h1->GetNbinsX(); i++){
    if(h1->GetBinContent(i+1) + h2->GetBinContent(i+2) > 1)
      chi2 += pow(h1->GetBinContent(i+1) - h2->GetBinContent(i+1),2)/(h1->GetBinContent(i+1) + h2->GetBinContent(i+2));
    else
      chi2 += pow(h1->GetBinContent(i+1) - h2->GetBinContent(i+1),2);
  }

  return chi2;
}

void PamUnfold::Unfold(){


  //Smoothing is applied to the spectrum, not to the counts histogram!!!
  // ------------------------------------------------------
  //   Prior initialization
  // ------------------------------------------------------

  for(Int_t i=0; i<_prior->GetNbinsX(); i++) 
    _prior->SetBinContent(i+1, _prior->GetBinContent(i+1)/_prior->GetBinWidth(i+1) );

  if(_smooth){
    if(_smooth_opt.Contains("WMA")){
      cout << " --- Using Weighted Moving Average smoothing\n";
      WMovAvSmooth(_prior, _excluded_bins);
    }
    else if(_smooth_opt.Contains("ROOT")){
      cout << " --- Using standard ROOT smoothing\n";
      _prior->Smooth();
    }
    else
      cout << " --- WARNING: No valid smoothing option specified" << endl;
  }

  for(Int_t i=0; i<_prior->GetNbinsX(); i++) 
    _prior->SetBinContent(i+1, _prior->GetBinContent(i+1)*_prior->GetBinWidth(i+1) );

  _prior->Scale(1./_prior->GetSumOfWeights()); //The prior is a 'probability' so it has to be normalized at the very last step

  // ------------------------------------------------------
  //   Unfolding
  // ------------------------------------------------------

  cout << " --- UNFOLDING!!" << endl;

  TH2D* theta = (TH2D*) _smearing->Clone("theta");
  theta->Reset();

  for(Int_t i=0; i<theta->GetNbinsY(); i++){
    for(Int_t j=0; j<theta->GetNbinsX(); j++){

      Double_t s = _smearing->GetBinContent(_smearing->GetBin(i+1,j+1))*_prior->GetBinContent(i+1);
      Double_t ls = 0;
      for(Int_t k=0; k<_smearing->GetNbinsX(); k++)
        ls+=_smearing->GetBinContent(_smearing->GetBin(k+1,j+1))*_prior->GetBinContent(k+1);

      if(ls)
        theta->SetBinContent(theta->GetBin(j+1,i+1),s/ls);
    }
  }

  _unfolded->Reset();

  for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
    Double_t num = 0;
    Double_t den = 0;
    Double_t err = 0;
    for(Int_t j=0; j<theta->GetNbinsX(); j++){

      num+=theta->GetBinContent( theta->GetBin(j+1,i+1) )*_measured->GetBinContent(j+1);
      den+=_smearing->GetBinContent(_smearing->GetBin(i+1,j+1));
      err+=pow( theta->GetBinContent(theta->GetBin(j+1,i+1))*_measured->GetBinError(j+1), 1 );
      //      err+=pow( theta->GetBinContent(theta->GetBin(j+1,i+1))*sp->GetBinError(j+1), 2 );

    }

    if(den){
      _unfolded->SetBinContent(i+1, num/den);
      _unfolded->SetBinError(i+1, fabs(err)/den);
      //      _unfolded->SetBinError(i+1,sqrt(err)/den);
    }
  }

}


void PamUnfold::ImprovedUnfold(){


  vector<Int_t> mu(_measured->GetNbinsX());

  //Smoothing is applied to the spectrum, not to the counts histogram!!!
  // ------------------------------------------------------
  //   Prior initialization
  // ------------------------------------------------------

  for(Int_t i=0; i<_prior->GetNbinsX(); i++) 
    _prior->SetBinContent(i+1, _prior->GetBinContent(i+1)/_prior->GetBinWidth(i+1) );

  if(_smooth){
    if(_smooth_opt.Contains("WMA")){
      cout << " --- Using Weighted Moving Average smoothing\n";
      WMovAvSmooth(_prior, _excluded_bins);
    }
    else if(_smooth_opt.Contains("ROOT")){
      cout << " --- Using standard ROOT smoothing\n";
      _prior->Smooth();
    }
    else
      cout << " --- WARNING: No valid smoothing option specified" << endl;
  }

  for(Int_t i=0; i<_prior->GetNbinsX(); i++) 
    _prior->SetBinContent(i+1, _prior->GetBinContent(i+1)*_prior->GetBinWidth(i+1) );

  _prior->Scale(1./_prior->GetSumOfWeights()); //The prior is a 'probability' so it has to be normalized at the very last step

  // ------------------------------------------------------
  //   Unfolding
  // ------------------------------------------------------

  cout << " --- UNFOLDING!!" << endl;

  _bin_list->Clear();

  TMVA::Timer* timer = new TMVA::Timer(_nsamples, "unf");
  for(UInt_t sample=0; sample<_nsamples; sample++){
    
    cout << " ---- Sampling... " << timer->GetLeftTime(sample) << " left..." << endl << flush <<"\33[1A";

    SampleMatrix();
    
    TH2D* theta = (TH2D*) _smearing_sample->Clone("theta");
    theta->Reset();
    
    for(Int_t i=0; i<theta->GetNbinsY(); i++){
      for(Int_t j=0; j<theta->GetNbinsX(); j++){
        
        Double_t s = _smearing_sample->GetBinContent(_smearing_sample->GetBin(i+1,j+1))*_prior->GetBinContent(i+1);
        Double_t ls = 0;
        for(Int_t k=0; k<_smearing_sample->GetNbinsX(); k++)
          ls+=_smearing_sample->GetBinContent(_smearing_sample->GetBin(k+1,j+1))*_prior->GetBinContent(k+1);
        
        if(ls)
          theta->SetBinContent(theta->GetBin(j+1,i+1),s/ls);

      }
    }
    
    TVectorD* result = new TVectorD(_unfolded->GetNbinsX());
    
    //Sampling mu_j and rounding to nearest integer
    for(Int_t j=0; j<_measured->GetNbinsX(); j++){
      mu[j] = TMath::Nint( rangen->Gamma( 1 + _measured->GetBinContent(j+1), 1 ) );

      vector<Double_t> theta_vec_j(_unfolded->GetNbinsX());
      for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
        theta_vec_j[i] = theta->GetBinContent( theta->GetBin(j+1,i+1) );
        //      cout << theta_vec_j[i] << " ";
      }
      //      cout << endl;

      vector<Int_t> res_partial = rangen->Multinomial(mu[j], theta_vec_j);
      for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
        //      cout << res_partial[i] << " ";
        (*result)[i] += res_partial[i];
      }
      // cout << endl;
      // cout << _measured->GetBinContent(j+1) << " " << mu[j] << " " << res_partial[17] << " " << (*result)[17] << endl;

    }

    
    for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
      Double_t den = 0;
      for(Int_t j=0; j<_measured->GetNbinsX(); j++)
        den+=_smearing_sample->GetBinContent(_smearing_sample->GetBin(i+1,j+1));   

        (*result)[i] /= den;
    }

    _bin_list->Add(result);
    
  }

  BuildFlux();

}


Bool_t PamUnfold::IsBinningOK(){

  Int_t n_meas = _measured->GetNbinsX();
  Int_t n_smy = _smearing->GetNbinsY();


  if(n_meas != n_smy){
    cerr << " --- Different binning between measured and smearing" << endl;
    cerr << " ---- measured: " << n_meas << " bins;  smearing:" << n_smy << " bins on Y-axis" << endl;
    return kFALSE;
  }

  const Double_t* x_meas = _measured->GetXaxis()->GetXbins()->GetArray();
  const Double_t* y_sm = _smearing->GetYaxis()->GetXbins()->GetArray();

  for(Int_t ib=0; ib < n_meas+1; ib++){
    if( fabs(x_meas[ib] - y_sm[ib])/x_meas[ib] > 1e-4 ){
      cerr << " --- Different binning between measured and smearing" << endl;
      cerr << " ---- x_meas[" << ib << "] = " << x_meas[ib] << ";"
           <<"  y_sm[" << ib << "] = " << y_sm[ib] << ";" << endl;
      return kFALSE;
    }
    else
      continue;
  }


  if(_prior){

    Int_t n_prior = _prior->GetNbinsX();
    Int_t n_smx = _smearing->GetNbinsX();

    if(n_prior != n_smx){
      cerr << " --- Different binning between prior and smearing" << endl;
      return kFALSE;
    }

    const Double_t* x_prior = _prior->GetXaxis()->GetXbins()->GetArray();
    const Double_t* x_sm = _smearing->GetXaxis()->GetXbins()->GetArray();
    
    for(Int_t ib=0; ib < n_prior+1; ib++){
      if( fabs(x_prior[ib] - x_sm[ib])/x_sm[ib] > 1e-4 ){
        cerr << " --- Different binning between prior and smearing" << endl;
        return kFALSE;
      }
      else
        continue;
    }
  }


  return kTRUE;
  
}

void PamUnfold::IterativeUnfolding(UInt_t niter, TList* list){

  cout << " -- PamUnfold object " << this->GetName() << endl
       << "      Starting unfolding";

  if(!niter)
    cout << " with chi2 convergence check";
  else
    cout << " with " << niter << " iterations";

  if(_is_improved)
    cout << " and improved algorithm (may take a while...)";
  else
    cout << " and standard algorithm";

  cout << endl;

  if(list)
    cout << "      Saving iterations to TList object at " << list << endl;


  UInt_t iiter = 0;
  Double_t chi2 = 1e6;
  Bool_t kFlag = kTRUE;

  while( kFlag ){

    if(iiter > 0)
      SetPrior(_old_unfolded);

    if(_is_improved)
      ImprovedUnfold();
    else
      Unfold();

    if(iiter > 0){
      chi2 = GetChi2H( _unfolded, _old_unfolded );
      cout << "       Chi2 of change from iteration " << iiter-1 << " to " << iiter << " = " << chi2 << endl;
    }

    if(!niter)
      kFlag = chi2 < _min_chi2  ?  kFALSE : kTRUE;
    else
      kFlag = iiter == niter  ?  kFALSE : kTRUE;

    if(iiter >= _max_steps){
      kFlag = kFALSE;
      cerr << "WARNING: Unfolding procedure did not converge in less than " << _max_steps << " steps\n";
    }

    _old_unfolded = GetUnfolded();

    _old_unfolded->SetName( Form("%s_%03i", _old_unfolded->GetName(), iiter) );
    if(list){ 
      list->Add( _old_unfolded );
    }
    iiter++;

  }

  cout << "      Unfolding converged" << endl;

}

void PamUnfold::Init(){

  cout << " -- PamUnfold object " << this->GetName() << endl
       << "      Initializing unfolded histogram" << endl;

  _unfolded = new TH1D( Form("%s_%s_u", _measured->GetName(), this->GetName()), "", _smearing->GetNbinsX(), _smearing->GetXaxis()->GetXbins()->GetArray());
  _unfolded->Reset();

  _prior = (TH1D*) _unfolded->Clone( Form("_prior_%s",this->GetName()) );
  _prior->Reset();

  cout << "      Initializing flat Prior" << endl;
  //Initializing flat prior
  for(UInt_t ibin=0; ibin<_prior->GetNbinsX(); ibin++){
    _prior->SetBinContent(ibin+1, _prior->GetBinWidth(ibin+1)/(_prior->GetBinLowEdge(_prior->GetNbinsX()+1)-_prior->GetBinLowEdge(1)) );
  }

  cout << "      Initializing Random number generator" << endl;
  rangen = new RanGen();

  _bin_list = new TList();
  _bin_hist_list = new TList();
  _smearing_sample = (TH2D*) _smearing->Clone( Form("%s_sample", _smearing->GetName()) );

  cout << "      Initialization done" << endl;

}


void PamUnfold::NormalizeMatrix(){

  for(UInt_t ibiny=0; ibiny<_smearing->GetNbinsY(); ibiny++){
    for(UInt_t ibinx=0; ibinx<_smearing->GetNbinsX(); ibinx++){
      Int_t globalbin = _smearing->GetBin(ibinx+1, ibiny+1);
      if(_norm->GetBinContent(ibinx+1))
        _smearing->SetBinContent( globalbin, _smearing->GetBinContent(globalbin)/_norm->GetBinContent(ibinx+1) );
    }
  }

}

void PamUnfold::SampleMatrix(){

  _smearing_sample->Reset();

  vector<Int_t> alpha(_smearing->GetNbinsY()+1);
  vector<Double_t> sampled;

  for(UInt_t ibinx=0; ibinx<_smearing->GetNbinsX(); ibinx++){

    //    cout << _norm->GetBinContent(ibinx+1) << endl;
    Double_t snorm=0;
    for(UInt_t ibiny=0; ibiny<_smearing->GetNbinsY(); ibiny++){
      Int_t globalbin = _smearing->GetBin(ibinx+1, ibiny+1);
      snorm += _smearing->GetBinContent(globalbin);
      alpha[ibiny] = TMath::Nint( _smearing->GetBinContent(globalbin) * _norm->GetBinContent(ibinx+1) ); 
    }
    alpha[_smearing->GetNbinsY()] = TMath::Nint( (1 - snorm) * _norm->GetBinContent(ibinx+1) );
    //cout << "bin: " << ibinx << "  " << alpha[_smearing->GetNbinsY()] << endl;

    sampled = rangen->Dirichlet(_smearing->GetNbinsY() + 1, alpha);

    for(UInt_t ibiny=0; ibiny<_smearing_sample->GetNbinsY(); ibiny++){
      Int_t globalbin = _smearing_sample->GetBin(ibinx+1, ibiny+1);
      //      cout << sampled[ibiny] << " ";
      _smearing_sample->SetBinContent(globalbin, sampled[ibiny]);
    }
    //    cout << endl;
  }
  //  cout << endl;

}


void PamUnfold::BuildFlux(){

  _unfolded->Reset();
  _bin_hist_list->Delete();

  for(Int_t j=0; j<_unfolded->GetNbinsX(); j++){
    Double_t min=1e9;
    Double_t max=0;
    Double_t check=0;

    for(Int_t i=0; i<_bin_list->GetEntries(); i++){
      TVectorD* vv = (TVectorD*) _bin_list->At(i);
      check = (*vv)(j);
      if(check<min) min=check;
      if(check>max) max=check;
    }

    _bin_hist_list->Add( new TH1D( Form("hh_%03i", j), ";Events", 150, min-20, max+20 ) );
  }

  for(Int_t i=0; i<_bin_list->GetEntries(); i++){
    TVectorD* vv = (TVectorD*) _bin_list->At(i);
    //cout << "Bin " << i << " ";
    for(Int_t j=0; j<vv->GetNoElements(); j++){
      //cout << (*vv)(j) << " ";
      ((TH1D*) _bin_hist_list->At(j))->Fill( (*vv)(j) );
    }
    //    cout << endl;
  }

  Double_t qq = 0.682689492137;
  Double_t yq[2], xq[2];
  xq[0] = (1-qq)/2; xq[1] = xq[0] + qq;


  for(Int_t j=0; j<_unfolded->GetNbinsX(); j++){
    TH1D* mt = ((TH1D*) _bin_hist_list->At(j));
    mt->GetQuantiles(2, yq, xq);

    double mean = 0.5*(yq[0] + yq[1]);
    double err  = 0.5*(yq[1] - yq[0]);
    if(mt->ComputeIntegral()){
      _unfolded->SetBinContent(j+1, mean);
      _unfolded->SetBinError(j+1, err);
    }
    else{
      _unfolded->SetBinContent(j+1, 0);
      _unfolded->SetBinError(j+1, 0);
    }

  }

}

void PamUnfold::Draw(TString path){

  Int_t nx, ny;
  FindSplitting(_unfolded->GetNbinsX() , nx, ny);
  canv = dynamic_cast<TCanvas*>(gDirectory->FindObject( Form("canv_%s", this->GetName()) ));
  if(!canv) canv = new TCanvas(Form("canv_%s", this->GetName()), "Bin Distributions", 0, 0, 4 + 500*nx, 28 + 500*ny);
  canv->Divide(nx, ny);

  for(Int_t ip=0; ip<_unfolded->GetNbinsX(); ip++){
    canv->cd(ip+1);

    TH1D* mt = ((TH1D*) _bin_hist_list->At(ip));
    mt->Draw();

    Double_t mean = _unfolded->GetBinContent(ip+1);
    Double_t err  = _unfolded->GetBinError(ip+1);

    TLine* line1 = new TLine(mean-err, 0, mean-err, mt->GetMaximum());
    TLine* line2 = new TLine(mean+err, 0, mean+err, mt->GetMaximum());
    TLine* line3 = new TLine(mean, 0, mean, mt->GetMaximum());

    line1->SetLineWidth(2);
    line1->SetLineColor(2);
    line2->SetLineWidth(2);
    line2->SetLineColor(2);
    line3->SetLineWidth(2);
    line3->SetLineColor(4);

    line1->Draw("same");
    line2->Draw("same");
    line3->Draw("same");

  }

  if(path) canv->Print(path);

}

void PamUnfold::FindSplitting(Int_t n, Int_t &nx, Int_t &ny){

  Int_t x = TMath::Nint( sqrt(n) );
  if(x*x == n){ nx=ny=x; return; }

  x = floor( sqrt(n) ) + 1;
  for(Int_t y=0; y<=x; y++){
    if(y*x>=n){ ny = y; break;}
  }
  nx = x;

  if(ny>nx){ nx=ny; ny=x;}

}
1	mayorov	1.1	#include <iostream>
2			#include <fstream>
3			#include <stdio.h>
4			#include <math.h>
5			#include <TH1F.h>
6			#include <TH2F.h>
7			#include <TH1D.h>
8			#include <TH2D.h>
9			#include <TFile.h>
10			#include <TROOT.h>
11			#include <TList.h>
12			#include <TString.h>
13			//#include <TObjectString.h>
14			#include <TGraphAsymmErrors.h>
15			#include <TGraphErrors.h>
16			#include <TChain.h>
17			#include <TCutG.h>
18			#include <TF1.h>
19			#include <TCanvas.h>
20			#include <TObjString.h>
21			#include <TMath.h>
22
23			#include <PamUnfold.h>
24
25			using namespace std;
26
27			ClassImp(PamUnfold);
28
29
30			PamUnfold::PamUnfold(TString name, TString title) : TNamed(name, title){
31
32			cout << "WARNING:: entering PamUnfold::PamUnfold(TString name, TString title)" << endl;
33			cout << " empty constructor be sure to initialize measured and smearing" << endl;
34
35			_measured = NULL;
36			_smearing = NULL;
37			_prior = NULL;
38
39			}
40
41			PamUnfold::~PamUnfold(){
42			}
43
44			PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing) : TNamed(name, title){
45
46			_measured = measured;
47			_smearing = smearing;
48			_prior = NULL;
49
50			if( !IsBinningOK() )
51			cerr << " -- ERROR in PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing)" << endl;
52
53			if( !IsSmNormalized() ){
54			cerr << " -- WARNING in PamUnfold::PamUnfold(TString name, TString title, TH1D* measured, TH2D* smearing)" << endl;
55			cerr << " ---- Remember to provide the normalization histogram for the smearing matrix" << endl;
56			}
57
58			Init();
59
60			_is_improved = kFALSE;
61			_smooth = kFALSE;
62			_smooth_opt = "ROOT";
63
64			_nsamples = 500;
65			_max_steps = 50;
66			_min_chi2 = 1.;
67
68			}
69
70			void PamUnfold::AddExcludedBin(Int_t bin){
71			_excluded_bins.push_back(bin);
72			}
73
74			TH1D* PamUnfold::GetMeasured(){
75			return _measured;
76			}
77
78			TH2D* PamUnfold::GetSmearing(){
79			return _smearing;
80			}
81
82			TH1D* PamUnfold::GetUnfolded(){
83			TH1D* __unfolded = (TH1D*) _unfolded->Clone( Form("%s_%s_unf", _measured->GetName(), this->GetName()) );
84			return __unfolded;
85			}
86
87			TList* PamUnfold::GetBinHistList(){
88			return _bin_hist_list;
89			}
90
91			void PamUnfold::SetMeasured(TH1D* measured){
92			_measured = measured;
93			}
94
95			void PamUnfold::SetSmearing(TH2D* smearing){
96
97			_smearing = smearing;
98
99			if( !IsBinningOK() )
100			cerr << " -- ERROR in PamUnfold::SetSmearing(TH2D* smearing)" << endl;
101
102			Init();
103
104			}
105
106			void PamUnfold::SetPrior(TH1D* prior){
107
108			_prior = (TH1D*) prior->Clone( Form("_prior_%s", this->GetName()) );
109
110			if( !IsBinningOK() )
111			cerr << " -- ERROR in PamUnfold::SetPrior(TH1D* prior)" << endl;
112
113			}
114
115			void PamUnfold::SetNormalization(TH1D* norm){
116			_norm = norm;
117
118			//If the matrix is not normalized it is normalized now.
119			if( !IsSmNormalized() ){
120			cout << " Normalizing smearing matrix" << endl;
121
122			try{ if(!_norm) throw 1; }
123			catch(int e){
124			if(e==1){
125			cerr << " -- ERROR in PamUnfold::Init()" << endl;
126			cerr << " ---- Normalization histogram is NULL" << endl;
127			cerr << " execution is likely to crash." << endl;
128			}
129			}
130
131			NormalizeMatrix();
132
133			}
134
135			}
136
137			void PamUnfold::SetImproved(Bool_t improv){
138			_is_improved = improv;
139			}
140
141			void PamUnfold::SetSmoothing(Bool_t smooth, TString smooth_opt){
142
143			_smooth = smooth;
144			_smooth_opt = smooth_opt;
145
146			if(!_smooth_opt.CompareTo(""))
147			_smooth_opt = "ROOT";
148
149			}
150
151			void PamUnfold::SetNsamples(UInt_t nsamples){
152			_nsamples = nsamples;
153			}
154
155			void PamUnfold::SetMaxSteps(UInt_t max_steps){
156			_max_steps = max_steps;
157			}
158
159			void PamUnfold::SetMinChi2(Double_t min_chi2){
160			_min_chi2 = min_chi2;
161			}
162
163
164			Bool_t PamUnfold::IsSmNormalized(){
165
166			Double_t sum = 0;
167			for(UInt_t i=0; i<_smearing->GetNbinsX(); i++){
168			sum = 0;
169			for(UInt_t j=0; j<_smearing->GetNbinsY(); j++)
170			sum += _smearing->GetBinContent(_smearing->GetBin(i+1,j+1));
171			if(sum>1.0000001){
172			cout << "Bin: " << i+1 << " sum=" << sum << endl;
173			return kFALSE;
174			}
175			}
176
177			return kTRUE;
178
179			}
180
181			void PamUnfold::WMovAvSmooth(TH1D* input, vector<Int_t>&excl){
182
183			Double_t* xarr = new Double_t [input->GetNbinsX()];
184
185			for(UInt_t ib=0; ib<input->GetNbinsX(); ib++)
186			xarr[ib] = (Double_t) input->GetBinContent(ib+1);
187
188			Bool_t excluding;
189			for(UInt_t ib=1; ib<input->GetNbinsX()-1; ib++){
190
191			excluding = kFALSE;
192			for(UInt_t iex=0; iex<excl.size(); iex++)
193			if(
194			ib == excl[iex] \|\|
195			ib+1 == excl[iex] \|\|
196			ib+2 == excl[iex]
197			) excluding = kTRUE;
198			if(excluding) continue;
199
200			Double_t xc = xarr[ib];
201			Double_t xl = xarr[ib-1];
202			Double_t xh = xarr[ib+1];
203			Double_t x = 0.25(xh+xl) + 0.5xc;
204			input->SetBinContent(ib+1, (Float_t) x);
205			}
206
207			delete xarr;
208
209			return;
210
211			}
212
213			Double_t PamUnfold::GetChi2H( TH1D* h1, TH1D* h2 ){
214
215			if(h1->GetNbinsX() != h2->GetNbinsX()){
216			cout << " -- Warning in PamUnfold::GetChi2H(TH1D, TH1D) : Histograms have different number of bins" << endl;
217			return -1;
218			}
219
220			const Double_t* x_1 = h1->GetXaxis()->GetXbins()->GetArray();
221			const Double_t* x_2 = h2->GetXaxis()->GetXbins()->GetArray();
222
223			for(Int_t ib=0; ib < h1->GetNbinsX()+1; ib++){
224			if(x_1[ib] != x_2[ib]){
225			cout << " -- Warning in PamUnfold::GetChi2H(TH1D, TH1D) : Histograms have different number of bins" << endl;
226			return -1;
227			}
228			else
229			continue;
230			}
231
232
233
234			Double_t chi2 = 0;
235			for(UInt_t i=0; i<h1->GetNbinsX(); i++){
236			if(h1->GetBinContent(i+1) + h2->GetBinContent(i+2) > 1)
237			chi2 += pow(h1->GetBinContent(i+1) - h2->GetBinContent(i+1),2)/(h1->GetBinContent(i+1) + h2->GetBinContent(i+2));
238			else
239			chi2 += pow(h1->GetBinContent(i+1) - h2->GetBinContent(i+1),2);
240			}
241
242			return chi2;
243			}
244
245			void PamUnfold::Unfold(){
246
247
248			//Smoothing is applied to the spectrum, not to the counts histogram!!!
249			// ------------------------------------------------------
250			// Prior initialization
251			// ------------------------------------------------------
252
253			for(Int_t i=0; i<_prior->GetNbinsX(); i++)
254			_prior->SetBinContent(i+1, _prior->GetBinContent(i+1)/_prior->GetBinWidth(i+1) );
255
256			if(_smooth){
257			if(_smooth_opt.Contains("WMA")){
258			cout << " --- Using Weighted Moving Average smoothing\n";
259			WMovAvSmooth(_prior, _excluded_bins);
260			}
261			else if(_smooth_opt.Contains("ROOT")){
262			cout << " --- Using standard ROOT smoothing\n";
263			_prior->Smooth();
264			}
265			else
266			cout << " --- WARNING: No valid smoothing option specified" << endl;
267			}
268
269			for(Int_t i=0; i<_prior->GetNbinsX(); i++)
270			_prior->SetBinContent(i+1, _prior->GetBinContent(i+1)*_prior->GetBinWidth(i+1) );
271
272			_prior->Scale(1./_prior->GetSumOfWeights()); //The prior is a 'probability' so it has to be normalized at the very last step
273
274			// ------------------------------------------------------
275			// Unfolding
276			// ------------------------------------------------------
277
278			cout << " --- UNFOLDING!!" << endl;
279
280			TH2D* theta = (TH2D*) _smearing->Clone("theta");
281			theta->Reset();
282
283			for(Int_t i=0; i<theta->GetNbinsY(); i++){
284			for(Int_t j=0; j<theta->GetNbinsX(); j++){
285
286			Double_t s = _smearing->GetBinContent(_smearing->GetBin(i+1,j+1))*_prior->GetBinContent(i+1);
287			Double_t ls = 0;
288			for(Int_t k=0; k<_smearing->GetNbinsX(); k++)
289			ls+=_smearing->GetBinContent(_smearing->GetBin(k+1,j+1))*_prior->GetBinContent(k+1);
290
291			if(ls)
292			theta->SetBinContent(theta->GetBin(j+1,i+1),s/ls);
293			}
294			}
295
296			_unfolded->Reset();
297
298			for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
299			Double_t num = 0;
300			Double_t den = 0;
301			Double_t err = 0;
302			for(Int_t j=0; j<theta->GetNbinsX(); j++){
303
304			num+=theta->GetBinContent( theta->GetBin(j+1,i+1) )*_measured->GetBinContent(j+1);
305			den+=_smearing->GetBinContent(_smearing->GetBin(i+1,j+1));
306			err+=pow( theta->GetBinContent(theta->GetBin(j+1,i+1))*_measured->GetBinError(j+1), 1 );
307			// err+=pow( theta->GetBinContent(theta->GetBin(j+1,i+1))*sp->GetBinError(j+1), 2 );
308
309			}
310
311			if(den){
312			_unfolded->SetBinContent(i+1, num/den);
313			_unfolded->SetBinError(i+1, fabs(err)/den);
314			// _unfolded->SetBinError(i+1,sqrt(err)/den);
315			}
316			}
317
318			}
319
320
321			void PamUnfold::ImprovedUnfold(){
322
323
324			vector<Int_t> mu(_measured->GetNbinsX());
325
326			//Smoothing is applied to the spectrum, not to the counts histogram!!!
327			// ------------------------------------------------------
328			// Prior initialization
329			// ------------------------------------------------------
330
331			for(Int_t i=0; i<_prior->GetNbinsX(); i++)
332			_prior->SetBinContent(i+1, _prior->GetBinContent(i+1)/_prior->GetBinWidth(i+1) );
333
334			if(_smooth){
335			if(_smooth_opt.Contains("WMA")){
336			cout << " --- Using Weighted Moving Average smoothing\n";
337			WMovAvSmooth(_prior, _excluded_bins);
338			}
339			else if(_smooth_opt.Contains("ROOT")){
340			cout << " --- Using standard ROOT smoothing\n";
341			_prior->Smooth();
342			}
343			else
344			cout << " --- WARNING: No valid smoothing option specified" << endl;
345			}
346
347			for(Int_t i=0; i<_prior->GetNbinsX(); i++)
348			_prior->SetBinContent(i+1, _prior->GetBinContent(i+1)*_prior->GetBinWidth(i+1) );
349
350			_prior->Scale(1./_prior->GetSumOfWeights()); //The prior is a 'probability' so it has to be normalized at the very last step
351
352			// ------------------------------------------------------
353			// Unfolding
354			// ------------------------------------------------------
355
356			cout << " --- UNFOLDING!!" << endl;
357
358			_bin_list->Clear();
359
360			TMVA::Timer* timer = new TMVA::Timer(_nsamples, "unf");
361			for(UInt_t sample=0; sample<_nsamples; sample++){
362
363			cout << " ---- Sampling... " << timer->GetLeftTime(sample) << " left..." << endl << flush <<"\33[1A";
364
365			SampleMatrix();
366
367			TH2D* theta = (TH2D*) _smearing_sample->Clone("theta");
368			theta->Reset();
369
370			for(Int_t i=0; i<theta->GetNbinsY(); i++){
371			for(Int_t j=0; j<theta->GetNbinsX(); j++){
372
373			Double_t s = _smearing_sample->GetBinContent(_smearing_sample->GetBin(i+1,j+1))*_prior->GetBinContent(i+1);
374			Double_t ls = 0;
375			for(Int_t k=0; k<_smearing_sample->GetNbinsX(); k++)
376			ls+=_smearing_sample->GetBinContent(_smearing_sample->GetBin(k+1,j+1))*_prior->GetBinContent(k+1);
377
378			if(ls)
379			theta->SetBinContent(theta->GetBin(j+1,i+1),s/ls);
380
381			}
382			}
383
384			TVectorD* result = new TVectorD(_unfolded->GetNbinsX());
385
386			//Sampling mu_j and rounding to nearest integer
387			for(Int_t j=0; j<_measured->GetNbinsX(); j++){
388			mu[j] = TMath::Nint( rangen->Gamma( 1 + _measured->GetBinContent(j+1), 1 ) );
389
390			vector<Double_t> theta_vec_j(_unfolded->GetNbinsX());
391			for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
392			theta_vec_j[i] = theta->GetBinContent( theta->GetBin(j+1,i+1) );
393			// cout << theta_vec_j[i] << " ";
394			}
395			// cout << endl;
396
397			vector<Int_t> res_partial = rangen->Multinomial(mu[j], theta_vec_j);
398			for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
399			// cout << res_partial[i] << " ";
400			(*result)[i] += res_partial[i];
401			}
402			// cout << endl;
403			// cout << _measured->GetBinContent(j+1) << " " << mu[j] << " " << res_partial[17] << " " << (*result)[17] << endl;
404
405			}
406
407
408			for(Int_t i=0; i<_unfolded->GetNbinsX(); i++){
409			Double_t den = 0;
410			for(Int_t j=0; j<_measured->GetNbinsX(); j++)
411			den+=_smearing_sample->GetBinContent(_smearing_sample->GetBin(i+1,j+1));
412
413			(*result)[i] /= den;
414			}
415
416			_bin_list->Add(result);
417
418			}
419
420			BuildFlux();
421
422			}
423
424
425			Bool_t PamUnfold::IsBinningOK(){
426
427			Int_t n_meas = _measured->GetNbinsX();
428			Int_t n_smy = _smearing->GetNbinsY();
429
430
431			if(n_meas != n_smy){
432			cerr << " --- Different binning between measured and smearing" << endl;
433			cerr << " ---- measured: " << n_meas << " bins; smearing:" << n_smy << " bins on Y-axis" << endl;
434			return kFALSE;
435			}
436
437			const Double_t* x_meas = _measured->GetXaxis()->GetXbins()->GetArray();
438			const Double_t* y_sm = _smearing->GetYaxis()->GetXbins()->GetArray();
439
440			for(Int_t ib=0; ib < n_meas+1; ib++){
441			if( fabs(x_meas[ib] - y_sm[ib])/x_meas[ib] > 1e-4 ){
442			cerr << " --- Different binning between measured and smearing" << endl;
443			cerr << " ---- x_meas[" << ib << "] = " << x_meas[ib] << ";"
444			<<" y_sm[" << ib << "] = " << y_sm[ib] << ";" << endl;
445			return kFALSE;
446			}
447			else
448			continue;
449			}
450
451
452
453			if(_prior){
454
455			Int_t n_prior = _prior->GetNbinsX();
456			Int_t n_smx = _smearing->GetNbinsX();
457
458			if(n_prior != n_smx){
459			cerr << " --- Different binning between prior and smearing" << endl;
460			return kFALSE;
461			}
462
463			const Double_t* x_prior = _prior->GetXaxis()->GetXbins()->GetArray();
464			const Double_t* x_sm = _smearing->GetXaxis()->GetXbins()->GetArray();
465
466			for(Int_t ib=0; ib < n_prior+1; ib++){
467			if( fabs(x_prior[ib] - x_sm[ib])/x_sm[ib] > 1e-4 ){
468			cerr << " --- Different binning between prior and smearing" << endl;
469			return kFALSE;
470			}
471			else
472			continue;
473			}
474			}
475
476
477			return kTRUE;
478
479			}
480
481			void PamUnfold::IterativeUnfolding(UInt_t niter, TList* list){
482
483			cout << " -- PamUnfold object " << this->GetName() << endl
484			<< " Starting unfolding";
485
486			if(!niter)
487			cout << " with chi2 convergence check";
488			else
489			cout << " with " << niter << " iterations";
490
491			if(_is_improved)
492			cout << " and improved algorithm (may take a while...)";
493			else
494			cout << " and standard algorithm";
495
496			cout << endl;
497
498			if(list)
499			cout << " Saving iterations to TList object at " << list << endl;
500
501
502			UInt_t iiter = 0;
503			Double_t chi2 = 1e6;
504			Bool_t kFlag = kTRUE;
505
506			while( kFlag ){
507
508			if(iiter > 0)
509			SetPrior(_old_unfolded);
510
511			if(_is_improved)
512			ImprovedUnfold();
513			else
514			Unfold();
515
516			if(iiter > 0){
517			chi2 = GetChi2H( _unfolded, _old_unfolded );
518			cout << " Chi2 of change from iteration " << iiter-1 << " to " << iiter << " = " << chi2 << endl;
519			}
520
521			if(!niter)
522			kFlag = chi2 < _min_chi2 ? kFALSE : kTRUE;
523			else
524			kFlag = iiter == niter ? kFALSE : kTRUE;
525
526			if(iiter >= _max_steps){
527			kFlag = kFALSE;
528			cerr << "WARNING: Unfolding procedure did not converge in less than " << _max_steps << " steps\n";
529			}
530
531			_old_unfolded = GetUnfolded();
532
533			_old_unfolded->SetName( Form("%s_%03i", _old_unfolded->GetName(), iiter) );
534			if(list){
535			list->Add( _old_unfolded );
536			}
537			iiter++;
538
539			}
540
541			cout << " Unfolding converged" << endl;
542
543			}
544
545			void PamUnfold::Init(){
546
547			cout << " -- PamUnfold object " << this->GetName() << endl
548			<< " Initializing unfolded histogram" << endl;
549
550			_unfolded = new TH1D( Form("%s_%s_u", _measured->GetName(), this->GetName()), "", _smearing->GetNbinsX(), _smearing->GetXaxis()->GetXbins()->GetArray());
551			_unfolded->Reset();
552
553			_prior = (TH1D*) _unfolded->Clone( Form("_prior_%s",this->GetName()) );
554			_prior->Reset();
555
556			cout << " Initializing flat Prior" << endl;
557			//Initializing flat prior
558			for(UInt_t ibin=0; ibin<_prior->GetNbinsX(); ibin++){
559			_prior->SetBinContent(ibin+1, _prior->GetBinWidth(ibin+1)/(_prior->GetBinLowEdge(_prior->GetNbinsX()+1)-_prior->GetBinLowEdge(1)) );
560			}
561
562			cout << " Initializing Random number generator" << endl;
563			rangen = new RanGen();
564
565			_bin_list = new TList();
566			_bin_hist_list = new TList();
567			_smearing_sample = (TH2D*) _smearing->Clone( Form("%s_sample", _smearing->GetName()) );
568
569			cout << " Initialization done" << endl;
570
571			}
572
573
574			void PamUnfold::NormalizeMatrix(){
575
576			for(UInt_t ibiny=0; ibiny<_smearing->GetNbinsY(); ibiny++){
577			for(UInt_t ibinx=0; ibinx<_smearing->GetNbinsX(); ibinx++){
578			Int_t globalbin = _smearing->GetBin(ibinx+1, ibiny+1);
579			if(_norm->GetBinContent(ibinx+1))
580			_smearing->SetBinContent( globalbin, _smearing->GetBinContent(globalbin)/_norm->GetBinContent(ibinx+1) );
581			}
582			}
583
584			}
585
586			void PamUnfold::SampleMatrix(){
587
588			_smearing_sample->Reset();
589
590			vector<Int_t> alpha(_smearing->GetNbinsY()+1);
591			vector<Double_t> sampled;
592
593			for(UInt_t ibinx=0; ibinx<_smearing->GetNbinsX(); ibinx++){
594
595			// cout << _norm->GetBinContent(ibinx+1) << endl;
596			Double_t snorm=0;
597			for(UInt_t ibiny=0; ibiny<_smearing->GetNbinsY(); ibiny++){
598			Int_t globalbin = _smearing->GetBin(ibinx+1, ibiny+1);
599			snorm += _smearing->GetBinContent(globalbin);
600			alpha[ibiny] = TMath::Nint( _smearing->GetBinContent(globalbin) * _norm->GetBinContent(ibinx+1) );
601			}
602			alpha[_smearing->GetNbinsY()] = TMath::Nint( (1 - snorm) * _norm->GetBinContent(ibinx+1) );
603			//cout << "bin: " << ibinx << " " << alpha[_smearing->GetNbinsY()] << endl;
604
605			sampled = rangen->Dirichlet(_smearing->GetNbinsY() + 1, alpha);
606
607			for(UInt_t ibiny=0; ibiny<_smearing_sample->GetNbinsY(); ibiny++){
608			Int_t globalbin = _smearing_sample->GetBin(ibinx+1, ibiny+1);
609			// cout << sampled[ibiny] << " ";
610			_smearing_sample->SetBinContent(globalbin, sampled[ibiny]);
611			}
612			// cout << endl;
613			}
614			// cout << endl;
615
616			}
617
618
619			void PamUnfold::BuildFlux(){
620
621			_unfolded->Reset();
622			_bin_hist_list->Delete();
623
624			for(Int_t j=0; j<_unfolded->GetNbinsX(); j++){
625			Double_t min=1e9;
626			Double_t max=0;
627			Double_t check=0;
628
629			for(Int_t i=0; i<_bin_list->GetEntries(); i++){
630			TVectorD* vv = (TVectorD*) _bin_list->At(i);
631			check = (*vv)(j);
632			if(check<min) min=check;
633			if(check>max) max=check;
634			}
635
636			_bin_hist_list->Add( new TH1D( Form("hh_%03i", j), ";Events", 150, min-20, max+20 ) );
637			}
638
639			for(Int_t i=0; i<_bin_list->GetEntries(); i++){
640			TVectorD* vv = (TVectorD*) _bin_list->At(i);
641			//cout << "Bin " << i << " ";
642			for(Int_t j=0; j<vv->GetNoElements(); j++){
643			//cout << (*vv)(j) << " ";
644			((TH1D) _bin_hist_list->At(j))->Fill( (vv)(j) );
645			}
646			// cout << endl;
647			}
648
649			Double_t qq = 0.682689492137;
650			Double_t yq[2], xq[2];
651			xq[0] = (1-qq)/2; xq[1] = xq[0] + qq;
652
653
654			for(Int_t j=0; j<_unfolded->GetNbinsX(); j++){
655			TH1D* mt = ((TH1D*) _bin_hist_list->At(j));
656			mt->GetQuantiles(2, yq, xq);
657
658			double mean = 0.5*(yq[0] + yq[1]);
659			double err = 0.5*(yq[1] - yq[0]);
660			if(mt->ComputeIntegral()){
661			_unfolded->SetBinContent(j+1, mean);
662			_unfolded->SetBinError(j+1, err);
663			}
664			else{
665			_unfolded->SetBinContent(j+1, 0);
666			_unfolded->SetBinError(j+1, 0);
667			}
668
669			}
670
671			}
672
673			void PamUnfold::Draw(TString path){
674
675			Int_t nx, ny;
676			FindSplitting(_unfolded->GetNbinsX() , nx, ny);
677			canv = dynamic_cast<TCanvas*>(gDirectory->FindObject( Form("canv_%s", this->GetName()) ));
678			if(!canv) canv = new TCanvas(Form("canv_%s", this->GetName()), "Bin Distributions", 0, 0, 4 + 500nx, 28 + 500ny);
679			canv->Divide(nx, ny);
680
681			for(Int_t ip=0; ip<_unfolded->GetNbinsX(); ip++){
682			canv->cd(ip+1);
683
684			TH1D* mt = ((TH1D*) _bin_hist_list->At(ip));
685			mt->Draw();
686
687			Double_t mean = _unfolded->GetBinContent(ip+1);
688			Double_t err = _unfolded->GetBinError(ip+1);
689
690			TLine* line1 = new TLine(mean-err, 0, mean-err, mt->GetMaximum());
691			TLine* line2 = new TLine(mean+err, 0, mean+err, mt->GetMaximum());
692			TLine* line3 = new TLine(mean, 0, mean, mt->GetMaximum());
693
694			line1->SetLineWidth(2);
695			line1->SetLineColor(2);
696			line2->SetLineWidth(2);
697			line2->SetLineColor(2);
698			line3->SetLineWidth(2);
699			line3->SetLineColor(4);
700
701			line1->Draw("same");
702			line2->Draw("same");
703			line3->Draw("same");
704
705			}
706
707			if(path) canv->Print(path);
708
709			}
710
711			void PamUnfold::FindSplitting(Int_t n, Int_t &nx, Int_t &ny){
712
713			Int_t x = TMath::Nint( sqrt(n) );
714			if(x*x == n){ nx=ny=x; return; }
715
716			x = floor( sqrt(n) ) + 1;
717			for(Int_t y=0; y<=x; y++){
718			if(y*x>=n){ ny = y; break;}
719			}
720			nx = x;
721
722			if(ny>nx){ nx=ny; ny=x;}
723
724			}