//std libs:
#include <vector>
#include <typeinfo>
//#include <cxxabi.h>//!!! VIP this is the one we need!!!
#include <iostream>
#include <typeinfo>
#include <cstdlib>
#include "math.h"
#include <map>
//#include <unordered_map>
#include <string>
#include <utility>

//ROOT
#include "TH1F.h"
#include "TH2F.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TKey.h"
#include "TCollection.h"
#include "TVectorT.h"
#include "TString.h"
#include "TSystemDirectory.h"
#include "TSystemFile.h"
#include "TSystem.h"
//#include "TIter.h"
//Unfolding:
#include "TRandom.h"
#include "TH1D.h"
#include "RooUnfoldResponse.h"
#include "RooUnfoldBayes.h"
//preventing quitting:
#include "TApplication.h"
#include "TH1D.h"
//poisson:
#include "TRandom.h"
#include "TStyle.h"
#include "TMatrix.h"
#include "TMatrixD.h"


void PutDirContentinVector(TString Directory, std::vector<TString>& FileNameVector);
TString* GenerateTreeNames(TFile* f, int& nkeys);

//void InitializeParameters(int argc,char **argv, int & DataorMonteCarlo, int & dopAorAp, int & CastorPosition, int & SelectedTriggerBit,int & TriggerBitForEfficiency, double & CastorRecoJetEnergyCutoff,int & BinningKinematicScenario, int & lowerbound,int & upperbound,int & NBins, double & CastorGenJetEnergyCutoff, int & UnfoldingScenario, double & genjetetalow, double & genjetetahigh, double & deltaphicut, int & NumberFilesToRead, int & nev);

//init in order given in argv
void InitializeParameters(int argc,char **argv, int & DataorMonteCarlo, int & dopAorAp, int & SelectedTriggerBit,int & TriggerBitForEfficiency, double & CastorRecoJetEnergyCutoff, int & CastorPosition, int & UnfoldingScenario, int & BinningKinematicScenario, int & NumberFilesToRead, int & nev, int & lowerbound,int & upperbound,int & NBins, double & CastorGenJetEnergyCutoff, double & genjetetalow, double & genjetetahigh, double & deltaphicut);

  
TString BuildNTupleDescriptor(int DataorMonteCarlo, int dopAorAp, int CastorPosition);
//TString BuildAnalysisDescriptor(int SelectedTriggerBit, int CastorRecoJetEnergyCutoff, int lowerbound,int NBins, int CastorGenJetEnergyCutoff, int genjetetalow,int genjetetahigh, int deltaphicut, int DataorMonteCarlo);
TString BuildAnalysisDescriptor(int SelectedTriggerBit, double CastorRecoJetEnergyCutoff, int lowerbound,int NBins, double CastorGenJetEnergyCutoff, double genjetetalow, double genjetetahigh, double deltaphicut, int DataorMonteCarlo);

void PerturbResponseHisto(TH2F* ResponseHisto);
void Perturb1DHisto(TH1* ResponseHisto);
void FillMatrixV(TMatrixD & V, TH1* mean , TH1* Perturbed);

using namespace std;

