csander/ToyMC/toyMC_LVMini.cpp

#include "TROOT.h"
#include "TSystem.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TStyle.h"
#include "TObjArray.h"
#include "TGraph.h"
#include "TMath.h"
#include "TSpline.h"

#include <iostream>
#include <vector>
#include <string>

#include "external.h"

using namespace::std;

////////////////////////////////////////////////////////////
class TDataPoint{

private:
  double weight;
  double weight_orig;
  double truth;
  double meas;
  double error;

public:
  double GetTruth(){return truth;};
  double GetMeas(){return meas;};
  double GetError(){return error;};
  double GetWeight(){return weight;};
  double GetWeight_orig(){return weight_orig;};
  void SetWeight(double w){weight=w;};
  void SetWeight_orig(double w){weight_orig=w;};

  //constructor
  TDataPoint(double t, double m, double e, double w){
    weight = w;
    weight_orig = w;
    truth = t;
    meas = m;
    error = e;
  };

  //destructor
  ~TDataPoint(){};

};
////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////
class TData{

private:
  int Nevts;
  double minX, maxX;
  double minY, maxY;
  double dY;
  vector<TDataPoint> data;
  vector<double> sp;

  //----------------------------------------------------------
  double chi2(vector<double>* PAR){
    
    TGraph cgr(sp.size(), &sp.front(), &PAR->front());
    TSpline3 *cspl = new TSpline3("cspl",&cgr);
    double x2 = 0.;
    for (vector<TDataPoint>::iterator i = data.begin(); i < data.end(); ++i){
      double mess=i->GetMeas();
      double truth=i->GetTruth();
      /*
      double c  = calib(mess,cspl);
      double dc = dcalib(mess,cspl);
      double d  = c  * mess - truth;
      double e  = dc * i->GetError();
      */
      /*
      double c  = calib(truth,cspl);
      double d  = c  * mess - truth;
      double e  = c * i->GetError();
      */
      double cmess = calib(mess,cspl);
      double cold = 1;
      double cinv  = calibinv(truth,cmess,cold,cspl);
      double d  = mess - truth/cinv;
      double e  = i->GetError();
      x2 += i->GetWeight()*pow(d/e, 2);
    }
    return x2;
  }
  //----------------------------------------------------------

  //----------------------------------------------------------
  double calib(double &x, TSpline3* calibSpline){
    return calibSpline->Eval(x);
  }
  //----------------------------------------------------------
  
  //----------------------------------------------------------
  //might be usefull if calibration constant is defined as function of truth
  double calibinv(double& x, double &cmess, double &cold, TSpline3* calibSpline){
    double cinv = calibSpline->Eval(x/cmess);
    //cout<<"cold: "<<cold<<" cnew: "<<cinv<<endl;
    double delta = fabs((cinv-cold)/cold);
    if (delta<1.e-4){
      return cinv;
    } else {
      return calibSpline->Eval(x/calibinv(x,cmess,cinv,calibSpline));
    }
  }
  //----------------------------------------------------------
  
  //----------------------------------------------------------
  double dcalib(double &mess, TSpline3* calibSpline){
    double eps=1.;
    double dc=((mess+eps)*calibSpline->Eval(mess+eps)-
               (mess-eps)*calibSpline->Eval(mess-eps))/2./eps;
    return dc;
  }
  //----------------------------------------------------------
  
  //----------------------------------------------------------
  //error parametrization, width of measured distribution
  double error(double x){
    double error = -1;
    if ( x >= 0 ) error = 1.25*sqrt(x);
    return error;
  }
  //----------------------------------------------------------

