PWGLF/totEt/macros/hadEt/Unfoldpp.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <iostream>
   3 using std::cout;
   4 using std::endl;
   5
   6 #include "TRandom.h"
   7 #include "TH1D.h"
   8
   9 #include "RooUnfoldResponse.h"
  10 #include "RooUnfoldBayes.h"
  11 //#include "RooUnfoldDagostini.h"
  12 #include "RooUnfoldSvd.h"
  13 //#include "RooUnfoldTUnfold.h"
  14 #include "RooUnfoldErrors.h"
  15 #endif
  16
  17 TF1 *myGaussian = new TF1("Gaussian","1/[0]/TMath::Sqrt(2*TMath::Pi())*exp(-0.5*(x/[0])**2)",-5,5);
  18 //gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2) and (0) means start numbering parameters at 0.
  19 TH1F *hSim;
  20 TH2F *resolutionFull;
  21 TH2F *resolution;
  22 //TClonesArray *histoarrayCopy = new TClonesArray("TH1D",200);
  23 TClonesArray histoarray("TH1D",200);
  24 TString *ytitle;
  25 TH1F *GetReconstructedEt(TF1 *inputfunc, int nevents, char *name, float lowbound, float highbound, int nbins);
  26 TH1F *GetReconstructedEt(TH1 *input, int nevents, char *name, float lowbound, float highbound, int nbins, bool smooth);
  27 TH1F *GetTrueEt(TF1 *inputfunc, int nevents, char *name, float lowbound, float highbound, int nbins);
  28 Float_t GetChi2(TH1F *reconstructed, TH1F *simulated);
  29 TH1F *RestrictAxes(TH1F *old,float minEt,float maxEt);
  30 TH2F *RestrictAxes(TH1F *old,float minEt,float maxEt);
  31 void Smooth(TH1F *, TF1 *,bool);
  32 //variables I had included in arguments but which are now not used
  33 bool zerolowetbins = false;
  34 int minbin = 4;
  35 bool its = true;
  36 int difftype = 0;
  37 int niter = 3;
  38 bool rescaleguess = false;
  39 //For both training and testing cases the code is:
  40 //0 = MC
  41 //1 = reweighted MC
  42 //2 = exponential
  43 //3 = flat distribution
  44 //4 = Straight line
  45 //5 = Half Gaussian
  46 //6 = double exponential
  47 //7 = Levy function
  48 //8 = Power law
  49 //for had =
  50 //0 = tot et
  51 //1 = had et
  52 //2 = pikp et
  53 //testmean is
  54 //the mean for an exponential
  55 //A for a Levy function and for a power function
  56 //testn is the input n for the Levy and Power functions
  57 void Unfoldpp(int simdataset = 20111,int recodataset = 20111, char *outputfilename = "junk", int had = 0, int trainingcase= 7,int neveUsed = 1e6,float minEt = 0.00, float maxEt = 50.0, int varymean=0, bool unfoldData = true, float testmean = 5,float testn = -4,bool smooth = true){
  58
  59 int unfoldcase = trainingcase;
  60 myGaussian->SetParameter(0,1.0);
  61 #ifdef __CINT__
  62   gSystem->Load("~/alicework/RooUnfold-1.1.1/libRooUnfold");
  63   //gSystem->Load("libRooUnfold");
  64 #endif
  65   gStyle->SetOptTitle(0);
  66   gStyle->SetOptStat(0);
  67   gStyle->SetOptFit(0);
  68   char *infilename = NULL;
  69   char *datainfilename = NULL;
  70   switch(simdataset){
  71   case 20111:
  72     infilename="rootFiles/LHC11b10a/Et.ESD.new.sim.LHC11b10a.root";
  73     break;
  74   case 2009:
  75     infilename="rootFiles/LHC11b1a/Et.ESD.new.sim.LHC11b1a.root";
  76     break;
  77   case 2010:
  78     infilename="rootFiles/LHC10e20/Et.ESD.new.sim.LHC10e20.root";
  79     break;
  80   }
  81   switch(recodataset){
  82   case 20111:
  83     datainfilename="rootFiles/LHC11a/Et.ESD.new.sim.LHC11a.root";
  84     break;
  85   case 2009:
  86     datainfilename="rootFiles/LHC10c/Et.ESD.new.sim.LHC10c.root";
  87     break;
  88   case 2010:
  89     datainfilename="rootFiles/LHC10e/Et.ESD.new.sim.LHC10e.root";
  90     break;
  91   }
  92
  93   //================READING IN HISTOGRAMS===========================
  94   TFile *outfile = new TFile("junk.root","RECREATE");
  95   TFile *datafile = new TFile(datainfilename);
  96   TFile *file = new TFile(infilename);
  97   TList *list = file->FindObject("out2");
  98   char histoname[200];
  99   TString *hadStr;
 100   TString *longHadStr;
 101   TString *tail;
 102   TString *detector;
 103   ytitle = new TString("1/N_{eve}dN/dE_{T}");
 104   if(had==1){
 105     hadStr = new TString("Had");
 106     longHadStr = new TString("E_{T}^{had}");
 107   }
 108   if(had==0){
 109     hadStr = new TString("Tot");
 110     longHadStr = new TString("E_{T}^{tot}");
 111   }
 112   if(had==2){
 113     hadStr = new TString("PiKP");
 114     longHadStr = new TString("E_{T}^{#pi,K,p}");
 115   }
 116   switch(difftype){
 117   case 0:
 118     tail = new TString("ND");//Non-diffractive
 119     break;
 120   case 1:
 121     tail = new TString("SD");//Singly-diffractive
 122     break;
 123   case 2:
 124     tail = new TString("DD");//Doubly-diffractive
 125     break;
 126   default:
 127     tail = new TString("");//none
 128   }
 129   if(its) detector = new TString("ITS");
 130   else{ detector = new TString("TPC");}
 131   sprintf(histoname,"Sim%sEt%s",hadStr->Data(),tail->Data());
 132   //sprintf(histoname,"Sim%sEt",hadStr->Data());
 133   file->cd();
 134   hTemp = (TH1F*) out2->FindObject(histoname);
 135   outfile->cd();
 136   hSim = (TH1F*) hTemp->Clone(Form("%sCopy",histoname));
 137
 138   file->cd();
 139   sprintf(histoname,"Reco%sEtFullAcceptanceITS%s",hadStr->Data(),tail->Data());
 140   hTemp = (TH1F*) out2->FindObject(histoname);
 141   outfile->cd();
 142   hITS = (TH1F*) hTemp->Clone(Form("%sSim",histoname));
 143
 144   datafile->cd();
 145   sprintf(histoname,"Reco%sEtFullAcceptanceITS",hadStr->Data());
 146   hTemp = (TH1F*) out2->FindObject(histoname);
 147   if(unfoldData && !hTemp){ cerr<<"no histogram "<<histoname<<endl; return;}
 148   //   outfile->cd();
 149   hITSData = (TH1F*) hTemp->Clone(Form("%sData",histoname));
 150
 151   file->cd();
 152   sprintf(histoname,"Reco%sEtFullAcceptanceTPC%s",hadStr->Data(),tail->Data());
 153   //sprintf(histoname,"Reco%sEtFullAcceptanceTPC",hadStr->Data());
 154   //cout<<"Numerator "<<histoname<<endl;
 155   hTemp = (TH1F*)  out2->FindObject(histoname);
 156   outfile->cd();
 157   hTPC = (TH1F*) hTemp->Clone(Form("%sSim",histoname));
 158
 159
 160   datafile->cd();
 161   sprintf(histoname,"Reco%sEtFullAcceptanceTPC",hadStr->Data());
 162   hTemp = (TH1F*) out2->FindObject(histoname);
 163   if(unfoldData && !hTemp){ cerr<<"no histogram "<<histoname<<endl; return;}
 164   //   outfile->cd();
 165   hTPCData = (TH1F*) hTemp->Clone(Form("%sData",histoname));
 166   file->cd();
 167
 168
 169   hSim->Sumw2();
 170   hITS->Sumw2();
 171   hTPC->Sumw2();
 172   hITSData->Sumw2();
 173   hTPCData->Sumw2();
 174   int rebin = 1;
 175   //   if(had==2){
 176   //     //hSim->Rebin(rebin);
 177   //     hTPC->Rebin(rebin);
 178   //     hITS->Rebin(rebin);
 179   //   }
 180
 181   outfile->cd();
 182   //TH1F *hITSClone = (TH1F*)hITS->Clone("hITSClone");
 183   file->cd();
 184
 185   //Ex SimTotEtMinusRecoEtFullAcceptanceITS
 186   //sprintf(histoname,"Sim%sEtMinusReco%sEtFullAcceptance%s",hadStr->Data(),hadStr->Data(),detector->Data());
 187   sprintf(histoname,"Sim%sEtVsReco%sEtFullAcceptance%s",hadStr->Data(),hadStr->Data(),detector->Data());
 188   hTemp2 = (TH2F*)out2->FindObject(histoname);
 189   outfile->cd();
 190   resolutionFull = (TH2F*) hTemp2->Clone(Form("%sFull",histoname));
 191
 192   resolutionFull->Rebin2D(rebin,rebin);
 193   hSim->Rebin(rebin);
 194   hITS->Rebin(rebin);
 195   if(hITSData) hITSData->Rebin(rebin);
 196   hTPC->Rebin(rebin);
 197   if(hTPCData) hTPCData->Rebin(rebin);
 198   cout<<"Rebinning "<<rebin<<"x"<<endl;
 199   //float minEt = 0.000;
 200   //float maxEt = 100.00;
 201   //   hSim->GetXaxis()->SetRange(minbin,maxbin);
 202   //   hITS->GetXaxis()->SetRange(minbin,maxbin);
 203   //   hTPC->GetXaxis()->SetRange(minbin,maxbin);
 204   //   hITSData->GetXaxis()->SetRange(minbin,maxbin);
 205   //   hTPCData->GetXaxis()->SetRange(minbin,maxbin);
 206   //   resolution->GetXaxis()->SetRange(minbin,maxbin);
 207   //   resolution->GetYaxis()->SetRange(minbin,maxbin);
 208
 209   //cout<<"I run from "<<hSim->GetBinLowEdge(1)<<" to "<<hSim->GetBinLowEdge(hSim->GetNbinsX()+1)<<endl;
 210
 211
 212   //cout<<"Histo "<<histoname<<endl;
 213   TCanvas *canvas0 = new TCanvas("canvas0","canvas0",600,600);
 214   canvas0->SetTopMargin(0.020979);
 215   canvas0->SetRightMargin(0.0184564);
 216   canvas0->SetLeftMargin(0.0989933);
 217   canvas0->SetBottomMargin(0.101399);
 218   canvas0->SetBorderSize(0);
 219   canvas0->SetFillColor(0);
 220   canvas0->SetFillColor(0);
 221   canvas0->SetBorderMode(0);
 222   canvas0->SetFrameFillColor(0);
 223   canvas0->SetFrameBorderMode(0);
 224 //   myGaussian->Draw();
 225 //   return;
 226   canvas0->SetLogz();
 227   if(smooth)resolutionFull->Smooth(5,"R");
 228   resolutionFull->Draw("colz");
 229   //return;
 230   //if(had==2)  resolution->Rebin2D(rebin,rebin);
 231   int nbins =  resolutionFull->GetXaxis()->GetNbins();
 232   if(zerolowetbins){
 233     for(int j=0;j<minbin;j++){//loop over bins in y=reconstructed et
 234       for(int i=0;i<nbins;i++){
 235         resolutionFull->SetBinContent(i,j,0.0);
 236       }
 237     }
 238   }
 239   for(int i=0;i<nbins;i++){//loop over bins in x
 240     histoarray[i]=(TH1D*)resolutionFull->ProjectionY(Form("tmp%i",i+1),i+1,i+1);
 241   }
 242   float lowrange = 0.0;
 243   float highrange = 100.0;
 244   file->cd();
 245   file->Close();
 246   outfile->cd();
 247
 248   //================FILLING TRAINING HISTOGRAMS===========================
 249   int neveUnused = 1e3;
 250   //int neveUsed = 1e7;
 251   //~~~~~~~~~~~~~~~~~~~~~~~~~~~EXPONENTIAL~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 252   //TF1 *func = new TF1("func","([0]+x*[2])/[1]*exp(-x/[1])",0,100);
 253   TF1 *func = new TF1("func","([0])/[1]*exp(-x/[1])",0.0,100);
 254   //TF1 *func = new TF1("func","([0])/[1]*exp(-x/[1])",0,100);
 255   func->SetParameter(0,1.0);
 256   func->SetParameter(0,1);
 257   //func->SetParameter(1,1.0/2.23876e-01);
 258   float mymean = hSim->GetMean();
 259   if(unfoldData) mymean = testmean;
 260   func->SetParameter(1,(1.0+varymean*0.3)* mymean);
 261   func->SetParameter(2,1);
 262   func->SetLineColor(4);
 263   //   TF1 *funcLong = new TF1("funcLong","([0]+[1]*x+[2]*x*x+[3]*x*x*x+x^[4])/[5]*exp(-(x**[6])/[5])",lowrange,highrange);
 264   //   funcLong->SetParameter(0,1.00467e-01);
 265   //   funcLong->SetParameter(1,-2.82339e-01);
 266   //   funcLong->SetParameter(2,-7.10366e-02);
 267   //   funcLong->SetParameter(3,1.22634e-02);
 268   //   funcLong->SetParameter(4,9.25757e-01);
 269   //   funcLong->SetParameter(5,6.77688e-01);
 270   //   funcLong->SetParameter(6,6.30298e-01);
 271   int nevents = neveUnused;//1e7;
 272   if(trainingcase==2 || unfoldcase==2) nevents = neveUsed;
 273   float lowbound = hSim->GetXaxis()->GetBinLowEdge(1);
 274   float highbound = hSim->GetXaxis()->GetBinLowEdge(nbins+1);
 275   //cout<<"Maxing histograms with ranges "<<lowbound<<" - "<<highbound<<endl;
 276   TH1F *hSimExponential = GetTrueEt(func,nevents,Form("testtrue%i",i),lowbound,highbound,hITS->GetNbinsX());
 277   TH1F *hMeasuredExponential = GetReconstructedEt(func,nevents,Form("testsmeared%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 278
 279   //~~~~~~~~~~~~~~~~~~~~~~~~~~~FLAT~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 280   TF1 *funcFlat = new TF1("funcFlat","[0]",0,100);
 281   funcFlat->FixParameter(0,0.01);
 282   nevents = neveUnused;//1e7;
 283   if(trainingcase==3 || unfoldcase==3) nevents = neveUsed;
 284   TH1F *hSimFlat = GetTrueEt(funcFlat,nevents,Form("testtrueflat%i",i),lowbound,highbound,hITS->GetNbinsX());
 285   TH1F *hMeasuredFlat = GetReconstructedEt(funcFlat,nevents,Form("testsmearedflat%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 286
 287   //~~~~~~~~~~~~~~~~~~~~~~~~~~~STRAIGHT LINE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 288   TF1 *funcStraight = new TF1("funcStraight","[0]-[1]*x",0.0,maxEt*2);
 289   funcStraight->FixParameter(0,1.0);
 290   funcStraight->FixParameter(1,0.01);
 291   nevents = neveUnused;//1e7;
 292   if(trainingcase==4 || unfoldcase==4) nevents = neveUsed;
 293   TH1F *hSimStraight = GetTrueEt(funcStraight,nevents,Form("testtruestraight%i",i),lowbound,highbound,hITS->GetNbinsX());
 294   TH1F *hMeasuredStraight = GetReconstructedEt(funcStraight,nevents,Form("testsmearedstraight%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 295
 296
 297   //~~~~~~~~~~~~~~~~~~~~~~~~~~~GAUS LINE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 298   TF1 *funcGaus = new TF1("funcGaus","[0]*exp(-x*x/[1]/[1]/TMath::Pi())",0.0,maxEt*2);//[1] is the mean x
 299   funcGaus->FixParameter(0,1.0);
 300   funcGaus->FixParameter(1,15);
 301   nevents = neveUnused;//1e7;
 302   if(trainingcase==5 || unfoldcase==5) nevents = neveUsed;
 303   TH1F *hSimGaus = GetTrueEt(funcGaus,nevents,Form("testtruegaus%i",i),lowbound,highbound,hITS->GetNbinsX());
 304   TH1F *hMeasuredGaus = GetReconstructedEt(funcGaus,nevents,Form("testsmearedgaus%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 305
 306   //~~~~~~~~~~~~~~~~~~~~~~~~~~~DOUBLE EXPONENT LINE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 307   TF1 *funcDoubleExponential = new TF1("funcDoubleExponential","[0]/[1]*exp(-x/[1])+[2]/[3]*exp(-x/[3])",0.0,maxEt*2);
 308   //TF1 *funcDoubleExponential = new TF1("funcDoubleExponential","[0]/[1]*exp(-x/[1])",0,100);
 309   funcDoubleExponential->SetParameter(0,1.0);
 310   funcDoubleExponential->SetParameter(1,testmean);
 311   funcDoubleExponential->FixParameter(2,0.1);
 312   funcDoubleExponential->FixParameter(3,1.0);
 313   funcDoubleExponential->SetLineColor(5);
 314   nevents = neveUnused;//1e7;
 315   if(trainingcase==6 || unfoldcase==6) nevents = neveUsed;
 316   TH1F *hSimDoubleExponential = GetTrueEt(funcDoubleExponential,nevents,Form("testtruedoubleexponent%i",i),lowbound,highbound,hITS->GetNbinsX());
 317   TH1F *hMeasuredDoubleExponential = GetReconstructedEt(funcDoubleExponential,nevents,Form("testsmeareddoubleexponent%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 318
 319
 320   //~~~~~~~~~~~~~~~~~~~~~~~~~~~TRUE SIMULATED~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 321   //trying to make something which is not identical to hTPC/hITS
 322   //if(trainingcase==0 || unfoldcase==0) nevents = neveUsed;
 323   //TH1F *hMeasuredSimMCSmeared =  GetReconstructedEt(hSim,nevents,Form("testtruemcsmeared%i",i),lowbound,highbound,hITS->GetNbinsX());
 324   if(trainingcase==0 || unfoldcase==0) TH1F *hMeasuredSimMCSmeared =  hITS;
 325   //~~~~~~~~~~~~~~~~~~~~~~~~~~~REWEIGHTED SIMULATED~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 326   TH1F *hSimReweighted = (TH1F*) hSim->Clone("hSimReweighted");
 327   TF1 *fWeight = new TF1("fWeight","1-[0]*([1]-x)+[2]*([3]-x)**2",0.0,2.0*maxEt);
 328   fWeight->SetParameter(0,-0.0005);
 329   fWeight->SetParameter(1,3);
 330   fWeight->SetParameter(2,-0.0005);
 331   fWeight->SetParameter(3,3);
 332   hSimReweighted->Divide(fWeight);
 333   nevents = neveUnused;//1e7;
 334   if(trainingcase==1 || unfoldcase==1) nevents = neveUsed;
 335   TH1F *hMeasReweighted =  GetReconstructedEt(hSimReweighted,nevents,Form("testreweighted%i",i),lowbound,highbound,hITS->GetNbinsX());
 336   //float chi2 = GetChi2(hITS,test);
 337
 338
 339   //~~~~~~~~~~~~~~~~~~~~~~~~~~~LEVY~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 340   TF1 *funcLevy = new TF1("funcLevy","[0]*x*(1+x/[1])**[2]",0,100);
 341   //TF1 *funcLevy = new TF1("funcLevy","[0]*x*([1]+x+[3]*x*x)**[2]",0,100);
 342   funcLevy->SetParameter(0,1.0);
 343   float n = testn;//-4;
 344   float A = 2;//- (1.0+varymean*0.3)* mymean *(n+3)/(n*n-4*n+5)*50;
 345   //cout<<"Mean "<<mymean<<" A "<<A<<" n "<<n<<endl;
 346   funcLevy->SetParameter(1,(1.0+varymean*0.3)* mymean);
 347   //funcLevy->SetParameter(1,A);
 348   //funcLevy->SetParameter(1,1.0);
 349   //funcLevy->SetParameter(2,-5.96110e+00);
 350   funcLevy->SetParameter(2,n);
 351   //funcLevy->SetParameter(3,1e-1);
 352   nevents = neveUnused;//1e7;
 353   if(trainingcase==7 || unfoldcase==7) nevents = neveUsed;
 354   //   hSim->Fit(funcLevy,"","",1,100);
 355   //return;
 356   TH1F *hSimLevy = GetTrueEt(funcLevy,nevents,Form("testtruelevy%i",i),lowbound,highbound,hITS->GetNbinsX());
 357   TH1F *hMeasuredLevy = GetReconstructedEt(funcLevy,nevents,Form("testsmearedlevy%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 358   //~~~~~~~~~~~~~~~~~~~~~~~~~~~POWER~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 359   TF1 *funcPower = new TF1("funcPower","[0]*(1+x/[1])**[2]",0,100);
 360   funcPower->SetParameter(0,1.0);
 361   float n = testn;//-100;
 362   funcPower->SetParameter(1,(1.0+varymean*0.3)* mymean*(-n+2));
 363   funcPower->SetParameter(2,n);
 364   nevents = neveUnused;//1e7;
 365   //cout<<"effective scale for power function "<<funcPower->GetParameter(1)<<endl;
 366   if(trainingcase==8 || unfoldcase==8) nevents = neveUsed;
 367   //   hSim->Fit(funcLevy,"","",1,100);
 368   //   return;
 369   TH1F *hSimPower = GetTrueEt(funcPower,nevents,Form("testtruepower%i",i),lowbound,highbound,hITS->GetNbinsX());
 370   TH1F *hMeasuredPower = GetReconstructedEt(funcPower,nevents,Form("testsmearedpower%i",i),lowbound,highbound,hITS->GetNbinsX(),smooth);
 371   //================TRAINING===========================
 372   //  RooUnfoldResponse (const TH1* measured, const TH1* truth, const TH2* response, const char* name= 0, const char* title= 0);  // create from already-filled histograms
 373   TH1F *hTrainingTruth;
 374   TH1F *hTrainingReconstructed;
 375   switch(trainingcase){
 376   case 0:
 377     cout<<"Training with pure MC"<<endl;
 378     hTrainingTruth = hSim;
 379     hTrainingReconstructed = hMeasuredSimMCSmeared;
 380     //     if(its) hTrainingReconstructed = hITS;
 381     //     else{hTrainingReconstructed = hTPC;}
 382     break;
 383   case 1:
 384     cout<<"Training with reweighted MC"<<endl;
 385     hTrainingTruth = hSimReweighted;
 386     hTrainingReconstructed = hMeasReweighted;
 387     break;
 388   case 2:
 389     cout<<"Training with exponential"<<endl;
 390     hTrainingTruth = hSimExponential;
 391     hTrainingReconstructed = hMeasuredExponential;
 392     break;
 393   case 3:
 394     cout<<"Training with flat distribution"<<endl;
 395     hTrainingTruth = hSimFlat       ;
 396     hTrainingReconstructed = hMeasuredFlat       ;
 397     break;
 398   case 4:
 399     cout<<"Training with straight line distribution"<<endl;
 400     hTrainingTruth = hSimStraight       ;
 401     hTrainingReconstructed = hMeasuredStraight    ;
 402     break;
 403   case 5:
 404     cout<<"Training with half Gaussian distribution"<<endl;
 405     hTrainingTruth = hSimGaus       ;
 406     hTrainingReconstructed = hMeasuredGaus    ;
 407     break;
 408   case 6:
 409     cout<<"Training with double exponential"<<endl;
 410     hTrainingTruth = hSimDoubleExponential;
 411     hTrainingReconstructed = hMeasuredDoubleExponential;
 412     break;
 413   case 7:
 414     cout<<"Training with Levy distribution"<<endl;
 415     hTrainingTruth = hSimLevy       ;
 416     hTrainingReconstructed = hMeasuredLevy       ;
 417     break;
 418   case 8:
 419     cout<<"Training with Power law distribution"<<endl;
 420     hTrainingTruth = hSimPower       ;
 421     hTrainingReconstructed = hMeasuredPower       ;
 422     break;
 423   }
 424   hSim = RestrictAxes(hSim,minEt,maxEt);
 425   hITS = RestrictAxes(hITS,minEt,maxEt);
 426   hTPC = RestrictAxes(hTPC,minEt,maxEt);
 427   hITSData = RestrictAxes(hITSData,minEt,maxEt);
 428   hTPCData = RestrictAxes(hTPCData,minEt,maxEt);
 429   resolution = RestrictAxes(resolutionFull,minEt,maxEt);
 430   hTrainingReconstructed = RestrictAxes(hTrainingReconstructed,minEt,maxEt);
 431   hTrainingTruth = RestrictAxes(hTrainingTruth,minEt,maxEt);
 432 //   if(rescaleguess && unfoldData){
 433 //     cout<<"Rescaling the guess by MC reconstructed/data reconstructed"<<endl;
 434 //     TH1F *hScale = hTrainingReconstructed->Clone("hScale");
 435 //     if(its) hScale->Divide(hITSData);
 436 //     else{hScale->Divide(hTPCData);}
 437 //     //hTrainingReconstructed is now ~MC/data
 438 //     //TH1F *hTrainingTruthCopy = (TH1F*)hTrainingTruth->Clone("hTrainingTruthCopy");
 439 //     hTrainingTruth->Divide(hScale);//This is now MC/(MC/data) and in principle is the best guess at data?
 440 //     //hTrainingTruth->Add(hTrainingTruthCopy);
 441 //     float scale = hTrainingTruth->Integral();
 442 //     cout<<"Integral "<<scale<<endl;
 443 //     hTrainingTruth->Scale(1.0/scale);
 444
 445 //     hTrainingReconstructed = GetReconstructedEt(hTrainingTruth,nevents,Form("testMCreweightedbydata%i",i),lowbound,highbound,hITS->GetNbinsX());
 446 //   }
 447
 448   RooUnfoldResponse response(hTrainingReconstructed,hTrainingTruth,resolution,"UnfoldResponseFromHistograms","MyRooUnfoldResponse");
 449   //RooUnfoldResponse response(hTPC,hSimReweighted,resolution,"UnfoldResponseFromHistograms","MyRooUnfoldResponse");
 450   //RooUnfoldResponse response(hSmeared,hTrue,resolution,"UnfoldResponseFromHistograms","MyRooUnfoldResponse");
 451
 452   //================NORMALIZING===========================
 453   //cout<<"Normalizing..."<<endl;
 454
 455   Float_t binwidth = hSim->GetBinWidth(1);
 456   Float_t neve = hSim->GetEntries();
 457   if(neve<=1e-5){
 458     cerr<<"Warning!! simulation histo not filled!  Stopping!"<<endl;
 459     return;
 460   }
 461   if(binwidth<=1e-5){
 462     cerr<<"Warning!! binwidth!  Stopping!"<<endl;
 463     return;
 464   }
 465   //true MC
 466   hSim->Scale(1.0/neve/binwidth);
 467   hSim->Scale(1.0/hSim->Integral());
 468   neve = hTPC->GetEntries();
 469   if(neve<=1e-5){
 470     cerr<<"Warning!! Histo not filled!  Stopping!"<<endl;
 471     return;
 472   }
 473   hTPC->Scale(1.0/neve/binwidth);
 474   neve = hITS->GetEntries();
 475   hITS->Scale(1.0/neve/binwidth);
 476   //Reweighted MC
 477   neve = hMeasReweighted->GetEntries();
 478   hMeasReweighted->Scale(1.0/neve/binwidth);
 479   neve = hSimReweighted->GetEntries();
 480   hSimReweighted->Scale(1.0/neve/binwidth);
 481   //Exponential MC
 482   neve = hMeasuredExponential->GetEntries();
 483   hMeasuredExponential->Scale(1.0/neve/binwidth);
 484   neve = hSimExponential->GetEntries();
 485   hSimExponential->Scale(1.0/neve/binwidth);
 486   //Flat MC
 487   neve = hMeasuredFlat->GetEntries();
 488   hMeasuredFlat->Scale(1.0/neve/binwidth);
 489   neve = hSimFlat->GetEntries();
 490   hSimFlat->Scale(1.0/neve/binwidth);
 491   //Straight MC
 492   neve = hMeasuredStraight->GetEntries();
 493   hMeasuredStraight->Scale(1.0/neve/binwidth);
 494   neve = hSimStraight->GetEntries();
 495   hSimStraight->Scale(1.0/neve/binwidth);
 496   //Gaus MC
 497   neve = hMeasuredGaus->GetEntries();
 498   hMeasuredGaus->Scale(1.0/neve/binwidth);
 499   neve = hSimGaus->GetEntries();
 500   hSimGaus->Scale(1.0/neve/binwidth);
 501   //DoubleExponential MC
 502   neve = hMeasuredDoubleExponential->GetEntries();
 503   hMeasuredDoubleExponential->Scale(1.0/neve/binwidth);
 504   neve = hSimDoubleExponential->GetEntries();
 505   hSimDoubleExponential->Scale(1.0/neve/binwidth);
 506   //Levy MC
 507   neve = hMeasuredLevy->GetEntries();
 508   hMeasuredLevy->Scale(1.0/neve/binwidth);
 509   neve = hSimLevy->GetEntries();
 510   hSimLevy->Scale(1.0/neve/binwidth);
 511   if(hITSData){
 512     neve = hITSData->GetEntries();
 513     hITSData->Scale(1.0/neve/binwidth);
 514     hITSData->Scale(1.0/hITSData->Integral());
 515   }
 516   if(hITSData){
 517     neve = hTPCData->GetEntries();
 518     hTPCData->Scale(1.0/neve/binwidth);
 519   }
 520   if(rescaleguess){
 521     neve = hTrainingReconstructed->GetEntries();
 522     hTrainingReconstructed->Scale(1.0/neve/binwidth);
 523     hTPCData->Scale(1.0/hTPCData->Integral());
 524   }
 525
 526   //TF1 *funcScale = new TF1("func","([0])/[1]*exp(-x/[1])",0,100);
 527 //   cout<<"integrals plain "<<hSim->Integral("width")<<" "<<hITS->Integral("width")<<" "<<hTPC->Integral("width")<<endl;
 528 //   cout<<"integrals plain "<<hSimReweighted->Integral("width")<<" "<<hMeasReweighted->Integral("width")<<endl;
 529 //   cout<<"integrals plain "<<hSimExponential->Integral("width")<<" "<<hMeasuredExponential->Integral("width")<<endl;
 530 //   cout<<"integrals plain "<<hSimFlat->Integral("width")<<" "<<hMeasuredFlat->Integral("width")<<endl;
 531 //   cout<<"integrals plain "<<hSimStraight->Integral("width")<<" "<<hMeasuredStraight->Integral("width")<<endl;
 532 //   cout<<"integrals plain "<<hSimGaus->Integral("width")<<" "<<hMeasuredGaus->Integral("width")<<endl;
 533 //   cout<<"integrals plain "<<hSimDoubleExponential->Integral("width")<<" "<<hMeasuredDoubleExponential->Integral("width")<<endl;
 534   //================SETTING TEST CASE========================
 535
 536   TH1F *hTruth;
 537   TH1F *hMeasured;
 538   switch(unfoldcase){
 539   case 0:
 540     cout<<"Test case is pure MC"<<endl;
 541     hTruth = hSim;
 542     if(its) hMeasured = hITS;
 543     else{hMeasured = hTPC;}
 544     break;
 545   case 1:
 546     cout<<"Test case is a reweighted MC"<<endl;
 547     hTruth = hSimReweighted;
 548     hMeasured = hMeasReweighted;
 549     break;
 550   case 2:
 551     cout<<"Test case is a exponential"<<endl;
 552     hTruth = hSimExponential;
 553     hMeasured = hMeasuredExponential;
 554     break;
 555   case 3:
 556     cout<<"Test case is a flat distribution"<<endl;
 557     hTruth = hSimFlat       ;
 558     hMeasured = hMeasuredFlat       ;
 559     break;
 560   case 4:
 561     cout<<"Test case is a straight line distribution"<<endl;
 562     hTruth = hSimStraight       ;
 563     hMeasured = hMeasuredStraight    ;
 564     break;
 565   case 5:
 566     cout<<"Test case is a half Gaussian distribution"<<endl;
 567     hTruth = hSimGaus       ;
 568     hMeasured = hMeasuredGaus    ;
 569     break;
 570   case 6:
 571     cout<<"Test case is a double exponential"<<endl;
 572     hTruth = hSimDoubleExponential;
 573     hMeasured = hMeasuredDoubleExponential;
 574     break;
 575   case 7:
 576     cout<<"Test case is a Levy distribution"<<endl;
 577     hTruth = hSimLevy       ;
 578     hMeasured = hMeasuredLevy       ;
 579     break;
 580   case 8:
 581     cout<<"Test case is a Power law"<<endl;
 582     hTruth = hSimPower       ;
 583     hMeasured = hMeasuredPower       ;
 584     break;
 585   }
 586   if(unfoldData){
 587     if(its) hMeasured = hITSData;
 588     else{hMeasured = hTPCData;}
 589   }
 590   //cout<<"MC truth: "<<hTruth->GetMean()<<" Mean observed "<<hMeasured->GetMean()<<endl;
 591   hTruth->GetXaxis()->SetTitle(hSim->GetTitle());
 592   hMeasured->GetXaxis()->SetTitle(hSim->GetTitle());
 593
 594
 595
 596
 597   //cout<<"Histo "<<histoname<<endl;
 598   TCanvas *canvasn1 = new TCanvas("canvasn1","canvas -1",600,600);
 599   canvasn1->SetTopMargin(0.020979);
 600   canvasn1->SetRightMargin(0.0184564);
 601   canvasn1->SetLeftMargin(0.0989933);
 602   canvasn1->SetBottomMargin(0.101399);
 603   canvasn1->SetBorderSize(0);
 604   canvasn1->SetFillColor(0);
 605   canvasn1->SetFillColor(0);
 606   canvasn1->SetBorderMode(0);
 607   canvasn1->SetFrameFillColor(0);
 608   canvasn1->SetFrameBorderMode(0);
 609   //canvasn1->Divide(2);
 610   canvasn1->SetLogy();
 611   hMeasured->SetMarkerStyle(22);
 612   hTrainingReconstructed->SetMarkerStyle(26);
 613   hTrainingReconstructed->SetMarkerColor(2);
 614   hTrainingReconstructed->SetLineColor(2);
 615   hMeasured->Draw();
 616   hTrainingReconstructed->Draw("same");
 617   TLegend *legendn1 = new TLegend(0.437919,0.800699,0.966443,0.958042);
 618   legendn1->AddEntry(hMeasured,"Measured");
 619   legendn1->AddEntry(hTrainingReconstructed,"Simulated reconstructed");
 620   legendn1->Draw();
 621   //cout<<"Reconstructed mean "<<hMeasured->GetMean()<<" Simulated reconstructed mean "<<hTrainingReconstructed->GetMean()<<endl;
 622   //return;
 623   if(unfoldData){
 624     float matchval = 10.0;
 625     int binsim = hTrainingReconstructed->FindBin(matchval+.01);
 626     float simval = hTrainingReconstructed->GetBinContent(binsim);
 627     //cerr<<__FILE__<<" "<<__LINE__<<endl;
 628     int binreco = hMeasured->FindBin(matchval+0.01);
 629     float recoval = hMeasured->GetBinContent(binreco);
 630     cout<<"matching at "<<matchval<<" simval "<<simval<<" recoval "<<recoval<<" bins "<<binsim<<" "<<binreco<<endl;
 631     if(recoval>0.0) hMeasured->Scale(simval/recoval);
 632     //else{cerr<<"Uh-oh!  Can't rescale.  :("<<endl;}
 633   }
 634   //   for(int i=1;i<=hMeasured->GetNbinsX();i++){
 635   //     cout<<hMeasured->GetBinContent(i);
 636   //     if(i!=hMeasured->GetNbinsX())cout<<", ";
 637   //   }
 638   //   cout<<endl;
 639   //hMeasured->Fit(funcPower,"N","",1,100);
 640   //funcLevy->Draw("same");
 641   //hMeasured->Fit(funcLevy,"","",1,100);
 642   //return;
 643   //return;
 644   if(simdataset==2009){
 645     cout<<"NOT doing full unfolding because this is 2009.  Exiting."<<endl;
 646     cout<<"Mean is "<<hTruth->GetMean()<<" ";
 647     GetChi2(hMeasured,hTrainingReconstructed);
 648     //return;
 649   }
 650
 651   //================UNFOLDING===========================
 652   if(zerolowetbins){
 653     for(int i=0;i<minbin;i++){
 654       hMeasured->SetBinContent(i,0.0);
 655     }
 656   }
 657   RooUnfoldBayes   unfold (&response, hMeasured, niter);    // OR
 658   unfold.SetSmoothing(true);
 659
 660   //================PLOTTING===========================
 661
 662
 663   TH1D* hReco1= (TH1D*) unfold.Hreco();
 664   //cout<<"MEAN WITHOUT ANY MASSAGING "<<hReco1->GetMean()<<endl;
 665
 666   TCanvas *canvas7 = new TCanvas("canvas7","Reconstructed Unfolded E_{T}",600,600);
 667   canvas7->SetTopMargin(0.020979);
 668   canvas7->SetRightMargin(0.0184564);
 669   canvas7->SetLeftMargin(0.0989933);
 670   canvas7->SetBottomMargin(0.101399);
 671   canvas7->SetBorderSize(0);
 672   canvas7->SetFillColor(0);
 673   canvas7->SetFillColor(0);
 674   canvas7->SetBorderMode(0);
 675   canvas7->SetFrameFillColor(0);
 676   canvas7->SetFrameBorderMode(0);
 677   hReco1->Draw();
 678
 679   TCanvas *canvas1 = new TCanvas("canvas1","Results",600,600);
 680   canvas1->SetTopMargin(0.020979);
 681   canvas1->SetRightMargin(0.0184564);
 682   canvas1->SetLeftMargin(0.0989933);
 683   canvas1->SetBottomMargin(0.101399);
 684   canvas1->SetBorderSize(0);
 685   canvas1->SetFillColor(0);
 686   canvas1->SetFillColor(0);
 687   canvas1->SetBorderMode(0);
 688   canvas1->SetFrameFillColor(0);
 689   canvas1->SetFrameBorderMode(0);
 690   //canvas1->Divide(2);
 691   canvas1->SetLogy();
 692
 693   //rescale just in case
 694   //   float myscale = hTruth->Integral("width") / hReco1->Integral("width");
 695   //   cout<<"Rescaling result by "<<myscale<<endl;
 696   //   hReco1->Scale(myscale);
 697   hReco1->Fit(func,"NI","",0.0,1.0);
 698   funcDoubleExponential->SetParameter(1,func->GetParameter(1));
 699   hReco1->Fit(funcDoubleExponential,"NI","",0.0,1.0);
 700   float xgoal1 = 1.0;
 701   float xgoal2 = 0.01;
 702   if(recodataset ==2010){
 703     float xgoal1 = 1.0;
 704     float xgoal2 = 2.0;
 705   }
 706   float x1 = hReco1->GetBinCenter(hReco1->FindBin(xgoal1));
 707   float y1 = hReco1->GetBinContent(hReco1->FindBin(xgoal1));
 708   float x2 = hReco1->GetBinCenter(hReco1->FindBin(xgoal2));
 709   float y2 = hReco1->GetBinContent(hReco1->FindBin(xgoal2));
 710 //   float x2 = hReco1->GetBinCenter(1);
 711 //   float y2 = hReco1->GetBinContent(1);
 712   float myb = -TMath::Log(y1/y2)/(x1-x2);
 713   float myA = y1/myb/exp(-myb*x1);
 714   func->SetParameter(0,myA);
 715   func->SetParameter(1,1.0/myb);
 716   //cout<<"x1 "<<x1<<" y1 "<<y1<<" x2 "<<x2<<" y2 "<<y2<<" b "<<1/myb<<" A "<<myA/myb<<endl;
 717
 718
 719   hTruth->SetMarkerStyle(22);
 720   hReco1->SetMarkerStyle(30);
 721   hReco1->SetLineColor(2);
 722   hReco1->SetMarkerColor(2);
 723   hITS->SetMarkerColor(3);
 724   hITS->SetLineColor(3);
 725   hTPC->SetMarkerColor(7);
 726   hTPC->SetLineColor(7);
 727   //hTruth->GetXaxis()->SetRange(0,hTruth->FindBin(10.0));
 728   hReco1->Draw();
 729   if(!unfoldData) hTruth->Draw("same");
 730   //hITS->Draw("same");
 731   //hTPC->Draw("same");
 732   func->Draw("same");
 733   funcDoubleExponential->Draw("same");
 734   hTruth->GetXaxis()->SetRange(0,hTruth->GetNbinsX());
 735   //cout<<" Means "<<hTruth->GetMean()<<" "<<hReco1->GetMean()<<endl;
 736 //   hTruth->GetXaxis()->SetRange(hTruth->FindBin(1.),hTruth->GetNbinsX());
 737 //   hReco1->GetXaxis()->SetRange(hTruth->FindBin(1.),hReco1->GetNbinsX());
 738   hTruth->GetXaxis()->SetRange(0,hTruth->FindBin(20.));
 739   hReco1->GetXaxis()->SetRange(0,hTruth->FindBin(20.));
 740   //cout<<" Truncated Means "<<hTruth->GetMean()<<" "<<hReco1->GetMean()<<endl;
 741   hTruth->GetXaxis()->SetRange(0,hTruth->GetNbinsX());
 742   hReco1->GetXaxis()->SetRange(0,hReco1->GetNbinsX());
 743
 744   TLegend *leg = new TLegend(0.505034,0.744755,0.605705,0.952797);
 745   leg->SetFillStyle(0);
 746   leg->SetFillColor(0);
 747   leg->SetBorderSize(0);
 748   //leg->AddEntry(hSim,"Simulated E_{T}");
 749   leg->AddEntry(hReco1,"Unfolded result");
 750   leg->AddEntry(func,"Exponential");
 751   leg->AddEntry(funcDoubleExponential,"Double exponential");
 752   hSim->SetLineColor(1);
 753   hSim->GetYaxis()->SetTitleOffset(1.3);
 754   hSim->GetYaxis()->SetTitle("1/N_{eve} dN/dE_{T}");
 755   leg->SetTextSize(0.0437063);
 756   leg->Draw();
 757   canvas1->SaveAs(Form("%s%s.png",outputfilename,"Extrap"));
 758   canvas1->SaveAs(Form("%s%s.eps",outputfilename,"Extrap"));
 759   canvas1->SaveAs(Form("%s%s.C",outputfilename,"Extrap"));
 760
 761   //canvas1->cd(2);
 762   TCanvas *canvas6 = new TCanvas("canvas6","Unfolded data/simulated MC truth",600,600);
 763   canvas6->SetTopMargin(0.020979);
 764   canvas6->SetRightMargin(0.0184564);
 765   canvas6->SetLeftMargin(0.0989933);
 766   canvas6->SetBottomMargin(0.101399);
 767   canvas6->SetBorderSize(0);
 768   canvas6->SetFillColor(0);
 769   canvas6->SetFillColor(0);
 770   canvas6->SetBorderMode(0);
 771   canvas6->SetFrameFillColor(0);
 772   canvas6->SetFrameBorderMode(0);
 773
 774   TH1D *hReco1Clone = (TH1D*) hReco1->Clone("hReco1Clone");
 775   if(unfoldData){
 776     float matchval = 1.0;
 777     int binsim = hTruth->FindBin(matchval+.01);
 778     float simval = hTruth->GetBinContent(binsim);
 779     //cerr<<__FILE__<<" "<<__LINE__<<endl;
 780     int binreco = hReco1Clone->FindBin(matchval+0.01);
 781     float recoval = hReco1Clone->GetBinContent(binreco);
 782     //cout<<"matching at "<<matchval<<" simval "<<simval<<" recoval "<<recoval<<" bins "<<binsim<<" "<<binreco<<endl;
 783     if(recoval>0.0) hReco1Clone->Scale(simval/recoval);
 784   }
 785   hReco1Clone->Divide(hTruth);
 786   hReco1Clone->Draw("");
 787   hReco1Clone->GetYaxis()->SetTitle("N_{eve}^{reco}/N_{eve}^{true}");
 788   hReco1Clone->GetXaxis()->SetTitle("E_{T}");
 789   hReco1Clone->SetMinimum(0.0);
 790   hReco1Clone->SetMaximum(2.0);
 791   hReco1Clone->GetXaxis()->SetRange(1,hReco1Clone->FindBin(15.0));
 792   //   canvas1->cd(3);
 793
 794   //   hReco1->Draw("");
 795   //cout<<"Bin 1: true :"<< hSim->GetBinContent(1)<<"+/-"<<hSim->GetBinError(1)<<" reco:"<<hReco1->GetBinContent(1)<<"+/-"<<hSim->GetBinError(1)<<endl;
 796
 797   //  unfold.PrintTable(cout,hSim);
 798
 799   if(trainingcase==1 || unfoldcase==1){
 800     TCanvas *canvas2 = new TCanvas("WeightingFunction","WeightingFunction",600,600);
 801     canvas2->SetTopMargin(0.020979);
 802     canvas2->SetRightMargin(0.0184564);
 803     canvas2->SetLeftMargin(0.0989933);
 804     canvas2->SetBottomMargin(0.101399);
 805     canvas2->SetBorderSize(0);
 806     canvas2->SetFillColor(0);
 807     canvas2->SetFillColor(0);
 808     canvas2->SetBorderMode(0);
 809     canvas2->SetFrameFillColor(0);
 810     canvas2->SetFrameBorderMode(0);
 811     fWeight->Draw();
 812   }
 813   //   TCanvas *canvas3 = new TCanvas("canvas3","canvas3",1200,600);
 814   //   canvas3->SetTopMargin(0.020979);
 815   //   canvas3->SetRightMargin(0.0184564);
 816   //   canvas3->SetLeftMargin(0.0989933);
 817   //   canvas3->SetBottomMargin(0.101399);
 818   //   canvas3->SetBorderSize(0);
 819   //   canvas3->SetFillColor(0);
 820   //   canvas3->SetFillColor(0);
 821   //   canvas3->SetBorderMode(0);
 822   //   canvas3->SetFrameFillColor(0);
 823   //   canvas3->SetFrameBorderMode(0);
 824   //   hReco1->Draw();
 825 //   TCanvas *canvas4 = new TCanvas("canvas4","canvas4",1200,600);
 826 //   canvas4->SetTopMargin(0.020979);
 827 //   canvas4->SetRightMargin(0.0184564);
 828 //   canvas4->SetLeftMargin(0.0989933);
 829 //   canvas4->SetBottomMargin(0.101399);
 830 //   canvas4->SetBorderSize(0);
 831 //   canvas4->SetFillColor(0);
 832 //   canvas4->SetFillColor(0);
 833 //   canvas4->SetBorderMode(0);
 834 //   canvas4->SetFrameFillColor(0);
 835 //   canvas4->SetFrameBorderMode(0);
 836 //   hTrainingReconstructed->Draw();
 837 //   hTrainingTruth->Draw("same");
 838 //   hTrainingTruth->SetLineColor(1);
 839 //   hTrainingReconstructed->SetLineColor(2);
 840
 841
 842   //   TCanvas *canvas5 = new TCanvas("canvas5","canvas5",1200,600);
 843   //   canvas5->SetTopMargin(0.020979);
 844   //   canvas5->SetRightMargin(0.0184564);
 845   //   canvas5->SetLeftMargin(0.0989933);
 846   //   canvas5->SetBottomMargin(0.101399);
 847   //   canvas5->SetBorderSize(0);
 848   //   canvas5->SetFillColor(0);
 849   //   canvas5->SetFillColor(0);
 850   //   canvas5->SetBorderMode(0);
 851   //   canvas5->SetFrameFillColor(0);
 852   //   canvas5->SetFrameBorderMode(0);
 853
 854   //Calculating mean/mean error with covariant matrix
 855   float cutoffLow = hReco1->GetBinLowEdge(hReco1->FindBin(1.0+0.01));
 856   float cutoffHigh = hReco1->GetBinLowEdge(hReco1->FindBin(30+0.01));
 857   if(had==1) cutoffHigh = hReco1->GetBinLowEdge(hReco1->FindBin(20+0.01));
 858   int cutoffbin = hReco1->FindBin(1.01);//Get the bin at 1 GeV, with a little wiggle room to be sure it returns the bin I mean
 859   int cutoffbinHigh = hReco1->FindBin(30.01);//Get the bin at 1 GeV, with a little wiggle room to be sure it returns the bin I mean
 860   if(had==1) cutoffbinHigh = hReco1->FindBin(20.01);//Get the bin at 1 GeV, with a little wiggle room to be sure it returns the bin I mean
 861   if(had==2) cutoffbinHigh = hReco1->FindBin(17.01);//Get the bin at 1 GeV, with a little wiggle room to be sure it returns the bin I mean
 862   int nbins = hReco1->GetNbinsX();
 863   TMatrixD cov = unfold.Ereco();
 864   TVectorD result = unfold.Vreco();
 865   TVectorD resultError = unfold.ErecoV();
 866   //   TMatrixD cov = unfold.GetMeasuredCov();
 867   //   TVectorD result = unfold.Vmeasured();
 868   //   TVectorD resultError = unfold.Emeasured();
 869   float cutoff = hReco1->GetBinLowEdge(cutoffbin);
 870   float binwidth = hReco1->GetBinWidth(1);
 871   float mean = 0;
 872   float meanerr = 0;//this is actually the error squared but I'm just trying a different way of calculating it
 873   float meanvar = 0;
 874   float total = 0;
 875   //truncated values
 876   float meantrunc = 0;
 877   float meanerrtrunc = 0;//this is actually the error squared but I'm just trying a different way of calculating it
 878   float meanvartrunc = 0;
 879   float totaltrunc = 0;
 880   for(int i=0;i<cutoffbinHigh;i++){
 881     total += result(i)*binwidth;
 882     float energyi = hReco1->GetBinCenter(i+1);
 883     mean += result(i)*binwidth*energyi;
 884     meanerr += resultError(i)*resultError(i)*binwidth*energyi*binwidth*energyi;
 885     if(i>=cutoffbin){
 886       totaltrunc += result(i)*binwidth;
 887       meantrunc += result(i)*binwidth*energyi;
 888       meanerrtrunc += resultError(i)*resultError(i)*binwidth*energyi*binwidth*energyi;
 889     }
 890     for(int j=0;j<nbins;j++){
 891       float energyj = hReco1->GetBinCenter(j+1);
 892       meanvar += energyi*binwidth  *  energyj*binwidth  * cov(i,j);
 893       if(i>=cutoffbin && j>=cutoffbin){
 894         meanvartrunc += energyi*binwidth  *  energyj*binwidth  * cov(i,j);
 895       }
 896       //       if(i==j && i<40 && j<40 && cov(i,j)>0.0){
 897       //        cout<<"i "<<i<<" bin center "<<hReco1->GetBinCenter(i+1)<<" j "<<j<<" cov "<<cov(i,j)<<" err "<< result(i)*binwidth*energyi  *  result(j)*binwidth*energyj  * cov(i,j);
 898       //        if(i==j) cout<<" error "<<resultError(i)<<" rel. err "<< resultError(i)/result(i)*100.0<<"%";
 899       //        cout<< endl;
 900       //       }
 901     }
 902     //     cout<<"i "<<i;
 903     //     cout<<" result "<<result(i);
 904     //     cout<<" +/- "<<resultError(i);
 905     //     cout<<endl;
 906   }
 907   cout<<"Mean "<<mean<<" total "<<total<<" mean err ";
 908   if(meanvar>0) cout<<TMath::Sqrt(meanvar);
 909   cout<<" or ";
 910   if(meanerr>0)cout<<TMath::Sqrt(meanerr);
 911   cout<<endl;
 912   cout<<"Mean trunc";
 913   if(totaltrunc>0)cout<<meantrunc/totaltrunc;
 914   cout<<endl;
 915   TF1 *funcMean = new TF1("funcMean","x*([0])/[1]*exp(-x/[1])",0,100);
 916   for(int i=0; i<2;i++){funcMean->SetParameter(i,func->GetParameter(i));}
 917   float expoExtrap = -1;
 918   if(funcMean->GetParameter(1)!=0) funcMean->Integral(0.0,cutoffLow);
 919   float expoExtrapTotal = -1;
 920   if(func->GetParameter(1)!=0) func->Integral(0.0,cutoffLow);
 921   TF1 *funcDoubleExponentialMean = new TF1("funcDoubleExponentialMean","x*[0]/[1]*exp(-x/[1])+x*[2]/[3]*exp(-x/[3])",0,100);
 922   for(int i=0; i<4;i++){funcDoubleExponentialMean->SetParameter(i,funcDoubleExponential->GetParameter(i));}
 923   float dblExpoExtrap = funcDoubleExponentialMean->Integral(0.0,cutoffLow);
 924   float dblExpoExtrapTotal = funcDoubleExponential->Integral(0.0,cutoffLow);
 925   float total = hReco1->Integral(1,hReco1->FindBin(cutoffHigh-0.01),"width");
 926   float truncated = hReco1->Integral(hReco1->FindBin(cutoffLow+0.01),hReco1->FindBin(cutoffHigh-0.01),"width");
 927   hReco1->GetXaxis()->SetRange(1,hReco1->FindBin(cutoffHigh));
 928   hReco1->GetXaxis()->SetRange(hReco1->FindBin(cutoffLow),hReco1->FindBin(cutoffHigh));
 929   cout<<"Truncated Mean "<<hReco1->GetMean()<<" total "<<total<<" mean err ";
 930   if(meanvar>0) cout<<TMath::Sqrt(meanvar);
 931   cout<<endl;
 932   hReco1->GetXaxis()->SetRange(1,hReco1->GetNbinsX());
 933   cout<<"Extrapolations:  exponential "<<expoExtrap<<" double exponential "<<dblExpoExtrap<<endl;
 934   cout<<"Total: high ";
 935   if((totaltrunc+dblExpoExtrapTotal)>0) cout<<(dblExpoExtrap+meantrunc)/(totaltrunc+dblExpoExtrapTotal);
 936   cout<<" best ";
 937   if((totaltrunc+expoExtrapTotal)>0) cout<<(expoExtrap+meantrunc)/(totaltrunc+expoExtrapTotal);
 938   cout<<" low "<<mean<<endl;
 939   cout<<"Raw mean "<<hReco1->GetMean();
 940   cout<<" Total scaled: low ";
 941   if((totaltrunc+dblExpoExtrapTotal)>0) cout<<(dblExpoExtrap+meantrunc)/(totaltrunc+dblExpoExtrapTotal);
 942   cout<<" best ";
 943   if(totaltrunc+expoExtrapTotal>0) cout<<(expoExtrap+meantrunc)/(totaltrunc+expoExtrapTotal);
 944   cout<<" high ";
 945   if(total>0) cout<<mean/total;
 946   cout<<endl;
 947   canvas1->cd();
 948   funcDoubleExponential->SetLineStyle(2);
 949   func->Draw("same");
 950   funcDoubleExponential->Draw("same");
 951   canvas1->SaveAs(Form("%s.C",outputfilename));
 952   canvas1->SaveAs(Form("%s.eps",outputfilename));
 953   canvas1->SaveAs(Form("%s.png",outputfilename));
 954
 955
 956
 957   TCanvas *canvas5 = new TCanvas("RecovsSimReco","RecovsSimReco",600,600);
 958   canvas5->SetTopMargin(0.020979);
 959   canvas5->SetRightMargin(0.0184564);
 960   canvas5->SetLeftMargin(0.0989933);
 961   canvas5->SetBottomMargin(0.101399);
 962   canvas5->SetBorderSize(0);
 963   canvas5->SetFillColor(0);
 964   canvas5->SetFillColor(0);
 965   canvas5->SetBorderMode(0);
 966   canvas5->SetFrameFillColor(0);
 967   canvas5->SetFrameBorderMode(0);
 968   canvas5->SetLogy();
 969   hMeasured->Draw();
 970   nevents = neveUsed;
 971   TH1F *hMeasuredSim = GetReconstructedEt((TH1*)hReco1,(int) nevents,(char *)"SmearedTruth",(float) hReco1->GetBinLowEdge(1),(float) hReco1->GetBinLowEdge(hReco1->GetNbinsX()),(int) hReco1->GetNbinsX());//(TH1D*) unfold.Hmeasured();
 972   float myscale2 = hMeasured->GetBinContent(hMeasured->FindBin(10.0)) / hMeasuredSim->GetBinContent(hMeasuredSim->FindBin(10.0));
 973   cout<<"Scaling by "<<hMeasured->GetBinContent(hMeasured->FindBin(10.0))<<" / "<<hMeasuredSim->GetBinContent(hMeasuredSim->FindBin(10.0))<<" = "<<myscale2<<endl;
 974   GetChi2(hMeasured,hMeasuredSim);
 975   //cout<<"Chi^2 of Smeared compared to measured "<<GetChi2(hMeasured,hMeasuredSim)<<endl;
 976   hMeasuredSim->Scale(myscale2);
 977   hMeasuredSim->Draw("same");
 978   hMeasuredSim->SetMarkerStyle(21);
 979   hMeasuredSim->SetLineColor(2);
 980   hMeasuredSim->SetMarkerColor(2);
 981   TLegend *legend5 = new TLegend(0.437919,0.800699,0.966443,0.958042);
 982   legend5->AddEntry(hMeasured,"Measured");
 983   legend5->AddEntry(hMeasuredSim,"Simulated Measured");
 984   legend5->Draw();
 985   canvas5->SaveAs(Form("%s%s.png",outputfilename,"Reco"));
 986   canvas5->SaveAs(Form("%s%s.eps",outputfilename,"Reco"));
 987   canvas5->SaveAs(Form("%s%s.C",outputfilename,"Reco"));
 988
 989 //   func->Draw();
 990   //funcDoubleExponential->Draw();
 991 }
 992
 993 TH1F *GetReconstructedEt(TF1 *inputfunc, int nevents, char *name, float lowbound, float highbound, int nbins,bool smooth){
 994   TH1F *hResult = new TH1F(name,ytitle->Data(),nbins,lowbound,highbound);
 995   hResult->Sumw2();
 996   //cerr<<__FILE__<<" "<<__LINE__<<" Creating "<<hResult->GetName()<<" with "<<nevents<<" events"<<endl;
 997   for(int i=0;i<nevents;i++){
 998     float et = inputfunc->GetRandom(lowbound,highbound);
 999     int mybin = resolutionFull->GetXaxis()->FindBin(et);
1000     if(mybin>0 && mybin<=nbins){//then this is in our range...
1001       float myres = ((TH1D*)histoarray.At(mybin-1))->GetRandom();
1002       //float etreco = (1-myres)*et;
1003       float etreco = myres;
1004       //cout<<"et true "<<et<<" reco "<<etreco<<endl;
1005       hResult->Fill(etreco);
1006     }
1007     //if(i%100000==0) Smooth(hResult,inputfunc);
1008   }
1009   //if(smooth)
1010 Smooth(hResult,inputfunc);
1011   //float binwidth = hResult->GetBinWidth(1);
1012   //hResult->Scale(1.0/nevents/binwidth);
1013   return hResult;
1014 }
1015 TH1F *GetReconstructedEt(TH1 *input, int nevents, char *name, float lowbound, float highbound, int nbins){
1016
1017   TH1F *hResult = new TH1F(name,ytitle->Data(),nbins,lowbound,highbound);
1018   hResult->Sumw2();
1019   for(int i=0;i<nevents;i++){
1020     float et = input->GetRandom();
1021     int mybin = resolutionFull->GetXaxis()->FindBin(et);
1022     if(mybin>0 && mybin<=nbins){//then this is in our range...
1023       float myres = ((TH1D*)histoarray.At(mybin-1))->GetRandom();
1024       //float etreco = (1-myres)*et;
1025       float etreco = myres;
1026       //cout<<"et true "<<et<<" reco "<<etreco<<endl;
1027       hResult->Fill(etreco);
1028     }
1029   }
1030   hResult->GetYaxis()->SetTitle("1/N_{eve} dN/dE_{T}");
1031   hResult->GetXaxis()->SetTitle("E_{T}");
1032   //float binwidth = hResult->GetBinWidth(1);
1033   //hResult->Scale(1.0/nevents/binwidth);
1034   return hResult;
1035 }
1036
1037 TH1F *GetTrueEt(TF1 *inputfunc, int nevents, char *name, float lowbound, float highbound, int nbins){
1038
1039   //cerr<<__FILE__<<" "<<__LINE__<<" Throwing simulation "<<name<<" with "<<nevents<<" events "<<endl;
1040   TH1F *hResult = new TH1F(name,ytitle->Data(),nbins,lowbound,highbound);
1041   hResult->Sumw2();
1042   //cerr<<"event ";
1043   for(int i=0;i<nevents;i++){
1044     //if(i%1000) cerr<<i<<" ";
1045     float et = inputfunc->GetRandom(lowbound,highbound);
1046     hResult->Fill(et);
1047   }
1048   //cerr<<endl;
1049   float binwidth = hResult->GetBinWidth(1);
1050   //hResult->Scale(1.0/nevents/binwidth);
1051   hResult->GetYaxis()->SetTitle("1/N_{eve} dN/dE_{T}");
1052   hResult->GetXaxis()->SetTitle("E_{T}");
1053   return hResult;
1054 }
1055
1056
1057 Float_t GetChi2(TH1F *reconstructed, TH1F *simulated){
1058   Float_t chi2 = 0;
1059   int nbins = reconstructed->FindBin(15.0);
1060   int lowbin =  reconstructed->FindBin(1.01);
1061   float ndf=0;
1062   for(int i=4;i<=nbins;i++){
1063     float y = reconstructed->GetBinContent(i);
1064     float yerr = reconstructed->GetBinError(i);
1065     float f = simulated->GetBinContent(i);
1066     float ferr = simulated->GetBinError(i);
1067     if((yerr*yerr)>0){
1068       //       cout<<"i "<< TMath::Power(y-f,2)/(yerr*yerr+ferr*ferr)<<endl;
1069       //       chi2+= TMath::Power(y-f,2)/(yerr*yerr+ferr*ferr);
1070       //cout<<"i "<< TMath::Power(y-f,2)/(yerr*yerr)<<endl;
1071       chi2+= TMath::Power(y-f,2)/(yerr*yerr);
1072       //float relerry = yerr/y;
1073       //float relerrf = ferr/f;
1074       //cout<<"i "<<i<<" chi2 increment "<< TMath::Power(y-f,2)/(yerr*yerr)<<" y "<<y<<" f "<<f<<" yerr/y "<<relerry<<" ferr/f "<<relerrf<<endl;
1075       ndf++;
1076     }
1077   }
1078   ndf--;//There is one parameter
1079   cout<<"Chi^2/ndf "<<chi2<<"/"<<ndf<<" = "<<chi2/ndf<<endl;
1080   if(ndf>0)return chi2/ndf;
1081   else{return 0;}
1082 }
1083
1084 TH1F *RestrictAxes(TH1F *old,float minEt,float maxEt){
1085   int minbin = old->FindBin(minEt+0.01);
1086   int maxbin = old->FindBin(maxEt-0.01);
1087   //cout<<"Restricting range.  Old bin width:"<<old->GetBinWidth(1);
1088   int totnbins = maxbin-minbin+1;
1089   TH1F *hNew = new TH1F(Form("%sNew",old->GetName()),old->GetTitle(),totnbins,minEt,maxEt);
1090   hNew->GetYaxis()->SetTitle(old->GetYaxis()->GetTitle());
1091   hNew->GetXaxis()->SetTitle(old->GetXaxis()->GetTitle());
1092   //cout<<old->GetName()<<" New range "<<minEt<<"-"<<maxEt<<" bins "<<minbin<<"-"<<maxbin<<endl;
1093   for(int i=minbin;i<=maxbin;i++){
1094     //cout<<hNew->FindBin(old->GetBinCenter(i))<<":"<<old->GetBinContent(i);
1095     //if(i!=maxbin) cout<<",";
1096     int oldbin = i;
1097     float newbin = hNew->FindBin(old->GetBinCenter(oldbin));
1098     //cout<<"old "<<oldbin<<" new "<<newbin<<" center "<<hNew->GetBinCenter(newbin)<<"="<<old->GetBinCenter(oldbin)<<" content "<<old->GetBinContent(oldbin)<<" +/- "<<old->GetBinError(oldbin)<<endl;
1099     hNew->SetBinContent(newbin,old->GetBinContent(oldbin));
1100     hNew->SetBinError(newbin,old->GetBinError(oldbin));
1101   }
1102   //cout<<endl;
1103   hNew->Sumw2();
1104   //delete old;
1105   //old = hNew;
1106   //cout<<" new bin width "<<hNew->GetBinWidth(1)<<endl;
1107   return hNew;
1108 }
1109 TH2F *RestrictAxes(TH2F *old,float minEt,float maxEt){
1110   int minbin = old->GetYaxis()->FindBin(minEt+0.01);
1111   int maxbin = old->GetYaxis()->FindBin(maxEt-0.01);
1112   //cout<<"Restricting range.  Old bin width:"<<old->GetBinWidth(1);
1113   int totnbins = maxbin-minbin+1;
1114   TH2F *hNew = new TH2F(Form("%sNew",old->GetName()),old->GetTitle(),totnbins,minEt,maxEt,totnbins,minEt,maxEt);
1115   hNew->GetYaxis()->SetTitle(old->GetYaxis()->GetTitle());
1116   hNew->GetXaxis()->SetTitle(old->GetXaxis()->GetTitle());
1117   for(int i=minbin;i<=maxbin;i++){
1118     int newbinx = hNew->GetXaxis()->FindBin(old->GetXaxis()->GetBinCenter(i));
1119     int oldbinx = i;//hNew->GetXaxis()->FindBin(old->GetXaxis()->GetBinCenter(i));
1120     for(int j=minbin;j<=maxbin;j++){
1121       int newbiny = hNew->GetYaxis()->FindBin(old->GetYaxis()->GetBinCenter(j));
1122       int oldbiny = j;//hNew->GetYaxis()->FindBin(old->GetYaxis()->GetBinCenter(j));
1123       //int mybinNew = hNew->FindBin(old->GetXaxis()->GetBinCenter(i),old->GetYaxis()->GetBinCenter(j));
1124       //int mybinOld = old->FindBin(old->GetXaxis()->GetBinCenter(i),old->GetYaxis()->GetBinCenter(j));
1125       //cout<<"i "<<i<<" j "<<j<<" ";
1126       //cout<<"Old bin ("<<oldbinx<<","<<oldbiny<<") is becoming ("<<newbinx<<","<<newbiny<<") with content "<<old->GetBinContent(oldbinx,oldbiny)<<" +/- "<<old->GetBinError(oldbinx,oldbiny)<<endl;
1127       hNew->SetBinContent(newbinx,newbiny,old->GetBinContent(oldbinx,oldbiny));
1128       hNew->SetBinError(newbinx,newbiny,old->GetBinError(oldbinx,oldbiny));
1129       //hNew->SetBinError(mybinNew,old->GetBinError(mybinOld));
1130     }
1131   }
1132   hNew->Sumw2();
1133   return hNew;
1134   //delete old;
1135   //old = hNew;
1136   //cout<<" new bin width "<<hNew->GetBinWidth(1)<<endl;
1137 }
1138
1139 void Smooth(TH1F *histo,TF1 *func){
1140   cout<<"Smoothing..."<<endl;
1141   //histo->Smooth(1,"R");
1142   histo->Fit(func,"N","",15,35);
1143   for(int i=histo->FindBin(15.0);i<histo->GetNbinsX();i++){
1144     float thisval = histo->GetBinContent(i);
1145     float thiscenter = histo->GetBinCenter(i);
1146     float thiserr = histo->GetBinError(i);
1147     float expectation = func->Eval(thiscenter);
1148     //if the number is more than 10 sigma away from the expectation and we would have expected at least 1 entry in that bin, recalculate it assuming Gaussian smoothing
1149     if(thisval<1.0) break;
1150     if(expectation>1.0 && TMath::Abs(thisval-expectation)/thiserr>5){
1151       float nSigma = myGaussian->GetRandom();
1152       histo->SetBinContent(i,expectation+nSigma*TMath::Sqrt(expectation));
1153       //cout<<"Resetting bin "<<i<<" at "<<histo->GetBinCenter(i)<<" from "<<thisval<<" which was "<<TMath::Abs(thisval-expectation)/thiserr<<" sigma away to "<<expectation+nSigma*TMath::Sqrt(expectation)<<" nSigma "<<nSigma<<" away from "<<expectation<<endl;
1154     }
1155 //     float nextval = histo->GetBinContent(i+1);
1156 //     float lastval = histo->GetBinContent(i-1);
1157 //     float nextcenter = histo->GetBinCenter(i+1);
1158 //     float lastcenter = histo->GetBinCenter(i-1);
1159 //     float nexterr = histo->GetBinError(i+1);
1160 //     float lasterr = histo->GetBinError(i-1);
1161 //     if(nextval>thisval){//We have a problem!  This function should decrease!
1162 //       //impose y=A e(-x/b) and use previous two values to calculate b and A
1163 //       float b = TMath::Log(thisval/lastval)/(lastcenter-thiscenter);
1164 //       float A = thisval/exp(-b*thiscenter);
1165 //       float newval = A*exp(-b*nextcenter);
1166 //       cout<<"A "<<A<<" b "<<" center "<<nextcenter<<" old "<<histo->GetBinContent(i+1)<<" new "<<newval<<endl;
1167 //       //if(newval)histo->SetBinContent(i+1,newval);
1168 //     }
1169   }
1170 }