//int main(int argc, char **argv){
//int UnfoldMC_RedChiBackSmeared(int argc, char **argv){
int UnfoldMC_RedChiBackSmeared(){

  int argc=1;
  char **argv;
  //argv[0]="";
  
  int NBayesIts=11;
  int NPerturbIterations=10;
  //  int NPerturbIterations=1;
  int NumberGoodEventsSpectrum;
  int NumberGoodEventsResponseObject;
  double RatioGoodEvtsResponseObjectandSpectrum;

  gStyle->SetOptStat(0);
  cout<<"start UnfoldMC"<<endl;
  // TApplication appl("appl",NULL,NULL); //take out since it seems to make memory issues..

  //a 100.0000 Evts sample, correct binning with 2bins less that were empty for 100K evts:
  //  TFile oldf("../Exec_ResponseMatrix_WithRooUnf/OutputFilesResponseObject/ReggeGribbov/2016_1_21_AlgoTrigger35_generator1_CollsionTypeAp_RecoJetEnergyCut150_GenJetEnergyCut150_genjetetalow-52_genjetetahigh-66_deltaphicut0.5_NumberOfFiles1_UpperBoundHisto2475_NBins31.root","READ");

  /*old root files:
  //a 100.0000 Evts sample, correct binning with 2bins less that were empty for 100K evts:
  TFile f("../Exec_ResponseMatrix_WithRooUnf/OutputFilesResponseObject/ReggeGribbov/2016_1_21_AlgoTrigger35_generator1_CollsionTypeAp_RecoJetEnergyCut150_GenJetEnergyCut150_genjetetalow-52_genjetetahigh-66_deltaphicut0.5_NumberOfFiles1_UpperBoundHisto2475_NBins31.root","READ");
  //4 bins only
  TFile f("../Exec_ResponseMatrix_WithRooUnf/OutputFilesResponseObject/ReggeGribbov/2016_1_20_AlgoTrigger35_generator1_CollsionTypeAp_RecoJetEnergyCut150_GenJetEnergyCut150_genjetetalow-52_genjetetahigh-66_deltaphicut0.5_NumberOfFiles1.root","READ");
  //2 bins only
  TFile f("../Exec_ResponseMatrix_WithRooUnf/OutputFilesResponseObject/ReggeGribbov/2016_1_22_AlgoTrigger35_generator1_CollsionTypeAp_RecoJetEnergyCut150_GenJetEnergyCut150_genjetetalow-52_genjetetahigh-66_deltaphicut0.5_NumberOfFiles1_UpperBoundHisto450_NBins2.root","READ");
  */
    
  //VARIABLE DECALRATIONS:  Merijn 2016 Jan 25: here we need to set pointers (super carefully!) for the response matrix we're interested in and the reco spectrum we want to unfold. we copy from jetcorrelation fastjet.cc. WE ASSUME HERE WE JUST PICK THE FILE IN MOSTRECENT, OTHERWISE SPECIFY!
  
  //DECLARATION and assignemnt of parameters for filename of spectrum that is to be unfolded: very important the defaults match the declaration in JetCorrelation_FastJet.cc
  int DataorMonteCarlo, dopAorAp,CastorPosition,SelectedTriggerBit,TriggerBitForEfficiency,BinningKinematicScenario,lowerbound, upperbound,NBins,UnfoldingScenario,NumberFilesToRead,nev;
  double CastorGenJetEnergyCutoff,genjetetalow,genjetetahigh,deltaphicut, CastorRecoJetEnergyCutoff;
  //  InitializeParameters( argc,argv, DataorMonteCarlo, dopAorAp, CastorPosition, SelectedTriggerBit, TriggerBitForEfficiency, CastorRecoJetEnergyCutoff, BinningKinematicScenario, lowerbound,upperbound,NBins, CastorGenJetEnergyCutoff, UnfoldingScenario, genjetetalow, genjetetahigh, deltaphicut, NumberFilesToRead, nev);
  //new variable assignment. initialized in the order argv is initialized
  cout<<"here 1"<<endl;
  InitializeParameters( argc,argv, DataorMonteCarlo, dopAorAp, SelectedTriggerBit, TriggerBitForEfficiency, CastorRecoJetEnergyCutoff, CastorPosition, UnfoldingScenario, BinningKinematicScenario, NumberFilesToRead, nev, lowerbound,upperbound,NBins, CastorGenJetEnergyCutoff, genjetetalow, genjetetahigh, deltaphicut);
  //~VARIABLE ASSIGNMENT spectra
  cout<<"here 2"<<endl;

  //DECLARATION and assignemnt of parameters for input response object.very important the defaults match the declaration in JetCorrelation_FastJet.cc
  int DataorMonteCarloForResponseObject, dopAorApForResponseObject, CastorPositionForResponseObject,SelectedTriggerBitForResponseObject,TriggerBitForEfficiencyForResponseObject,BinningKinematicScenarioForResponseObject,lowerboundForResponseObject, upperboundForResponseObject,NBinsForResponseObject,UnfoldingScenarioForResponseObject,NumberFilesToReadForResponseObject,nevForResponseObject;
  double CastorGenJetEnergyCutoffForResponseObject,genjetetalowForResponseObject,genjetetahighForResponseObject,deltaphicutForResponseObject, CastorRecoJetEnergyCutoffForResponseObject;
  // InitializeParameters(argc,argv, DataorMonteCarloForResponseObject, dopAorApForResponseObject, CastorPositionForResponseObject, SelectedTriggerBitForResponseObject, TriggerBitForEfficiencyForResponseObject, CastorRecoJetEnergyCutoffForResponseObject, BinningKinematicScenarioForResponseObject, lowerboundForResponseObject,upperboundForResponseObject,NBinsForResponseObject, CastorGenJetEnergyCutoffForResponseObject, UnfoldingScenarioForResponseObject, genjetetalowForResponseObject, genjetetahighForResponseObject, deltaphicutForResponseObject, NumberFilesToReadForResponseObject,nevForResponseObject);
  //new variable assignment. initialized in the order argv is initialized
  InitializeParameters(argc,argv, DataorMonteCarloForResponseObject, dopAorApForResponseObject, SelectedTriggerBitForResponseObject, TriggerBitForEfficiencyForResponseObject, CastorRecoJetEnergyCutoffForResponseObject, CastorPositionForResponseObject, UnfoldingScenarioForResponseObject, BinningKinematicScenarioForResponseObject, NumberFilesToReadForResponseObject,nevForResponseObject, lowerboundForResponseObject,upperboundForResponseObject,NBinsForResponseObject, CastorGenJetEnergyCutoffForResponseObject, genjetetalowForResponseObject, genjetetahighForResponseObject, deltaphicutForResponseObject);
  
    //~VARIABLE ASSIGNMENT:
  cout<<"deltaphicut "<<deltaphicut <<endl;
  cout<<"deltaphicutForResponseObject "<<deltaphicutForResponseObject <<endl;

  //BUILD THE INPUT ROOOT FILE PATHS CASTOR JET SPECTRUM:
  TString NTupleDescriptorCastorJetSpectrum=BuildNTupleDescriptor(DataorMonteCarlo, dopAorAp, CastorPosition);
  TString AnalysisDescriptorCastorJetSpectrum=BuildAnalysisDescriptor( SelectedTriggerBit,  CastorRecoJetEnergyCutoff,  lowerbound, NBins,  CastorGenJetEnergyCutoff,  genjetetalow, genjetetahigh,  deltaphicut,DataorMonteCarlo);
  TString DirectoryStringCastorJetSpectrum="../OutputFiles/ResponseObject/"; //here we left the root files with our response objects and particle spectra..
  DirectoryStringCastorJetSpectrum+=NTupleDescriptorCastorJetSpectrum;// dirstring=/user.. add the NTuple description
  DirectoryStringCastorJetSpectrum+=AnalysisDescriptorCastorJetSpectrum;
  cout<<"NTupleDescriptorCastorJetSpectrum  "<<NTupleDescriptorCastorJetSpectrum <<" AnalysisDescriptorCastorJetSpectrum  "<<AnalysisDescriptorCastorJetSpectrum <<endl;
  DirectoryStringCastorJetSpectrum+="/MostRecent/";//needs extra slash..
  vector<TString> FileNameVectorCastorJetSpectrum;
  //  PutDirContentinVector(DirectoryStringCastorJetSpectrum,FileNameVectorCastorJetSpectrum);
  cout<<"dir for spectrum to be unfolded: "<<DirectoryStringCastorJetSpectrum<<endl;
  //~INPUT ROOOT PATHS CASTOR JET SPECTRUM:

  //BUILD THE INPUT ROOOT FILE PATHS RESPONSE OBJECT:
  TString NTupleDescriptorResponseObject=BuildNTupleDescriptor(DataorMonteCarloForResponseObject, dopAorApForResponseObject, CastorPositionForResponseObject);
  TString AnalysisDescriptorResponseObject=BuildAnalysisDescriptor(SelectedTriggerBitForResponseObject,CastorRecoJetEnergyCutoffForResponseObject,lowerboundForResponseObject,NBinsForResponseObject,CastorGenJetEnergyCutoffForResponseObject,genjetetalowForResponseObject,genjetetahighForResponseObject,deltaphicutForResponseObject,DataorMonteCarloForResponseObject);
  TString DirectoryStringResponseObject="../OutputFiles/ResponseObject/"; //here we left the root files with our response objects and particle spectra..
  DirectoryStringResponseObject+=NTupleDescriptorResponseObject;// dirstring=/user.. add the NTuple description
  DirectoryStringResponseObject+=AnalysisDescriptorResponseObject;
  DirectoryStringResponseObject+="/MostRecent/";//needs extra slash..
  vector<TString> FileNameVectorResponseObject;
  //  PutDirContentinVector(DirectoryStringResponseObject,FileNameVectorResponseObject);
  cout<<"dir for response object: "<<DirectoryStringResponseObject<<endl;
  //~INPUT ROOOT PATHS RESPONSE OBJECT:
  
  //BUILD THE OUTPUT ROOOT FILE PATHS:
  TString OutputDirectoryString="../OutputFiles/UnfoldedDistributions/"; //here we left the root files with our response objects and particle spectra..
  //!! We might want to put an additional subdir here specifying what we used to unfold the spectrum
  OutputDirectoryString+=NTupleDescriptorCastorJetSpectrum;// dirstring=/user.. add the NTuple description
  OutputDirectoryString+=AnalysisDescriptorCastorJetSpectrum;
  OutputDirectoryString+="/NBayesIts_";
  OutputDirectoryString+=NBayesIts;
  OutputDirectoryString+="/MostRecent";
  gSystem->mkdir(OutputDirectoryString,kTRUE); //create the directory; kTRUE do recursively if needed
  TString outputdirstring=OutputDirectoryString;
  outputdirstring+="/UnfoldedDistribution.root";
  TFile f2(outputdirstring.Data(),"RECREATE");
  f2.Close();
  cout<<"dir for output files: "<<OutputDirectoryString<<endl;
  //~BUILD THE OUTPUT ROOOT FILE PATHS  

  //Merijn 2016 Feb 20: temporarily set some paths for data unfolding. pA RG nominal and pA data. this set gives reasonable results.
  /* Data and MC:
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/MostRecent/2016_2_16_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";
    FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/DATA/Ap/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/MostRecent/2016_2_18_NrFiles1_NrEvts118524_NrGoodEvts4096.root";
*/
    
  //  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/DATA/pA/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/MostRecent/2016_2_18_NrFiles3_NrEvts119586_NrGoodEvts1689.root";

  //one large MC with itself:
  /*  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/MostRecent/2016_2_29_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";
      FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/MostRecent/2016_2_29_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root"; */

  //small spectrum, large object:
  /*
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_2_29_NrFiles1_NrEvts200000_NrGoodEvts152436.root";
  //  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_2_29_NrFiles1_NrEvts200000_NrGoodEvts152436.root";
  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_2_29_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";
*/

  // different sized samples: Larger response object:
  /*
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles4_NrEvts800000_NrGoodEvts609665.root";

  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
  */
    // different sized samples: smaller response object:

  /*
  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles4_NrEvts800000_NrGoodEvts609665.root";
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
  */
  //same sized different samples:
  /*FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
    FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles4_NrEvts400000_NrGoodEvts304373.root";*/

  // same size:
  /*
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles4_NrEvts800000_NrGoodEvts609665.root";
  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles4_NrEvts800000_NrGoodEvts609665.root"; */
  
  //  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_1_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
  
  //~ 2016 3 4: above is for hadronic jets only

  
  // 2016 3 4: below is for all jets: 
  //same file, all jets: (BRU):
  //need to do trick..
  TString tempstring="";
  FileNameVectorCastorJetSpectrum.push_back(tempstring);
  FileNameVectorResponseObject.push_back(tempstring);
  
  /* ssame size:
  FileNameVectorCastorJetSpectrum[0]="2016_3_3_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";
  FileNameVectorResponseObject[0]="2016_3_3_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";*/
  //different size:

  //small spectrum:
  FileNameVectorCastorJetSpectrum[0]="2016_3_4_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
  FileNameVectorResponseObject[0]="2016_3_3_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root";

  
  /*
  //small spectrum:
  FileNameVectorCastorJetSpectrum[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_4_NrFiles2_NrEvts400000_NrGoodEvts304794.root";
  FileNameVectorResponseObject[0]="/Users/Merijn/PhD/CMSCastorJetAnalysis_Klundert_NB/Analysis/Software/OutputFiles/ResponseObject/ReggeGribov/pA/NominalCastor/TriggerBit35/CastorRecoJetEnergyCutoff150/lowerRECObound150_NBins31/CastorGenJetEnergyCutoff150_genjetetalow-52_genjetetahigh-66_deltaphicut5/2016_3_3_NrFiles24_NrEvts4800000_NrGoodEvts3687890.root"; */

  
  //GET THE RESPONSE OBJECT AND DISTRI WE WANT TO UNFOLD
  cout<<"try to read Castor Jet Spectum from file "<<FileNameVectorCastorJetSpectrum[0]<<endl;
  cout<<"try to read Response Object from file "<<FileNameVectorResponseObject[0]<<endl;
  
  TFile ResponseObjectFile(FileNameVectorResponseObject[0],"READ");
  TFile CastorJetSpectrumFile(FileNameVectorCastorJetSpectrum[0],"READ");
  
  cout<<"before get response object"<<endl;
  
  RooUnfoldResponse* ResponseMatrix=(RooUnfoldResponse*)ResponseObjectFile.Get("response");//response object for 1M evts
  //these are the original histograms
  TH1F* hTrue=(TH1F*)CastorJetSpectrumFile.Get("EnergySpectrumTrueGenJetsResponseObject");
  TH1F* hMeas=(TH1F*)CastorJetSpectrumFile.Get("EnergySpectrumCompleteCalibCastorJetsResponseObject");  

  cout<<"before starting:: hTrue->GetBinContent(1) "<<hTrue->GetBinContent(1)<<endl;
  cout<<"before starting:: hMeas->GetBinContent(1) "<<hMeas->GetBinContent(1)<<endl;

  
  //likely the misses from ResponseObjectFile
  TH1F* mis=(TH1F*)ResponseObjectFile.Get("EnergySpectrumMisses");
  
  //old file that used to work
  /*
  cout<<"note we're working still with an old default Nominal RG in this case.."<<endl;
  RooUnfoldResponse* ResponseMatrix=(RooUnfoldResponse*)oldf.Get("response");//response object for 1M evts
  TH1F* hTrue=(TH1F*)oldf.Get("EnergySpectrumTrueGenJetsResponseObject");
  TH1F* hMeas=(TH1F*)oldf.Get("EnergySpectrumCalibratedCastorJetsResponseObject"); */
  //GET THE RESPONSE OBJECT AND DISTRI WE WANT TO UNFOLD

  //seems not in use currently
  NumberGoodEventsSpectrum=((TH1F*)CastorJetSpectrumFile.Get("CCPlot_ChannelsperEvent"))->GetEntries(); // CCPlot_ChannelsperEvent
  NumberGoodEventsResponseObject=((TH1F*)ResponseObjectFile.Get("CCPlot_ChannelsperEvent"))->GetEntries(); // CCPlot_ChannelsperEvent
  RatioGoodEvtsResponseObjectandSpectrum=float(NumberGoodEventsResponseObject)/float(NumberGoodEventsSpectrum);
  cout<<"NumberGoodEventsSpectrum "<<NumberGoodEventsSpectrum <<" NumberGoodEventsResponseObject "<<NumberGoodEventsResponseObject <<" RatioGoodEvtsResponseObjectandSpectrum "<<RatioGoodEvtsResponseObjectandSpectrum <<endl;
  
  //TH1F* entryhisto=(TH1F*)CastorJetSpectrumFile.Get("CCPlot_ChannelsperEvent");
  // entryhisto->SetTitle("NumberGoodEvents");
  
  //Merijn 2016 Jan 26: we're going to loop over the number of bayesian iterations.
  TH1F* RedChiHisto=new TH1F("RedChiHisto","RedChiHisto",NBayesIts+1,-0.5,NBayesIts+0.5);

  vector<RooUnfoldResponse*> RooUnfoldResponsevector;
  TH1* scaledfakes=ResponseMatrix->Hfakes();
  cout<<"before scaling: scaledfakes->GetBinContent(1) "<<scaledfakes->GetBinContent(1)<<endl;
  scaledfakes->Scale(1./(RatioGoodEvtsResponseObjectandSpectrum));
  //scaledfakes->Scale(RatioGoodEvtsResponseObjectandSpectrum);
  cout<<"after scaling: scaledfakes->GetBinContent(1) "<<scaledfakes->GetBinContent(1)<<endl;
  
  //2016 march 1: compare scaled fakes responseobject to fakes from spectrum file..
  // TH1F* fakecheck=new TH1F();
  //scaledfakes->Copy(*fakecheck);  //(RooUnfoldResponse*)CastorJetSpectrumFile.Get("response")->Hfakes();
  TH1F* fakecheck=scaledfakes->Clone("fakecheck");
  fakecheck->Divide(((RooUnfoldResponse*)CastorJetSpectrumFile.Get("response"))->Hfakes());
  fakecheck->SetName("RatioScaledfakeResponseFakeSpectrum");
  //~2016 march 1: 

  
  for(int BayesianIterationIndex=1;BayesianIterationIndex<NBayesIts;BayesianIterationIndex++){

    RooUnfoldBayes unfold (ResponseMatrix, hMeas, BayesianIterationIndex);    // OR
    //get the unfolded spectrum from the reco spectrum
    TH1* hReco= (TH1*) unfold.Hreco();
    //backsmear the histogram we obtained above   normal way 
    TH1* UnperturbedBacksmeared=ResponseMatrix->ApplyToTruth(hReco); //original
    UnperturbedBacksmeared->Sumw2();
    UnperturbedBacksmeared->Add(scaledfakes);
    //get fakes to add to backsmeared data. Could scale if samples have different sizes.
    TH1* fakes=ResponseMatrix->Hfakes();
    //add the fakes

    //try something with strange histo's and normal response..
    //    RooUnfoldResponse* PerturbedResponsecheck=new RooUnfoldResponse(fakes,fakes,ResponseMatrix->Hresponse()); //check: does anything except the input spectrum matter?
    // RooUnfoldResponse* PerturbedResponsecheck=new RooUnfoldResponse(ResponseMatrix->Hmeasured(),ResponseMatrix->Htruth(),ResponseMatrix->Hresponse()); //they do..
    //   UnperturbedBacksmeared->Sumw2();
    // UnperturbedBacksmeared->Add(scaledfakes);
    
    //~these are the original histograms
    
    vector<TH1*> PerturbedHistoPtrVec;
    TH1F* MeanPerturbedHisto=new TH1F(); //should perhaps properly initialise here? -> YES! or things go crazy. init properly done below
    
    //only input used in this loop: ResponseHistogram, scaledfakes, hMeas,hReco
    TH2* ResponseHistogram=ResponseMatrix->Hresponse(); //   Response matrix as a 2D-histogram: (x,y)=(measured,truth)
    for(int pertIndex=0;pertIndex<NPerturbIterations;pertIndex++){
      
      //perturbed response matrix:
      TH2F* ResponseHistogramS=new TH2F();
      ResponseHistogram->Copy(*ResponseHistogramS); //must copy it otherwise we'll perturb one histo twice
      PerturbResponseHisto(ResponseHistogramS);
      
      //perturbed fakes:
      TH1F* FakeS=new TH1F();
      scaledfakes->Copy(*FakeS);
      Perturb1DHisto(FakeS);

      // Klundert V1 method to get perturbed data:
      //init perturbed response object: smeared matrix, unfolded gen level and measured reco level
      /*
      RooUnfoldResponse* PerturbedResponse=new RooUnfoldResponse(hMeas,hReco,ResponseHistogramS); //note we initialise with the unfolded spectrum
      //  RooUnfoldResponse* PerturbedResponse=new RooUnfoldResponse(FakeS,FakeS,ResponseHistogramS); //check: does anything except the input spectrum matter?
      TH1* DataS=PerturbedResponse->ApplyToTruth(hReco); //backsmear unfolded data with smeared matrix:
      DataS->Add(FakeS);//add unpertubred fares
      */
      
      //Merijn here starts new calculation (anticipate same result):
      //misses:
      TH1F* UnscaledPerturbedFakeS=new TH1F();
      fakes->Copy(*UnscaledPerturbedFakeS);
      Perturb1DHisto(UnscaledPerturbedFakeS);
      
      TH1F* MisS=new TH1F();
      mis->Copy(*MisS);
      Perturb1DHisto(MisS);

      TH1D* PerturbedMatchedData=ResponseHistogramS->ProjectionX(); //anticipate this is correct, i.e. y is gen level, x data level
      Perturb1DHisto(PerturbedMatchedData);
      TH1D* PerturbedMatchedGen=ResponseHistogramS->ProjectionY();
      Perturb1DHisto(PerturbedMatchedGen);
      /*
      cout<<"ResponseHistogramS->GetNbinsX()   "<<ResponseHistogramS->GetNbinsX() <<" ResponseHistogramS->GetNbinsY()   "<<ResponseHistogramS->GetNbinsY() <<endl;
      cout<<"PerturbedMatchedData->GetNbinsX()   "<<PerturbedMatchedData->GetNbinsX() <<" PerturbedMatchedGen->GetNbinsX()   "<<PerturbedMatchedGen->GetNbinsX() <<endl;
      cout<<"UnscaledPerturbedFakeS->GetNbinsX()   "<<UnscaledPerturbedFakeS->GetNbinsX() <<" PerturbedMatchedData->GetNbinsX()   "<<PerturbedMatchedData->GetNbinsX() <<endl;
      cout<<"make TotalPerturbedData"<<endl;
      */
      TH1F* TotalPerturbedData=new TH1F();
      UnscaledPerturbedFakeS->Copy(*TotalPerturbedData);
      TotalPerturbedData->Add(PerturbedMatchedData);

      //      cout<<"make TotalPerturbedTruth"<<endl;
      TH1F* TotalPerturbedTruth=new TH1F();
      MisS->Copy(*TotalPerturbedTruth);
      TotalPerturbedTruth->Add(PerturbedMatchedData);

      //Alex Method: now can make our response object:
      
      RooUnfoldResponse* PerturbedResponse=new RooUnfoldResponse(TotalPerturbedData,TotalPerturbedTruth,ResponseHistogramS); //note we initialise with the unfolded spectrum
      TH1* DataS=PerturbedResponse->ApplyToTruth(hReco); //backsmear unfolded data with smeared matrix:
      //for now add unscaled   DataS->Add(UnscaledPerturbedFakeS);
      DataS->Add(FakeS);//this are the scaled perturbed fakes..
      
      //here we have our backfolded, perturbed distribution. store it and add to a mean matrix
      PerturbedHistoPtrVec.push_back(DataS);
      if(pertIndex==0){DataS->Copy(*MeanPerturbedHisto);
	MeanPerturbedHisto->SetName("MeanPerturbedHisto");
      }  //this makes sure bins etc are properly initialised of histogram for mean 
      else{MeanPerturbedHisto->Add(DataS);}
      
      //display plot histogram and entries in terminal
      // uncomment to plots:
      /*
	TCanvas* aap1=new TCanvas("aap1","testcanvas1",500,500);
	hMeas->SetMarkerStyle(8);
	hMeas->SetMarkerColor(1);
	hMeas->SetDirectory(0);
	hMeas->DrawClone("p");
	DataS->SetMarkerStyle(8);
	DataS->SetMarkerColor(2);
	DataS->DrawClone("p same");
	
	TCanvas* aap2=new TCanvas("aap2","testcanvas2",500,500);
	TH1F* hMeasCopy=new TH1F();
	hMeas->Copy(*hMeasCopy);
	hMeasCopy->Divide(DataS);
	hMeasCopy->DrawClone("p");  
      */
      
      delete PerturbedResponse;
      //    RooUnfoldResponsevector.push_back(PerturbedResponse);
      //      delete PerturbedResponse;
      delete ResponseHistogramS;
      delete FakeS;

      //added histo's to delete:
      delete PerturbedMatchedData;
      delete PerturbedMatchedGen;
      delete TotalPerturbedData;
      delete TotalPerturbedTruth;
      delete MisS;
      delete UnscaledPerturbedFakeS;
    }
    
    
    MeanPerturbedHisto->Scale(1./float(NPerturbIterations)); //scale the mean histogram appropriately:
    cout<<" after filling perturbed and calculating mean histogram"<<endl;
    
    //start build V matrix. This sets input histograms and matrices to equal dimensions. Note that matrix indices start at 0, histogram at 1!
    int NBinsMatrixV=(MeanPerturbedHisto->GetXaxis()->GetNbins()); 
    vector<TMatrixD> VectorMatrices(NPerturbIterations,TMatrixD(NBinsMatrixV,NBinsMatrixV)); // http://stackoverflow.com/questions/6142830/how-do-i-initialize-a-stl-vector-of-objects-who-themselves-have-non-trivial-cons
    
    for(int pertIndex=0;pertIndex<NPerturbIterations;pertIndex++){
      FillMatrixV(VectorMatrices.at(pertIndex),MeanPerturbedHisto,PerturbedHistoPtrVec.at(pertIndex));}//we properly set bin 1 to entry 0 in matrix. We checked the matrix content
    
    //here we have our desired matrices. Now average over them.
    TMatrixD MatrixV(NBinsMatrixV,NBinsMatrixV);
    for(int pertIndex=0;pertIndex<NPerturbIterations;pertIndex++){
      MatrixV+=VectorMatrices.at(pertIndex);}   
    for(int i=0;i<NBinsMatrixV;i++){     
      for(int j=0;j<NBinsMatrixV;j++){
	MatrixV[i][j]/=float(NPerturbIterations);}} //checked our average matrix makes sense
    
    /*cout<<"mean matrix before adding data: "<<endl; for(int i=0;i<NBinsMatrixV;i++){for(int j=0;j<NBinsMatrixV;j++){cout<<MatrixV[i][j]<<"  ";}cout<<endl;} */
    
    //add data to diagonal, checked it works as expected. careful with indices!
    for(int binindex=1;binindex<NBinsMatrixV+1;binindex++){
      MatrixV[binindex-1][binindex-1]+=hMeas->GetBinContent(binindex);}
    
    /*cout<<"mean matrix after adding data: "<<endl; for(int i=0;i<NBinsMatrixV;i++){for(int j=0;j<NBinsMatrixV;j++){cout<<MatrixV[i][j]<<"  ";}cout<<endl;}*/
    
    cout<<"Start Inverting: "<<endl;
    MatrixV.Invert();
    cout<<"After Inverting: matrix "<<endl;
        
    cout<<"Start calc red. Chisquared: "<<endl;
    
    //this calculates reduced chisquared for E>300
    int MatrixBin300=hMeas->FindBin(300)-1;
    int NBinsRedChi=NBinsMatrixV-MatrixBin300;
    double redchiStatingBin300=0;
    for(int i=MatrixBin300;i<NBinsMatrixV;i++){     
      double ctr=0;
      for(int j=MatrixBin300;j<NBinsMatrixV;j++){ctr+=MatrixV[i][j]*(hMeas->GetBinContent(j+1)-UnperturbedBacksmeared->GetBinContent(j+1));}
      redchiStatingBin300+=ctr*(hMeas->GetBinContent(i+1)-UnperturbedBacksmeared->GetBinContent(i+1));}
    
    cout<<"self calculated rec chi E>300  "<<redchiStatingBin300<<endl;  //worked just perfect!
    RedChiHisto->Fill(BayesianIterationIndex,redchiStatingBin300/NBinsRedChi);

  cout<<"before filling:: hTrue->GetBinContent(1) "<<hTrue->GetBinContent(1)<<endl;
  cout<<"before filling:: hMeas->GetBinContent(1) "<<hMeas->GetBinContent(1)<<endl;

    // Fill some histograms of last iteration for fast cross checks..
    if(BayesianIterationIndex==NBayesIts-1){
    TFile f69(outputdirstring.Data(),"UPDATE");
    //unfolded unprturbed particle level
    hReco->SetName("UnperturbedUnfolded");
    hReco->Write();
   
    TH1F* UnperturbedUnfoldedScaledtoGoodEvts=new TH1F();
    hReco->Copy(*UnperturbedUnfoldedScaledtoGoodEvts);
    UnperturbedUnfoldedScaledtoGoodEvts->Sumw2();
    UnperturbedUnfoldedScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
    UnperturbedUnfoldedScaledtoGoodEvts->SetName("UnperturbedUnfoldedScaledtoGoodEvts");
    UnperturbedUnfoldedScaledtoGoodEvts->Write();

    UnperturbedBacksmeared->SetName("UnPerturbedBackSmeared");
    UnperturbedBacksmeared->Write();  //last unperturbed backsmeared distribution. is dependent on nr iterations so must write it here

    TH1F* UnperturbedBackSmearedScaledtoGoodEvts=new TH1F();
    UnperturbedBacksmeared->Copy(*UnperturbedBackSmearedScaledtoGoodEvts);
    UnperturbedBackSmearedScaledtoGoodEvts->Sumw2();
    UnperturbedBackSmearedScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
    UnperturbedBackSmearedScaledtoGoodEvts->SetName("UnperturbedBackSmearedScaledtoGoodEvts");
    UnperturbedBackSmearedScaledtoGoodEvts->Write();    
    
    //smeared distributions:
    PerturbedHistoPtrVec[PerturbedHistoPtrVec.size()-1]->SetName("PerturbedBackSmeared");
    PerturbedHistoPtrVec[PerturbedHistoPtrVec.size()-1]->Write();//write one perturbed backsmeared histogram for last iteration before it gets deleted

    TH1F* PerturbedBackSmearedScaledtoGoodEvts=new TH1F();
    PerturbedHistoPtrVec[PerturbedHistoPtrVec.size()-1]->Copy(*PerturbedBackSmearedScaledtoGoodEvts);
    PerturbedBackSmearedScaledtoGoodEvts->Sumw2();
    PerturbedBackSmearedScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
    PerturbedBackSmearedScaledtoGoodEvts->SetName("PerturbedBackSmearedScaledtoGoodEvts");
    PerturbedBackSmearedScaledtoGoodEvts->Write();
   
    
    TH1F* RatioPerturbedBackSmearedvsData=new TH1F();
    hMeas->Copy(*RatioPerturbedBackSmearedvsData);
    RatioPerturbedBackSmearedvsData->Sumw2();
    RatioPerturbedBackSmearedvsData->Divide(PerturbedHistoPtrVec[PerturbedHistoPtrVec.size()-1]);
    RatioPerturbedBackSmearedvsData->SetName("RatioPerturbedBackSmearedvsData");
    RatioPerturbedBackSmearedvsData->Write();//ratio plot last perturbed backsmeared w.r.t. real data
    
    TH1F* RatioUnPerturbedBackSmearedvsData=new TH1F();
    hMeas->Copy(*RatioUnPerturbedBackSmearedvsData);
    RatioUnPerturbedBackSmearedvsData->Sumw2();
    RatioUnPerturbedBackSmearedvsData->Divide(UnperturbedBacksmeared);
    RatioUnPerturbedBackSmearedvsData->SetName("RatioUnPerturbedBackSmearedvsData");
    RatioUnPerturbedBackSmearedvsData->Write();//ratio plot last perturbed backsmeared w.r.t. real data
    
    MeanPerturbedHisto->SetName("MeanPerturbedBacksmeared");
    MeanPerturbedHisto->Write(); //average backsmeared histogram
    
    TH1F* MeanPerturbedBacksmearedScaledtoGoodEvts=new TH1F();
    MeanPerturbedHisto->Copy(*MeanPerturbedBacksmearedScaledtoGoodEvts);
    MeanPerturbedBacksmearedScaledtoGoodEvts->Sumw2();
    MeanPerturbedBacksmearedScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
    MeanPerturbedBacksmearedScaledtoGoodEvts->SetName("MeanPerturbedBacksmearedScaledtoGoodEvts");
    MeanPerturbedBacksmearedScaledtoGoodEvts->Write();

    TH1F* RatioMeanPerturbedBacksmearedvsData=new TH1F();
    hMeas->Copy(*RatioMeanPerturbedBacksmearedvsData);
    RatioMeanPerturbedBacksmearedvsData->Sumw2();
    RatioMeanPerturbedBacksmearedvsData->Divide(MeanPerturbedHisto);
    RatioMeanPerturbedBacksmearedvsData->SetName("RatioMeanPerturbedBackSmearedvsData");
    RatioMeanPerturbedBacksmearedvsData->Write();//ratio plot last perturbed backsmeared w.r.t. real dat

    TH1F* RatioUnfoldedvsTruth=new TH1F();
    hReco->Copy(*RatioUnfoldedvsTruth);
    RatioUnfoldedvsTruth->Sumw2();
    RatioUnfoldedvsTruth->Divide(hTrue);
    RatioUnfoldedvsTruth->SetName("RatioUnfoldedvsTruth");
    RatioUnfoldedvsTruth->Write();//ratio plot last perturbed backsmeared w.r.t. real dat

    fakecheck->Write();
    f69.Close();
    //~~ Fill some histograms for fast cross checks..

    }
    
    //cleanup:
    cout<<"start deletion"<<endl;//this seems to work..
    for (int delindex=0;delindex<PerturbedHistoPtrVec.size();delindex++){
      delete PerturbedHistoPtrVec[delindex];}
    VectorMatrices.clear();
    
    delete MeanPerturbedHisto;  
    for(int i=0;i<RooUnfoldResponsevector.size();i++){
      delete RooUnfoldResponsevector[i];}
    delete hReco;
    delete UnperturbedBacksmeared;

    //
    // delete fakes; //not necessarily in loop .. seems it causes a crahs but could also be hReco?
    
    cout<<"end deletion"<<endl;    
  }//closes loop Bayesion its.

    

  
  TFile f69(outputdirstring.Data(),"UPDATE");
  //here all unperturbed things
  hMeas->Write();//measure distri
  
  TH1F* hMeasScaledtoGoodEvts=new TH1F();
  //hMeasScaledtoGoodEvts->Sumw2();
  hMeas->Copy(*hMeasScaledtoGoodEvts);
  hMeasScaledtoGoodEvts->Sumw2();
  hMeasScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
  hMeasScaledtoGoodEvts->SetName("EnergySpectrumCalibratedCastorJetsResponseObjectScaledtoGoodEvts");
  hMeasScaledtoGoodEvts->Write();
  RedChiHisto->Write();

  //write copy of particle level distri and scaled copy to nr events as reference:
  hTrue->Write();
  TH1F* hTrueScaledtoGoodEvts=new TH1F();
  hTrue->Copy(*hTrueScaledtoGoodEvts);
  hTrueScaledtoGoodEvts->Sumw2();
  hTrueScaledtoGoodEvts->Scale(1./float(NumberGoodEventsSpectrum));
  hTrueScaledtoGoodEvts->SetName("hTrueScaledtoGoodEvts");
  hTrueScaledtoGoodEvts->Write();
  
  f69.Close();
  cout<<" output writte to "<<outputdirstring<<endl;
  ResponseObjectFile.Close();
  CastorJetSpectrumFile.Close();
  //oldf.Close();
  
  delete ResponseMatrix;

  cout<<"program reached end type eneter to terminate"<<endl;
  //system("pause") might work as well http://stackoverflow.com/questions/21257544/c-wait-for-user-input
  char c;
  std::cin >> c;
  
  return 0;
} //closes macro

  
TString* GenerateTreeNames(TFile* f, int& nkeys){
  nkeys=f->GetListOfKeys()->GetSize();
  // int* CycleNumber[nkeys];
TString* TTreenames=new TString[nkeys];
  for(int i=0;i<nkeys;i++){ 
    TTreenames[i]=((TKey*)f->GetListOfKeys()->At(i))->GetName();
    TTreenames[i]+=";";
    TTreenames[i]+=((TKey*)f->GetListOfKeys()->At(i))->GetCycle();} 
    cout<<"nkeys"<<nkeys<<"TTreenames[0]"<<TTreenames[0]<<endl;
    return TTreenames;}

	    //root list content of directory cern  https://root.cern.ch/phpBB3/viewtopic.php?t=12793