  //----------------------------------------------------------
  //reweight events to have flat truth distribution
  void FlattenSpectra(){

    TH1F* hist1 = new TH1F("truth non weighted","True distribution of measurement", 50, 0, 1000);
    TH1F* hist2 = new TH1F("truth weighted","Weighted true distribution of measurement", 50, 0, 1000);
    TH1F* weights = new TH1F("weights","Weights to flatten distribution", 50, 0, 1000);
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      hist1->Fill(i->GetTruth(),i->GetWeight());
    }
    double Integral1=hist1->Integral();
    for (int i=0; i<hist1->GetNbinsX(); ++i){
      weights->SetBinContent(i,1./hist1->GetBinContent(i));
    }
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      double w=weights->GetBinContent(hist1->FindBin(i->GetTruth()));
      hist2->Fill(i->GetTruth(),w);
    }
    double Integral2=hist2->Integral();
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      double w=weights->GetBinContent(hist1->FindBin(i->GetTruth()));
      w*=i->GetWeight();
      w*=Integral1/Integral2;
      i->SetWeight(w);
      i->SetWeight_orig(w);
    }
  }
  //----------------------------------------------------------

  //----------------------------------------------------------
  //reweight spectra to correct for energy dependent resolution
  void ResolutionReweight(bool logscale = false, bool outlier=false, int sigmamax=3){

    cout<<"Events at start of ResolutionReweight: "<<data.size()<<endl;

    int Nbins=49;
    int minEntries=200;
    vector<TDataPoint> dataCut;
    vector<double> mean;
    vector<double> sigma;

    double d;
    if (logscale){
      d=(log(maxY)-log(minY))/Nbins;
    } else {
      d=(maxY-minY)/Nbins;
    }

    //book histos
    TH1F* hist[Nbins];
    TF1* f[Nbins];
    for (int n=0; n<Nbins; ++n){
      char hname[100];
      sprintf(hname,"hr%i",n);
      hist[n] = new TH1F(hname,"Truth for MeasBin", 1000, minX, maxX);
      hist[n]->Sumw2();
    }

    //Fill histos
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      int n;
      if (logscale){
        n=(int)(Nbins*(log(i->GetMeas())-log(minY))/(log(maxY)-log(minY)));
        if (n>=Nbins) n=Nbins-1;
      } else {
        n=(int)(Nbins*(i->GetMeas()-minY)/d);
        if (n>=Nbins) n=Nbins-1;
      };
      hist[n]->Fill(i->GetTruth(),i->GetWeight());
    }

    //Fit histos
    for (int n=0; n<Nbins; ++n){
      cout<<n<<" of "<<Nbins<<endl;
      if (hist[n]->GetEntries()<minEntries){
        mean.push_back(0);
        sigma.push_back(0);
        continue;
      }
      char fname[100];
      sprintf(fname,"fr%i",n);
      f[n] = new TF1(fname,"gausn",minX,maxX);
      //use a guess for starting values
      f[n]->SetParameters(hist[n]->Integral("width"),hist[n]->GetMean(),hist[n]->GetRMS());
      hist[n]->Fit(fname,"ILLRQN");
      double Mean=f[n]->GetParameter(1);
      double Sigma=f[n]->GetParameter(2);
      mean.push_back(Mean);
      sigma.push_back(Sigma);
    }

    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      int n;
      if (logscale){
        n=(int)(Nbins*(log(i->GetMeas())-log(minY))/(log(maxY)-log(minY)));
        if (n>=Nbins) n=Nbins-1;
      } else {
        n=(int)(Nbins*(i->GetMeas()-minY)/d);
        if (n>=Nbins) n=Nbins-1;
      };
      double delta = mean.at(n) - i->GetTruth();
      //Use only events with truth value in 3 sigma of mean truth... ???
      if (outlier && fabs(delta/sigma.at(n))>3.) continue;
      double A, B;
      double ERR = error(i->GetTruth()+2*delta);
      if (ERR > 0){
        A = 1./sqrt(2*TMath::Pi())/i->GetError()*exp(-pow(delta,2)/2./pow(i->GetError(),2));
        B = 1./sqrt(2*TMath::Pi())/ERR*exp(-pow(delta,2)/2./pow(ERR,2));
        double newweight = 2*B/(A+B);
        i->SetWeight(i->GetWeight_orig()*newweight);
      }
      dataCut.push_back(*i);
    }
    data.clear();
    data.assign(dataCut.begin(),dataCut.end());
    cout<<"Events at end of ResolutionReweight: "<<data.size()<<endl;

    minY=99999.;
    maxY=0;
    minX=99999.;
    maxX=0;
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      if (i->GetMeas() < minY) minY=i->GetMeas();
      if (i->GetMeas() > maxY) maxY=i->GetMeas();
      if (i->GetTruth() < minX) minX=i->GetTruth();
      if (i->GetTruth() > maxX) maxX=i->GetTruth();
    }

  }
  //----------------------------------------------------------

  //----------------------------------------------------------
  //symmetrize truth distribution for one measurement
  void CorrectCutoff(bool logscale = false){

    cout<<"Events at start of CorrectCutoff: "<<data.size()<<endl;
    int Nbins=49;
    int minEntries=200;
    vector<TDataPoint> dataCut;
    vector<double> mean;
    vector<double> dmean;

    TCanvas *c1 = new TCanvas("c1","Correct for cut off",0,0,1000,1000);
    c1->Divide(7,7);
    double d;
    if (logscale){
      d=(log(maxY)-log(minY))/Nbins;
    } else {
      d=(maxY-minY)/Nbins;
    }

    //book histos
    TH1F* hist[Nbins];
    TF1* f[Nbins];
    for (int n=0; n<Nbins; ++n){
      char hname[100];
      sprintf(hname,"hc%i",n);
      hist[n] = new TH1F(hname,"Truth for MeasBin", 1000, minX, maxX);
      hist[n]->Sumw2();
    }

    //Fill histos
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      int n;
      if (logscale){
        n=(int)(Nbins*(log(i->GetMeas())-log(minY))/(log(maxY)-log(minY)));
        if (n>=Nbins) n=Nbins-1;
      } else {
        n=(int)(Nbins*(i->GetMeas()-minY)/d);
        if (n>=Nbins) n=Nbins-1;
      };
      hist[n]->Fill(i->GetTruth(),i->GetWeight());
    }

    //Fit histos
    for (int n=0; n<Nbins; ++n){
      cout<<n<<" of "<<Nbins<<endl;
      if (hist[n]->GetEntries()<minEntries){
        mean.push_back(0);
        dmean.push_back(0);
        continue;
      }
      char fname[100];
      sprintf(fname,"fc%i",n);
      f[n] = new TF1(fname,"gausn",minX,maxX);
      //use a guess for starting values
      f[n]->SetParameters(hist[n]->Integral("width"),hist[n]->GetMean(),hist[n]->GetRMS());
      c1->cd(n+1);
      hist[n]->SetLineWidth(0);
      hist[n]->Draw();
      f[n]->SetLineWidth(0);
      f[n]->SetLineColor(2);
      hist[n]->Fit(fname,"ILLRQ");
      double Mean=f[n]->GetParameter(1);
      mean.push_back(Mean);
      double dMean=0;
      if (Mean<minX || Mean>maxX){
        dMean=0;
      } else {
        double dMeanMin=fabs(Mean-minX);
        double dMeanMax=fabs(Mean-maxX);
        dMeanMin<dMeanMax?dMean=dMeanMin:dMean=dMeanMax;
      }
      dmean.push_back(dMean);
    }
    c1->SaveAs("cut_LVMini.eps");

    //symmetrize distributions
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      int n;
      if (logscale){
        n=(int)(Nbins*(log(i->GetMeas())-log(minY))/(log(maxY)-log(minY)));
        if (n>=Nbins) n=Nbins-1;
      } else {
        n=(int)(Nbins*(i->GetMeas()-minY)/d);
        if (n>=Nbins) n=Nbins-1;
      };
      double delta = fabs(mean.at(n) - i->GetTruth());
      if (delta>dmean.at(n)) continue;
      dataCut.push_back(*i);
    }

    data.clear();
    data.assign(dataCut.begin(),dataCut.end());
    cout<<"Events at end of CorrectCutoff: "<<data.size()<<endl;

    minY=99999.;
    maxY=0;
    minX=99999.;
    maxX=0;
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      if (i->GetMeas() < minY) minY=i->GetMeas();
      if (i->GetMeas() > maxY) maxY=i->GetMeas();
      if (i->GetTruth() < minX) minX=i->GetTruth();
      if (i->GetTruth() > maxX) maxX=i->GetTruth();
    }
    
  }
  //----------------------------------------------------------

