//unfolding exmple for my analysis
//last modified nov 17, 2008

//decide whether to unfold mx(1) or mx2(2) und the regularization param
//(use values like 2,3,4,5 for playing, the largest possible k is the number
//of bins in the spectrum, i.e. 16 for mx and 18 for mx2)
//k=1 is maximal regularization (you get your MC input back) and k=max
//is no regularization
void unfold(int specdim, int k){

  //for the filename
  char basedir[300];
  char namebase[50];
  char var[10];
  char fspec[20];
  char fullname[300];

  sprintf(basedir, "/data/babar/kerstin/vubunf/Ibu/R14tests/Ibu");
  sprintf(fspec, "Fitres.root");
  if(specdim==1){
    sprintf(var, "mx");
  }else if(specdim==2){
    sprintf(var, "mx2");
  }else{
    cout << "Do not know about mx^" << specdim << " spectrum..." << endl;
    break;
  }

  //for the histo names
  char hgen[20];
  char hreco[20];
  char hdata[20];
  char hdetmat[20];
  sprintf(hgen, "mxhadgenwoph_20");
  sprintf(hreco, "mxhadfit_11");
  sprintf(hdata, "subdatachop");
  sprintf(hdetmat, "detmat_4000");

  //*the* spectrum (i.e. for the mean)
  sprintf(namebase, "ujun07spec");
  sprintf(fullname, "%s%s%s/%s%s%s", basedir, namebase, var, namebase, var, fspec);
  cout << "Reading from " << fullname << endl; 
  TFile *fdata = new TFile(fullname);
  //the true MC distribution, note that this needs to be completely unbiased, i.e.
  //you cannot take the truth after skimming or any cuts!
  TH1D *xdata = fdata->Get(hgen);

  //the reco'd MC spectrum
  TH1D *bdata = fdata->Get(hreco);

  //the detector response matrix, where the entries are numbers of events after
  //all cuts etc., make sure to have x and y correct!
  TH2D *Adata = fdata->Get(hdetmat);

  //the data spectrum
  TH1D *sdata = fdata->Get(hdata);

  //get the parameters for the histos
  double min = sdata->GetXaxis()->GetXmin();
  double max = sdata->GetXaxis()->GetXmax();
  int numbins = sdata->GetNbinsX();
  cout << numbins << " bins, min " << min << " and max " << max << endl;


  //the covariance matrices on the measured spectrum, stat only, 
  //i.e. no correlations in my case
  TMatrixD sstatcov(numbins,numbins);

  for(int i=0; i<numbins; i++){
    sstatcov(i,i) = sdata->GetBinError(i+1)*sdata->GetBinError(i+1);
    cout << sdata->GetBinContent(i+1) << " " << TMath::Sqrt(sdata->GetBinContent(i+1)) << "  " << sdata->GetBinError(i+1) << endl;
  }

  //the vectors and histos we need
  TVectorD unfres(numbins);
  TH1D* udata = new TH1D("udata","",numbins,min,max);

  //the matrices we need
  TMatrixD ustatcov(numbins,numbins), ucov(numbins,numbins);
  
  cout << "**************** NOW DO THE UNFOLDING ****************" << endl;
  TUnfHisto *unfhist = new TUnfHisto(numbins,numbins);
  unfhist->init(sdata, bdata, xdata, Adata);
  unfres = unfhist->Unfold(k);

  //  cout << "The unfolded spectrum" << endl;
  for(Int_t i=0; i<numbins; i++){
    udata->SetBinContent(i+1,unfres(i));
    //    cout << unfres(i) << endl;
  }

  //  now the statistical covmat
   ustatcov=unfhist->GetCov(sstatcov, sdata, 1000, k);

   //  the full covarianve matrix
   ucov = ustatcov;

   //  put the error bars on the unfolded spectrum
  for(int i=0; i<numbins; i++){
    udata->SetBinError(i+1, TMath::Sqrt(ucov(i,i)));
  }

  //  cosmetics
  udata->SetMarkerStyle(20);
  udata->SetLineWidth(2);
  if(specdim==1){
    udata->SetXTitle("M_{X} / GeV");}
  else if(specdim==2){
    udata->SetXTitle("M_X^{2} / GeV^{2}");
  }

  TCanvas *canv1 = new TCanvas("canv1","canv1",500,400);
  udata->Draw();


}