void PutDirContentinVector(TString Directory, vector<TString>& FileNameVector){
  cout<<"PutDirContentinVector "<<Directory <<endl;
  const char *ext=".root";
  TSystemDirectory dir(Directory, Directory);//why 2 times same name?
  TList *files = dir.GetListOfFiles();
  if (files) {
    TSystemFile *file;
    TString fname;
    TIter next(files);
    while ((file=(TSystemFile*)next())) {
      fname = file->GetName();
      //cout<<"files contains "<< fname<<endl;

      if (!file->IsDirectory() && fname.EndsWith(ext)) {
	TString DUmmyFileName=Directory;
	DUmmyFileName+=fname.Data();
	FileNameVector.push_back(DUmmyFileName);
	cout << fname.Data() << endl;
      }
    }
  }

  //try to clean up
    if (files) {
    TSystemFile *file;
    TString fname;
    TIter next(files);
    while ((file=(TSystemFile*)next())) {
      fname = file->GetName();
      //cout<<"about to delete "<< fname<<endl;
    }
  }
    files->Clear();
    delete files;
}

      
void PerturbResponseHisto(TH2F* ResponseHisto){
  
  for(int nbx=1;nbx<ResponseHisto->GetXaxis()->GetNbins()+1;nbx++){
    for(int nby=1;nby<ResponseHisto->GetYaxis()->GetNbins()+1;nby++){  
      // for(int nbx=10;nbx<12;nbx++){
      //for(int nby=10;nby<12;nby++){
      // for(int nb=0;nby<ResponseHisto->GetYaxis()->GetNbins();nby++){
      
      double response=ResponseHisto->GetBinContent(nbx,nby);
      TRandom* random1 = new TRandom;
      random1->SetSeed(0); //results will differ per iteration
      double newresponse=random1->Poisson(response);
      // cout<<"response "<<response <<" perturnbed response "<<newresponse<<endl;// seems to work
      //fill the histogram back:
      int bin=ResponseHisto->GetBin(nbx,nby);
      ResponseHisto->SetBinContent(bin,newresponse);
      // cout<<"matrix: response "<<response<<" newresponse "<<newresponse<<" response-newresponse ^2 "<<TMath::Power(response-newresponse,2)<<endl;
    
      //double check=ResponseHisto->GetBinContent(nbx,nby);
      //cout<<"check "<<check <<endl;
      //~ proceude checked..
      delete random1;
    }	    
  }	
}