public:

  vector<TDataPoint>* GetDataRef(){return &data;};

  //----------------------------------------------------------
  void DrawTruth(){
    TH1F* hist = new TH1F("truth","True distribution of measurement", 100, 0, maxX);
    hist->Sumw2();
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      hist->Fill(i->GetTruth(),i->GetWeight());
    }
    TCanvas *c1 = new TCanvas("c1","Truth",0,0,600,600);
    c1->Divide(1,1);
    c1->cd(1);
    c1->cd(1)->SetLogy();
    hist->SetXTitle("Truth");
    hist->Draw();
    c1->SaveAs("truth_LVMini.eps");
  };
  //----------------------------------------------------------

  //----------------------------------------------------------
  void DrawMeas(){
    TH1F* hist = new TH1F("meas","Measured distribution of measurement", 100, 0, maxY);
    hist->Sumw2();
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      hist->Fill(i->GetMeas(),i->GetWeight());
    }
    TCanvas *c1 = new TCanvas("c2","Measurement",0,0,600,600);
    c1->Divide(1,1);
    c1->cd(1);
    c1->cd(1)->SetLogy();
    hist->SetXTitle("Measurement");
    hist->Draw();
    c1->SaveAs("meas_LVMini.eps");
  };
  //----------------------------------------------------------
  
  //----------------------------------------------------------
  void DrawCalibConst(vector<double> p){
    TGraph cgr(sp.size(), &sp.front(), &p.front());
    TSpline3 *cspl = new TSpline3("cspl",&cgr);
    vector<double> vx, vy;
    int Nbins = 100;
    for (int i=0; i<=Nbins; ++i){
      double x=minY+i*(maxY-minY)/Nbins;
      vx.push_back(double(x));
      double c=calib(x,cspl);
      vy.push_back(c);
    }
    TGraph* CalibConst = new TGraph(vx.size(), &vx.front(), &vy.front());
    TCanvas *c1 = new TCanvas("c3","Calibration Constant vs. Measurement",0,0,600,600);
    c1->Divide(1,1);
    c1->cd(1);
    CalibConst->SetLineColor(2);
    CalibConst->SetLineWidth(4);
    CalibConst->SetTitle("Calibration Constant");
    CalibConst->GetXaxis()->SetTitle("Truth");
    CalibConst->GetYaxis()->SetTitle("C");
    CalibConst->SetMinimum(0.9);
    CalibConst->SetMaximum(1.3);
    CalibConst->Draw("ACP");
    //CalibConst->Draw("L");
    c1->SaveAs("calibconst_LVMini.eps");
  };
  //----------------------------------------------------------

  //----------------------------------------------------------
  void DrawRatio(vector<double> p){
    TGraph cgr(sp.size(), &sp.front(), &p.front());
    TSpline3 *cspl = new TSpline3("cspl",&cgr);
    TH2F* hist1 = new TH2F("r1","Measurement/Truth vs. Truth",                 100, 0, 1000, 100, 0.0, 2.0);
    TH2F* hist2 = new TH2F("r2","Measurement/Truth vs. Measurement",           100, 0, 1000, 100, 0.0, 2.0);
    TH2F* hist3 = new TH2F("r3","Measurement/Truth vs. CalibMeasurement",      100, 0, 1000, 100, 0.0, 2.0);
    TH2F* hist4 = new TH2F("r4","CalibMeasurement/Truth vs. Truth",            100, 0, 1000, 100, 0.0, 2.0);
    TH2F* hist5 = new TH2F("r5","CalibMeasurement/Truth vs. Measurement",      100, 0, 1000, 100, 0.0, 2.0);
    TH2F* hist6 = new TH2F("r6","CalibMeasurement/Truth vs. CalibMeasurement", 100, 0, 1000, 100, 0.0, 2.0);
    hist1->Sumw2();
    hist2->Sumw2();
    hist3->Sumw2();
    hist4->Sumw2();
    hist5->Sumw2();
    hist6->Sumw2();
    hist1->SetTitleOffset(1.2,"Y");
    hist2->SetTitleOffset(1.2,"Y");
    hist3->SetTitleOffset(1.2,"Y");
    hist4->SetTitleOffset(1.2,"Y");
    hist5->SetTitleOffset(1.2,"Y");
    hist6->SetTitleOffset(1.2,"Y");
    hist1->SetTitleOffset(1.2,"Y");
    hist1->SetYTitle("k=Meas/Truth");
    hist2->SetYTitle("k=Meas/Truth");
    hist3->SetYTitle("k=Meas/Truth");
    hist4->SetYTitle("k=calibMeas/Truth");
    hist5->SetYTitle("k=calibMeas/Truth");
    hist6->SetYTitle("k=calibMeas/Truth");
    hist1->SetXTitle("Truth");
    hist4->SetXTitle("Truth");
    hist2->SetXTitle("Measurement");
    hist5->SetXTitle("Measurement");
    hist3->SetXTitle("calib. Measurement");
    hist6->SetXTitle("calib. Measurement");
    for (vector<TDataPoint>::iterator i=data.begin(); i<data.end(); ++i){
      double mess=i->GetMeas();
      double truth=i->GetTruth();
      double c=calib(mess,cspl);
      hist1->Fill(truth,  mess/truth,     i->GetWeight());
      hist2->Fill(mess,   mess/truth,     i->GetWeight());
      hist3->Fill(c*mess, mess/truth,     i->GetWeight());
      hist4->Fill(truth,  (c*mess)/truth, i->GetWeight());
      hist5->Fill(mess,   (c*mess)/truth, i->GetWeight());
      hist6->Fill(c*mess, (c*mess)/truth, i->GetWeight());
    }
    gROOT->SetStyle("Plain");
    gStyle->SetPalette(51,0);
    TCanvas *c1 = new TCanvas("c4","Ration",0,0,900,900);
    TF1 *f = new TF1("f","gaus",minX,maxX);
    TObjArray aSlices;
    c1->Divide(3,3);
    c1->cd(1);
    gPad->SetRightMargin(0.15);
    gPad->SetLeftMargin(0.15);
    hist1->Draw("COLZ");
    hist1->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist1_1 = (TH1*)(aSlices[1]);
    TH1* hist1_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist1_1->GetNbinsX(); ++i){
     hist1_1->SetBinError(i,hist1_2->GetBinContent(i));
    }
    hist1_1->SetLineColor(2);
    hist1_1->Draw("same");
    hist1->SetLineColor(5);
    hist1->ProfileX()->Draw("same");

    c1->cd(2);
    gPad->SetRightMargin(0.12);
    hist2->Draw("COLZ");
    hist2->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist2_1 = (TH1*)(aSlices[1]);
    TH1* hist2_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist2_1->GetNbinsX(); ++i){
     hist2_1->SetBinError(i,hist2_2->GetBinContent(i));
    }
    hist2_1->SetLineColor(2);
    hist2_1->Draw("same");
    hist2->SetLineColor(5);
    hist2->ProfileX()->Draw("same");

    c1->cd(3);
    gPad->SetRightMargin(0.12);
    hist3->Draw("COLZ");
    hist3->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist3_1 = (TH1*)(aSlices[1]);
    TH1* hist3_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist3_1->GetNbinsX(); ++i){
      hist3_1->SetBinError(i,hist3_2->GetBinContent(i));
    }
    hist3_1->SetLineColor(2);
    hist3_1->Draw("same");
    hist3->SetLineColor(5);
    hist3->ProfileX()->Draw("same");

    c1->cd(4);
    gPad->SetRightMargin(0.12);
    hist4->Draw("COLZ");
    hist4->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist4_1 = (TH1*)(aSlices[1]);
    TH1* hist4_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist4_1->GetNbinsX(); ++i){
     hist4_1->SetBinError(i,hist4_2->GetBinContent(i));
    }
    hist4_1->SetLineColor(2);
    hist4_1->Draw("same");
    hist4->SetLineColor(5);
    hist4->ProfileX()->Draw("same");

    c1->cd(5);
    gPad->SetRightMargin(0.12);
    hist5->Draw("COLZ");
    hist5->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist5_1 = (TH1*)(aSlices[1]);
    TH1* hist5_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist5_1->GetNbinsX(); ++i){
     hist5_1->SetBinError(i,hist5_2->GetBinContent(i));
    }
    hist5_1->SetLineColor(2);
    hist5_1->Draw("same");
    hist5->SetLineColor(5);
    hist5->ProfileX()->Draw("same");

    c1->cd(6);
    gPad->SetRightMargin(0.12);
    hist6->Draw("COLZ");
    hist6->FitSlicesY(f, 0, -1, 0, "NQ", &aSlices);
    TH1* hist6_1 = (TH1*)(aSlices[1]);
    TH1* hist6_2 = (TH1*)(aSlices[2]);
    for (int i=0; i<hist6_1->GetNbinsX(); ++i){
     hist6_1->SetBinError(i,hist6_2->GetBinContent(i));
    }
    hist6_1->SetLineColor(2);
    hist6_1->Draw("same");
    hist6->SetLineColor(5);
    hist6->ProfileX()->Draw("same");

    c1->cd(7);
    hist4_1->SetTitleOffset(1.2,"Y");
    hist4_1->SetXTitle("Truth");
    hist4_1->SetYTitle("k");
    hist4_1->SetLineWidth(0);
    hist4_1->SetLineColor(2);
    hist4_1->SetMinimum(0.95);
    hist4_1->SetMaximum(1.05);
    hist4_1->Draw();

    c1->cd(8);
    hist5_1->SetTitleOffset(1.2,"Y");
    hist5_1->SetXTitle("Measurement");
    hist5_1->SetYTitle("k");
    hist5_1->SetLineWidth(0);
    hist5_1->SetLineColor(2);
    hist5_1->SetMinimum(0.95);
    hist5_1->SetMaximum(1.05);
    hist5_1->Draw();
    
    c1->cd(9);
    hist6_1->SetTitleOffset(1.2,"Y");
    hist6_1->SetXTitle("calib. Measurement");
    hist6_1->SetYTitle("k");
    hist6_1->SetLineWidth(0);
    hist6_1->SetLineColor(2);
    hist6_1->SetMinimum(0.95);
    hist6_1->SetMaximum(1.05);
    hist6_1->Draw();

    c1->SaveAs("ratio_LVMini.eps");
  };
  //----------------------------------------------------------

  //----------------------------------------------------------
  //constructor
  TData(int N=100, double min=300., double max=700., int npar=10, 
        bool flat=false, 
        bool cutoff=false, bool cutofflog=false, 
        bool reweight=false, bool reweightlog=false, 
        bool reco=false){
    Nevts = N;
    minX=99999;
    maxX=0;
    minY=99990;
    maxY=0;
    int Nnegative=0;
    TF1 *recoeff = new TF1("f","pow(x/30.,2)/(1.+pow(x/30.,2))");
    for (int i=0;i<Nevts;i++) {
      double x = gRandom->Rndm(i);
      //double truth = min+(max-min)*x;
      double truth = gRandom->Exp(150);
      double width = error(truth);
      double meas  = gRandom->Gaus(0.9*truth,width);
      double pr = gRandom->Rndm(i);
      //double pr = -9999.;
      if (meas<0 && truth>min){
        ++Nnegative;
        cout<<"WARNING: negative measurement (truth="<<truth<<",meas="<<meas<<") no. "<<Nnegative<<endl;
      }
      if (meas>0 && truth>min && pr<(*recoeff)(truth)){
        if (meas<minY) minY=meas;
        if (meas>maxY) maxY=meas;
        if (truth<minX) minX=truth;
        if (truth>maxX) maxX=truth;
        double error = width;
        double weight=1.;
        //reco eff correction
        if (reco) weight = 1./(*recoeff)(truth);
        TDataPoint d(truth, meas, error, weight);
        data.push_back(d);
      }
    }
    sp.clear();

    if (flat) FlattenSpectra();
    if (reweight){
      bool outlier=true;
      ResolutionReweight(reweightlog, outlier, 5);
      ResolutionReweight(reweightlog, outlier, 4);
      ResolutionReweight(reweightlog, outlier, 3);
    }
    if (cutoff) CorrectCutoff(cutofflog);

    dY=(maxY-minY)/(npar-3);
    for (int i=0; i<npar; ++i){
      sp.push_back(minY+(i-1)*dY);
    }

  };
  //----------------------------------------------------------
  
  //----------------------------------------------------------
  void evalF(int &NPAR, vector<double>* GRAD, double &FSUM, vector<double>* PAR){
    
    FSUM = chi2(PAR);
    // gradients
    double eps=1.e-5;
    for (int i=0; i<NPAR; ++i){
      vector<double> PAR1 = *PAR;
      vector<double> PAR2 = *PAR;
      PAR1.at(i) += eps;
      PAR2.at(i) -= eps;
      GRAD->at(i) = (chi2(&PAR1)-chi2(&PAR2))/2./eps;
      PAR1.at(i) += eps;
      PAR2.at(i) -= eps;
      GRAD->at(i+NPAR) = (chi2(&PAR1)-2*FSUM+chi2(&PAR2))/4./eps/eps;
    }
  }
  //----------------------------------------------------------

  //destructor
  ~TData(){};

};
////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////
int main(){

  int NPAR=11;

  TROOT simple("simple","Some calibration toyMC tests");
  bool flat        = true;
  bool cutoff      = true;
  bool cutofflog   = true;
  bool reweight    = true;
  bool reweightlog = true;
  bool reco        = true;
  TData* toy= new TData(1000000,20,99999,NPAR,flat,cutoff,cutofflog,reweight,reweightlog,reco);
  toy->DrawTruth();
  toy->DrawMeas();

  ////// LVMINI //////
  int NITER = 500;
  int MVEC;
  NPAR<29?MVEC=NPAR:MVEC=29;
  int NAUX = 10000;

  vector<double> AUX(NAUX, 0.);
  vector<double> parameter;
  double FSUM;
  double FOPT, FEDM, DUMMY;

  for (int i=0; i<NPAR; ++i){
    parameter.push_back(1./0.9); 
  }
  float eps =float(1.E-3);
  float wlf1=1.E-4;
  float wlf2=0.9;
  //default values
  //double wlf1=-1.;
  //double wlf2=-1.;

  lvmeps_(eps,wlf1,wlf2);

  int NPARNEG = -fabs(NPAR);
  if(lvmdim_(NPAR, MVEC) > NAUX)
    cout<<"Aux field too small: "<<NAUX<<"<"<<lvmdim_(NPAR, MVEC)<<" needed entries"<<endl;
  lvmini_(NPARNEG, MVEC, NITER, &AUX.front());
  int IFLAG;
  int IRET = 0;
  int ITER = 0;
  do {
    ++ITER;
    toy->evalF(NPAR, &AUX, FSUM, &parameter);
    lvmfun_(&parameter.front(), FSUM, IRET, &AUX.front());
    lvmprt_(2, &AUX.front(), 2); //print out
    /*
    if (ITER%10==0){
      toy->DrawRatio(parameter);
      toy->DrawCalibConst(parameter);
    }
    */
  } while (IRET<0);
  int error_index=2;
  error_index = lvmind_(error_index);

  toy->DrawRatio(parameter);
  toy->DrawCalibConst(parameter);

  return 0;

}
////////////////////////////////////////////////////////////

Revision:	1.1.1.1 (vendor branch)
Committed:	Tue Nov 4 13:41:08 2008 UTC (16 years, 6 months ago) by csander
Branch:	ToyMc, MAIN
CVS Tags:	start, HEAD
Changes since 1.1:	+0 -0 lines
Error occurred while calculating annotation data.
Log Message:	program using toyMC to studie gamma jet calibration (using root and LVMini)