void Perturb1DHisto(TH1* ResponseHisto){
  
  //  for(int nbx=0;nbx<ResponseHisto->GetXaxis()->GetNbins();nbx++){
  for(int nbx=1;nbx<ResponseHisto->GetXaxis()->GetNbins()+1;nbx++){
    // for(int nbx=10;nbx<12;nbx++){
    double response=ResponseHisto->GetBinContent(nbx);
    //cout<<"response "<<response <<endl;
    
    TRandom* random1 = new TRandom;
    random1->SetSeed(0);
    double newresponse=random1->Poisson(response);
    //cout<<"response "<<response<<" newresponse "<<newresponse<<" response-newresponse ^2 "<<TMath::Power(response-newresponse,2)<<endl;

    //fill the histogram back:
    
    int bin=ResponseHisto->GetBin(nbx);
    ResponseHisto->SetBinContent(bin,newresponse);
    //double check=ResponseHisto->GetBinContent(nbx);
    //cout<<"check "<<check <<endl;
    //~ procedure checked..
    delete random1;
  }	
}

void FillMatrixV(TMatrixD & V, TH1* mean , TH1* Perturbed){
  //  cout<<"start FillMatrixV"<<endl;
  // for(int nbx=0;nbx<Perturbed->GetXaxis()->GetNbins();nbx++){//Merijn 2016 Jan 19 go from [1,NBins] instead.
    //cout<<"nbx "<<nbx<<endl;
    //  for(int nby=0;nby<Perturbed->GetXaxis()->GetNbins();nby++){
  
  for(int nbx=1;nbx<Perturbed->GetXaxis()->GetNbins()+1;nbx++){
    for(int nby=1;nby<Perturbed->GetXaxis()->GetNbins()+1;nby++){
      
      //    cout<<"nby "<<nby<<endl;
      //  V[nbx][nby]=(Perturbed->GetBinContent(nbx) - mean->GetBinContent(nbx)) * (Perturbed->GetBinContent(nby) - mean->GetBinContent(nby));
      //we loop over histograms starting at 1 but fill matrix at 0. Careful with this..!
      V[nbx-1][nby-1]=(Perturbed->GetBinContent(nbx) - mean->GetBinContent(nbx)) * (Perturbed->GetBinContent(nby) - mean->GetBinContent(nby));


      //cout<<"V[nbx][nby]" <<V[nbx][nby];
      //      if(nby==10) {cout<<"V[nbx][nby]" <<V[nbx][nby];}
      // if(nby==10) { cout<<nbx<<"  Perturbed->GetBinContent(nbx) "<<Perturbed->GetBinContent(nbx)<<"mean->GetBinContent(nbx) "<<mean->GetBinContent(nbx)<<endl;cout<<"V[nbx][nby] "<<V[nbx][nby]<<endl;}
      //cout<<nbx<<"  Perturbed->GetBinContent(nby) "<<Perturbed->GetBinContent(nby)<<"mean->GetBinContent(nby) "<<mean->GetBinContent(nby)<<endl;
    }	    
  }
  //    cout<<"end FillMatrixV"<<endl;
  return ;
}
 
 TString BuildNTupleDescriptor(int DataorMonteCarlo, int dopAorAp, int CastorPosition){
   TString NTupleDescriptor="";
   if(DataorMonteCarlo==0)NTupleDescriptor="DATA/";
   if(DataorMonteCarlo==1)NTupleDescriptor="ReggeGribov/";
   if(DataorMonteCarlo==2)NTupleDescriptor="Hijing/";
   if(dopAorAp==0)NTupleDescriptor+="pA/";
   if(dopAorAp==1)NTupleDescriptor+="Ap/";
   if(DataorMonteCarlo>0&&CastorPosition==0)NTupleDescriptor+="NominalCastor/";
   if(DataorMonteCarlo>0&&CastorPosition==1)NTupleDescriptor+="DisplacedCastor/";  
   return NTupleDescriptor;
 }
 
//TString BuildAnalysisDescriptor(int SelectedTriggerBit, int CastorRecoJetEnergyCutoff, int lowerbound,int NBins, int CastorGenJetEnergyCutoff, int genjetetalow,int genjetetahigh, int deltaphicut, int DataorMonteCarlo){
TString BuildAnalysisDescriptor(int SelectedTriggerBit, double CastorRecoJetEnergyCutoff, int lowerbound,int NBins, double CastorGenJetEnergyCutoff, double genjetetalow, double genjetetahigh, double deltaphicut, int DataorMonteCarlo){
   TString AnalysisDescriptor="";   
   AnalysisDescriptor+="TriggerBit";  //data and MC:
   AnalysisDescriptor+=SelectedTriggerBit;
   AnalysisDescriptor+="/CastorRecoJetEnergyCutoff";
   AnalysisDescriptor+=CastorRecoJetEnergyCutoff;
   AnalysisDescriptor+="/lowerRECObound";
   AnalysisDescriptor+=lowerbound;
   AnalysisDescriptor+="_NBins";
   AnalysisDescriptor+=NBins;
   if(DataorMonteCarlo>0){  //MC specific
     AnalysisDescriptor+="/CastorGenJetEnergyCutoff";
     AnalysisDescriptor+=CastorGenJetEnergyCutoff;
     AnalysisDescriptor+="_genjetetalow";
     AnalysisDescriptor+=int(10*genjetetalow);
     AnalysisDescriptor+="_genjetetahigh";
     AnalysisDescriptor+=int(10*genjetetahigh);
     AnalysisDescriptor+="_deltaphicut";
     AnalysisDescriptor+=int(10*deltaphicut);
     cout<<"deltaphicut "<<deltaphicut<<endl;
     cout<<"int(10*deltaphicut) "<<int(10*deltaphicut)<<endl;
       }
   
   return AnalysisDescriptor;
 }
  

//void InitializeParameters(int argc,char **argv, int & DataorMonteCarlo, int & dopAorAp, int & CastorPosition, int & SelectedTriggerBit,int & TriggerBitForEfficiency, double & CastorRecoJetEnergyCutoff,int & BinningKinematicScenario, int & lowerbound,int & upperbound,int & NBins, double & CastorGenJetEnergyCutoff, int & UnfoldingScenario, double & genjetetalow, double & genjetetahigh, double & deltaphicut, int & NumberFilesToRead, int & nev){
void InitializeParameters(int argc,char **argv, int & DataorMonteCarlo, int & dopAorAp, int & SelectedTriggerBit,int & TriggerBitForEfficiency, double & CastorRecoJetEnergyCutoff, int & CastorPosition, int & UnfoldingScenario, int & BinningKinematicScenario, int & NumberFilesToRead, int & nev, int & lowerbound,int & upperbound,int & NBins, double & CastorGenJetEnergyCutoff, double & genjetetalow, double & genjetetahigh, double & deltaphicut){


  if(argc!=1){//non-default init
    DataorMonteCarlo=atoi(argv[1]); //per default we do MC (>0) and RG (==1)  for (int i = 0; i < argc; ++i)   std::cout << argv[i]<<"    "<<atoi(argv[i]) <<std::endl;
    
    dopAorAp=atoi(argv[2]);//def=pA, proton to Castor. 1 is for Ap.
    SelectedTriggerBit=atoi(argv[3]);
    TriggerBitForEfficiency=atoi(argv[4]);
    CastorRecoJetEnergyCutoff=atoi(argv[5]);
    //defs for set 2:
    CastorGenJetEnergyCutoff=CastorRecoJetEnergyCutoff; 
    CastorPosition=atoi(argv[6]); //def for MC is displaced
    UnfoldingScenario=atoi(argv[7]); //this are the best values according to Alex..
    BinningKinematicScenario=atoi(argv[8]);
    //defs for set 3:
    NumberFilesToRead=atoi(argv[9]);// -1 means take all files 
    nev=atoi(argv[10]);} //-1 takes all events     
  
  else{ //with these 10 parameters we can steer the complete macro if we don need to specify unfolding scenario and BinningKinematicScenario
    //defs for set 1:
    DataorMonteCarlo=1; //per default we do MC (>0) and RG (==1). 
    dopAorAp=0;//def=pA, proton to Castor. 1 is for Ap.
    SelectedTriggerBit=35;
    TriggerBitForEfficiency=0;
    CastorRecoJetEnergyCutoff=150;
    //defs for set 2:
    CastorGenJetEnergyCutoff=CastorRecoJetEnergyCutoff; 
    CastorPosition=1; //def for MC is displaced
    UnfoldingScenario=0;//this are the best values according to Alex..
    BinningKinematicScenario=0;
    //defs for set 3:
    NumberFilesToRead=1;// -1 means take all files 
    nev=10000; //-1 takes all events 
  }
  
  //these get specified if default is set.
  if(UnfoldingScenario==0){
    genjetetalow=-5.2; //only used in CutEnergyOrderedGenJetVectoronEta
    genjetetahigh=-6.6;
    deltaphicut=0.5;} //take hardest genjet in deltaphicut around Castor Jet
  else{
    genjetetalow=-5.2; //only used in CutEnergyOrderedGenJetVectoronEta
    genjetetahigh=-6.6;
    deltaphicut=0.5; }
  if(BinningKinematicScenario==0&&CastorRecoJetEnergyCutoff==150){  //binning, ranges could be extended with 2 bins and 150 GeV perhaps for larger statistics. 
    lowerbound=CastorRecoJetEnergyCutoff;
    upperbound=2475;
    NBins=31;}
  if(BinningKinematicScenario==0&&CastorRecoJetEnergyCutoff==0){
    lowerbound=CastorRecoJetEnergyCutoff;
    upperbound=2475;
    NBins=33;}
   if(BinningKinematicScenario==1&&CastorRecoJetEnergyCutoff==150){  //binning, ranges could be extended with 2 bins and 150 GeV perhaps for larger statistics. 
    lowerbound=CastorRecoJetEnergyCutoff;
    upperbound=2100;
    NBins=26;}
  if(BinningKinematicScenario==1&&CastorRecoJetEnergyCutoff==0){
    lowerbound=CastorRecoJetEnergyCutoff;
    upperbound=2100;
    NBins=28;}
  return;

}
 


//thgis are old matrix cross check calculations
    //here scale the matrix first!
    /*
    //~    
    TMatrixD DeltaColVect(NBinsMatrixV,1); //column vector
    TMatrixD DeltaRowVect(1,NBinsMatrixV); //row vector
    
    for(int binindex=1;binindex<NBinsMatrixV+1;binindex++){
      DeltaColVect[binindex-1][0]=(hMeas->GetBinContent(binindex)-UnperturbedBacksmeared->GetBinContent(binindex));
      DeltaRowVect[0][binindex-1]=(hMeas->GetBinContent(binindex)-UnperturbedBacksmeared->GetBinContent(binindex));}
    
    // TMatrixD Step1(NBinsMatrixV,1);
    TMatrixD Step1(NBinsMatrixV,1);
    Step1.Mult(MatrixV,DeltaColVect);
    
    //for(int binindex=1;binindex<NBinsMatrixV+1;binindex++){cout<<"Step1[binindex-1][0]"<<Step1[binindex-1][0]<<endl; } //check positive.
    TMatrixD Step2(1,1);
    Step2.Mult(DeltaRowVect,Step1);
    cout<<"Proper Reduced Chisquared: "<<Step2[0][0] <<endl;  //check positive.

    //Here we have a proper reduced Chisquared taking bin content into account.
    
    //here try to reproduce previous result with self built matrix multiplication:
    double totalredchi=0;
    for(int i=0;i<NBinsMatrixV;i++){     
      double ctr=0;
      for(int j=0;j<NBinsMatrixV;j++){ctr+=MatrixV[i][j]*(hMeas->GetBinContent(j+1)-UnperturbedBacksmeared->GetBinContent(j+1));}
      totalredchi+=ctr*(hMeas->GetBinContent(i+1)-UnperturbedBacksmeared->GetBinContent(i+1));}
    //    cout<<"self calculated rec chi "<<totalredchi<<endl;  //worked just perfect! note binning of histograms lags 1 behind w.r.t. matrix
    */