SchieferD/jetcalib/mcresponse_x.cpp

////////////////////////////////////////////////////////////////////////////////
//
// mcresponse_x
// ------------
//
//            06/16/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
////////////////////////////////////////////////////////////////////////////////


#include "utils/cmdline.h"
#include "utils/hist.h"

#include <TROOT.h>
#include <TRint.h>
#include <TFile.h>
#include <TTree.h>
#include <TEventList.h>
#include <TH1F.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TLegend.h>

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <string>
#include <vector>


using namespace std;


////////////////////////////////////////////////////////////////////////////////
// defaults
////////////////////////////////////////////////////////////////////////////////
string dflt_etbins ="7,9,12,15,20,25,30,45,65,90,120";
string dflt_etabins=".2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.3,2.5,2.7,3.0,3.3";
string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";


////////////////////////////////////////////////////////////////////////////////
// declare global functions
////////////////////////////////////////////////////////////////////////////////
void     replaceHistos(int nbins,TH1F**& hVar,TH1F**& hEtCorr,TTree**& tVar);
void     fitHistos    (TH1F**  histos,int n,double nsigma);
void     fitHistos    (TH1F*** histos,int n1,int n2,double nsigma);
void     saveHistos   (TH1F**  histos,int n);
void     saveHistos   (TH1F*** histos,int n1,int n2);


////////////////////////////////////////////////////////////////////////////////
// main
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
int main(int argc,char**argv)
{
  cmdline cl;
  cl.parse(argc,argv);
  
  vector<string> input     = cl.get_vector<string>("input");
  string         output    = cl.get_value<string> ("output","mcresponse.root");
  string         selection = cl.get_value<string> ("selection",       "njt>0");
  float          drmax     = cl.get_value<float>  ("drmax",               0.3);
  short          applyjes  = cl.get_value<short>  ("applyjes",              0);
  string         treename  = cl.get_value<string> ("treename",            "t");
  float          nsigma    = cl.get_value<float>  ("nsigma",              3.0);
  vector<float>  etbins    = cl.get_vector<float> ("etbins",      dflt_etbins);
  vector<float>  etabins   = cl.get_vector<float> ("etabins",    dflt_etabins);
  vector<float>  phibins   = cl.get_vector<float> ("phibins",    dflt_phibins);
  bool           abseta    = cl.get_value<bool>   ("abseta",             true);
  
  if (!cl.check()) return 0;
  cl.print();
  

  // etabins
  if (!abseta) {
    int neta=(int)etabins.size();
    std::reverse(etabins.begin(),etabins.end());
    for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
    for (int ieta=0;ieta<neta;ieta++)    etabins[ieta]*=-1;
  }
  
  
  // jetalgs
  vector<string> jetalgs;
  for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
    string jetalg = *it;
    string::size_type pos;
    while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
    jetalg = jetalg.substr(0,jetalg.find(".root"));
    jetalgs.push_back(jetalg);
  }
  
  // prepare branch variables and histograms / graphs
  float weight;
  char  njt;
  float jtet[100];
  float jtpt[100];
  float jteta[100];
  float jtphi[100];
  float jtgenet[100];
  float jtgenpt[100];
  float jtgendr[100];
  float jtjes[100][3];
  
  int netbins  = etbins.size()+1;
  int netabins = etabins.size()+1;
  int nphibins = phibins.size()+1;
  
  vector<TH1F**>         hGenEt;
  vector<TH1F**>         hEtCorrGenEt;
  vector<TH1F**>         hAbsRspGenEt;
  vector<TH1F**>         hRelRspGenEt;

  vector<TH1F**>         hEt;
  vector<TH1F**>         hEtCorrEt;
  vector<TH1F**>         hAbsRspEt;
  vector<TH1F**>         hRelRspEt;

  vector<TH1F**>         hEta;
  vector<TH1F**>         hEtCorrEta;
  vector<TH1F**>         hAbsRspEta;
  vector<TH1F**>         hRelRspEta;

  vector<TH1F**>         hPhi;
  vector<TH1F**>         hEtCorrPhi;
  vector<TH1F**>         hAbsRspPhi;
  vector<TH1F**>         hRelRspPhi;
  
  vector<TH1F***>        hEtEtaEt;
  vector<TH1F***>        hEtCorrEtaEt;
  vector<TH1F***>        hAbsRspEtaEt;
  vector<TH1F***>        hRelRspEtaEt;
  
  vector<TGraphErrors*>  gRspVsGenEt;
  vector<TGraphErrors*>  gSgmVsGenEt;
  vector<TGraphErrors*>  gRspVsEt;
  vector<TGraphErrors*>  gSgmVsEt;
  vector<TGraphErrors*>  gRspVsEta;
  vector<TGraphErrors*>  gSgmVsEta;
  vector<TGraphErrors*>  gRspVsPhi;
  vector<TGraphErrors*>  gSgmVsPhi;
  
  vector<TGraphErrors**> gRspVsEtaEt;
  vector<TGraphErrors**> gSgmVsEtaEt;
  
  argc=3; argv[1]="-l"; argv[2]="-b";
  TRint* app = new TRint(argv[0],&argc,argv);
  
  
  //
  // loop over all input files
  //
  TFile* plotfile = new TFile(output.c_str(),"RECREATE");
  
  for (unsigned int i=0;i<input.size();++i) {
    TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
    TTree* t = (TTree*)f->Get(treename.c_str());
    
    TEventList* el = new TEventList("sel","sel");
    t->Draw(">>sel",selection.c_str());

    t->SetBranchAddress("weight", &weight);
    t->SetBranchAddress("njt",    &njt);
    t->SetBranchAddress("jtet",   jtet);
    t->SetBranchAddress("jtpt",   jtpt);
    t->SetBranchAddress("jteta",  jteta);
    t->SetBranchAddress("jtphi",  jtphi);
    t->SetBranchAddress("jtgenet",jtgenet);
    t->SetBranchAddress("jtgenpt",jtgenpt);
    t->SetBranchAddress("jtgendr",jtgendr);

    if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
    else
      for (int i1=0;i1<100;i1++)
        for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
    
    
    string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
    string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
    
    hGenEt      .push_back(hist::initHistos("GenEt",jetalgs[i],
                                            100,"jet E_{T}^{gen} [GeV]",
                                            "GenEt",etbins));
    hEtCorrGenEt.push_back(hist::initHistos("EtCorrGenEt",jetalgs[i],
                                            100,"jet E_{T}^{gen} [GeV]",
                                            "GenEt",etbins));
    hAbsRspGenEt.push_back(hist::initHistos("AbsRspGenEt",jetalgs[i],
                                            100,-50,150,absrsp_xtitle,
                                            "GenEt",etbins));
    hRelRspGenEt.push_back(hist::initHistos("RelRspGenEt",jetalgs[i],
                                            100,0,2,relrsp_xtitle,
                                            "GenEt",etbins));
    
    hEt         .push_back(hist::initHistos("Et",jetalgs[i],
                                            100,"jet E_{T} [GeV]",
                                            "Et",etbins));
    hEtCorrEt   .push_back(hist::initHistos("EtCorrEt",jetalgs[i],
                                            100,"jet E_{T} [GeV]",
                                            "Et",etbins));
    hAbsRspEt   .push_back(hist::initHistos("AbsRspEt",jetalgs[i],
                                            100,-50,150,absrsp_xtitle,
                                            "Et",etbins));
    hRelRspEt   .push_back(hist::initHistos("RelRspEt",jetalgs[i],
                                            100,0,2,relrsp_xtitle,
                                            "Et",etbins));
    
    hEta        .push_back(hist::initHistos("Eta",jetalgs[i],
                                            100,"jet #eta","Eta",etabins));
    hEtCorrEta  .push_back(hist::initHistos("EtCorrEta",jetalgs[i],
                                            100,"jet E_{T} [GeV]",
                                            "Eta",etabins));
    hAbsRspEta  .push_back(hist::initHistos("AbsRspEta",jetalgs[i],
                                            100,-50,150,absrsp_xtitle,
                                            "Eta",etabins));
    hRelRspEta  .push_back(hist::initHistos("RelRspEta",jetalgs[i],
                                            100,0,2,relrsp_xtitle,
                                            "Eta",etabins));
    
    hPhi        .push_back(hist::initHistos("Phi",jetalgs[i],
                                            100,"jet #phi","Phi",phibins));
    hEtCorrPhi  .push_back(hist::initHistos("EtCorrPhi",jetalgs[i],
                                            100,"jet E_{T} [Gev]",
                                            "Phi",phibins));
    hAbsRspPhi  .push_back(hist::initHistos("AbsRspPhi",jetalgs[i],
                                            100,-50,150,absrsp_xtitle,
                                            "Phi",phibins));
    hRelRspPhi  .push_back(hist::initHistos("RelRspPhi",jetalgs[i],
                                            100,0,2,relrsp_xtitle,
                                            "Phi",phibins));
    
    hEtEtaEt    .push_back(hist::initHistos("Et",jetalgs[i],
                                            100,"jet E_{T} [GeV]",
                                            "Eta",etabins,"Et",etbins));
    hEtCorrEtaEt.push_back(hist::initHistos("EtCorrEtaEt",jetalgs[i],
                                            100,"jet E_{T} [GeV]",
                                            "Eta",etabins,"Et",etbins));
    hAbsRspEtaEt.push_back(hist::initHistos("AbsRspEtaEt",jetalgs[i],
                                            100,-50,150,absrsp_xtitle,
                                            "Eta",etabins,"Et",etbins));
    hRelRspEtaEt.push_back(hist::initHistos("RelRspEtaEt",jetalgs[i],
                                            100,0,2,relrsp_xtitle,
                                            "Eta",etabins,"Et",etbins));
    
    
    int nevts = el->GetN();
    cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;

    // prepare intermediate trees
    float genet,et,eta,phi,etcorr;
    TTree**  tGenEt   = new TTree*[netbins];
    TTree**  tEt      = new TTree*[netbins];
    TTree**  tEta     = new TTree*[netabins];
    TTree**  tPhi     = new TTree*[nphibins];
    TTree*** tEtEtaEt = new TTree**[netabins];

    for (int iet=0;iet<netbins;iet++) {
      tGenEt[iet]=new TTree("tGenEt","tGenEt");
      tGenEt[iet]->Branch("val",   &genet, "val/F");
      tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
      tGenEt[iet]->Branch("weight",&weight,"weight/F");
      tEt[iet]=new TTree("tEt","tEt");
      tEt[iet]->Branch("val",   &et,    "val/F");
      tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
      tEt[iet]->Branch("weight",&weight,"weight/F");
    }
    for (int ieta=0;ieta<netabins;ieta++) {
      tEta[ieta]=new TTree("tEta","tEta");
      tEta[ieta]->Branch("val",   &eta,   "val/F");
      tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
      tEta[ieta]->Branch("weight",&weight,"weight/F");
    }
    for (int iphi=0;iphi<nphibins;iphi++) {
      tPhi[iphi]=new TTree("tPhi","tPhi");
      tPhi[iphi]->Branch("val",   &phi,   "val/F");
      tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
      tPhi[iphi]->Branch("weight",&weight,"weight/F");
    }
    for (int ieta=0;ieta<netabins;ieta++) {
      tEtEtaEt[ieta] = new TTree*[netbins];
      for (int iet=0;iet<netbins;iet++) {
        tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
        tEtEtaEt[ieta][iet]->Branch("val",    &et,   "val/F");
        tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
        tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
      }
    }
    
    
    // loop over all events
    for (int ievt=0;ievt<nevts;ievt++) {
      
      t->GetEntry(el->GetEntry(ievt));
      
      // loop over all jets in the event
      for (int ijt=0;ijt<njt;ijt++) {
        
        if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;

        if (ijt>0) continue;

        et      = jtet[ijt];
        genet   = jtgenet[ijt];
        eta     = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
        phi     = jtphi[ijt];
        etcorr  = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
        
        float absrsp  = genet - etcorr;
        float relrsp  = etcorr/genet;

        // genet
        int igenet = hist::get_ibin(genet,etbins);
        hAbsRspGenEt[i][igenet] ->Fill(absrsp,weight);
        hRelRspGenEt[i][igenet] ->Fill(relrsp,weight);
        tGenEt[igenet]->Fill();
        
        // et
        int iet = hist::get_ibin(et,etbins);
        hAbsRspEt[i][iet]->Fill(absrsp,weight);
        hRelRspEt[i][iet]->Fill(relrsp,weight);
        tEt[iet]->Fill();
        
        // eta
        int ieta = hist::get_ibin(eta,etabins);
        hAbsRspEta[i][ieta]->Fill(absrsp,weight);
        hRelRspEta[i][ieta]->Fill(relrsp,weight);
        tEta[ieta]->Fill();
        
        // phi
        int iphi = hist::get_ibin(phi,phibins);
        hAbsRspPhi[i][iphi]->Fill(absrsp,weight);
        hRelRspPhi[i][iphi]->Fill(relrsp,weight);
        tPhi[iphi]->Fill();
        
        // et-eta
        hAbsRspEtaEt[i][ieta][iet]->Fill(absrsp,weight);
        hRelRspEtaEt[i][ieta][iet]->Fill(relrsp,weight);
        tEtEtaEt[ieta][iet]->Fill();
        
      } // jets
    } // evts
    
    cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
    

    // replace histograms
    replaceHistos(netbins,hGenEt[i],hEtCorrGenEt[i],tGenEt);
    replaceHistos(netbins,hEt[i],   hEtCorrEt[i],      tEt);
    replaceHistos(netabins,hEta[i], hEtCorrEta[i],    tEta);
    replaceHistos(nphibins,hPhi[i], hEtCorrPhi[i],    tPhi);
    for (int ieta=0;ieta<netabins;ieta++)
      replaceHistos(netbins,hEtEtaEt[i][ieta],hEtCorrEtaEt[i][ieta],tEtEtaEt[ieta]);
    
    TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
    d->cd();
    
    // fit response histos
    fitHistos(hAbsRspGenEt[i],netbins,nsigma);
    fitHistos(hRelRspGenEt[i],netbins,nsigma);

    fitHistos(hAbsRspEt[i],   netbins,nsigma);
    fitHistos(hRelRspEt[i],   netbins,nsigma);

    fitHistos(hAbsRspEta[i],  netabins,nsigma);
    fitHistos(hRelRspEta[i],  netabins,nsigma);

    fitHistos(hAbsRspPhi[i],  nphibins,nsigma);
    fitHistos(hRelRspPhi[i],  nphibins,nsigma);
    
    fitHistos(hAbsRspEtaEt[i],netabins,netbins,nsigma);
    fitHistos(hRelRspEtaEt[i],netabins,netbins,nsigma);

    
    // save histograms
    saveHistos(hGenEt[i],      netbins);
    saveHistos(hEtCorrGenEt[i],netbins);
    saveHistos(hAbsRspGenEt[i],netbins);
    saveHistos(hRelRspGenEt[i],netbins);

    saveHistos(hEt[i],         netbins);
    saveHistos(hEtCorrEt[i],   netbins);
    saveHistos(hAbsRspEt[i],   netbins);
    saveHistos(hRelRspEt[i],   netbins);

    saveHistos(hEta[i],        netabins);
    saveHistos(hEtCorrEta[i],  netabins);
    saveHistos(hAbsRspEta[i],  netabins);
    saveHistos(hRelRspEta[i],  netabins);

    saveHistos(hPhi[i],        nphibins);
    saveHistos(hEtCorrPhi[i],  nphibins);
    saveHistos(hAbsRspPhi[i],  nphibins);
    saveHistos(hRelRspPhi[i],  nphibins);
    
    saveHistos(hEtEtaEt[i],    netabins,netbins);
    saveHistos(hEtCorrEtaEt[i],netabins,netbins);
    saveHistos(hAbsRspEtaEt[i],netabins,netbins);
    saveHistos(hRelRspEtaEt[i],netabins,netbins);
  }
  
  plotfile->Write();
  plotfile->Close();
  delete plotfile;
  
  //app->Run();
  
  return 0;
}


////////////////////////////////////////////////////////////////////////////////
// implementation of global functions
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void replaceHistos(int nbins,TH1F**& hVar,TH1F**& hEtCorr,TTree**& tVar)
{
  TH1::SetDefaultSumw2();
  for (int ibin=0;ibin<nbins;ibin++) {

    TH1F* h(0);

    tVar[ibin]->Project("h(50)","val","weight*(1)");
    h = (TH1F*)gROOT->FindObject("h");
    h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
    h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
    h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
    delete hVar[ibin];
    hVar[ibin] = h;

    tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
    h = (TH1F*)gROOT->FindObject("h");
    h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
    h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
    h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
    delete hEtCorr[ibin];
    hEtCorr[ibin] = h;

    delete tVar[ibin];
  }
  
  delete [] tVar;
}

//______________________________________________________________________________
void fitHistos(TH1F** histos,int n,double nsigma)
{
  for (int i=0;i<n;i++) {
    float nrm  = histos[i]->Integral();
    float mean = histos[i]->GetMean();
    float rms  = histos[i]->GetRMS();
    TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
    f->SetParameter(0,nrm);
    f->SetParameter(1,mean);
    f->SetParameter(2,rms);
    f->SetLineWidth(1);
    f->SetLineColor(kRed);
    histos[i]->Fit(f,"QR");
  }
}


//______________________________________________________________________________
void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
{
  for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
}


//______________________________________________________________________________
void saveHistos(TH1F** histos,int n)
{
  for (int i=0;i<n;i++) histos[i]->Write();
}


//______________________________________________________________________________
void saveHistos(TH1F*** histos,int n1,int n2)
{
  for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
}
Revision:	1.3
Committed:	Thu Aug 23 07:01:26 2007 UTC (17 years, 8 months ago) by schiefer
Branch:	MAIN
Changes since 1.2:	+272 -677 lines
Log Message:	jes_mcresponse_x makes histos but no plots.
#	User	Rev	Content
1	schiefer	1.1	////////////////////////////////////////////////////////////////////////////////
2			//
3			// mcresponse_x
4			// ------------
5			//
6			// 06/16/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7			////////////////////////////////////////////////////////////////////////////////
8
9
10			#include "utils/cmdline.h"
11	schiefer	1.2	#include "utils/hist.h"
12	schiefer	1.1
13			#include <TROOT.h>
14			#include <TRint.h>
15			#include <TFile.h>
16			#include <TTree.h>
17			#include <TEventList.h>
18			#include <TH1F.h>
19			#include <TF1.h>
20			#include <TGraphErrors.h>
21			#include <TMultiGraph.h>
22			#include <TCanvas.h>
23			#include <TLegend.h>
24
25			#include <iostream>
26			#include <iomanip>
27			#include <fstream>
28			#include <cmath>
29			#include <string>
30			#include <vector>
31
32
33			using namespace std;
34
35
36			////////////////////////////////////////////////////////////////////////////////
37			// defaults
38			////////////////////////////////////////////////////////////////////////////////
39	schiefer	1.3	string dflt_etbins ="7,9,12,15,20,25,30,45,65,90,120";
40			string dflt_etabins=".2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.3,2.5,2.7,3.0,3.3";
41	schiefer	1.1	string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";
42
43
44			////////////////////////////////////////////////////////////////////////////////
45			// declare global functions
46			////////////////////////////////////////////////////////////////////////////////
47	schiefer	1.3	void replaceHistos(int nbins,TH1F& hVar,TH1F& hEtCorr,TTree**& tVar);
48			void fitHistos (TH1F** histos,int n,double nsigma);
49			void fitHistos (TH1F*** histos,int n1,int n2,double nsigma);
50			void saveHistos (TH1F** histos,int n);
51			void saveHistos (TH1F*** histos,int n1,int n2);
52	schiefer	1.1
53
54
55			////////////////////////////////////////////////////////////////////////////////
56			// main
57			////////////////////////////////////////////////////////////////////////////////
58
59			//______________________________________________________________________________
60			int main(int argc,char**argv)
61			{
62			cmdline cl;
63			cl.parse(argc,argv);
64	schiefer	1.2
65	schiefer	1.1	vector<string> input = cl.get_vector<string>("input");
66	schiefer	1.3	string output = cl.get_value<string> ("output","mcresponse.root");
67			string selection = cl.get_value<string> ("selection", "njt>0");
68			float drmax = cl.get_value<float> ("drmax", 0.3);
69			short applyjes = cl.get_value<short> ("applyjes", 0);
70			string treename = cl.get_value<string> ("treename", "t");
71			float nsigma = cl.get_value<float> ("nsigma", 3.0);
72			vector<float> etbins = cl.get_vector<float> ("etbins", dflt_etbins);
73			vector<float> etabins = cl.get_vector<float> ("etabins", dflt_etabins);
74			vector<float> phibins = cl.get_vector<float> ("phibins", dflt_phibins);
75			bool abseta = cl.get_value<bool> ("abseta", true);
76	schiefer	1.1
77			if (!cl.check()) return 0;
78			cl.print();
79
80
81			// etabins
82			if (!abseta) {
83			int neta=(int)etabins.size();
84			std::reverse(etabins.begin(),etabins.end());
85			for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
86			for (int ieta=0;ieta<neta;ieta++) etabins[ieta]*=-1;
87			}
88
89
90			// jetalgs
91			vector<string> jetalgs;
92			for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
93			string jetalg = *it;
94			string::size_type pos;
95			while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
96			jetalg = jetalg.substr(0,jetalg.find(".root"));
97			jetalgs.push_back(jetalg);
98			}
99
100			// prepare branch variables and histograms / graphs
101			float weight;
102			char njt;
103	schiefer	1.3	float jtet[100];
104	schiefer	1.1	float jtpt[100];
105			float jteta[100];
106			float jtphi[100];
107	schiefer	1.3	float jtgenet[100];
108	schiefer	1.1	float jtgenpt[100];
109			float jtgendr[100];
110			float jtjes[100][3];
111
112	schiefer	1.3	int netbins = etbins.size()+1;
113	schiefer	1.1	int netabins = etabins.size()+1;
114			int nphibins = phibins.size()+1;
115
116	schiefer	1.3	vector<TH1F**> hGenEt;
117			vector<TH1F**> hEtCorrGenEt;
118			vector<TH1F**> hAbsRspGenEt;
119			vector<TH1F**> hRelRspGenEt;
120
121			vector<TH1F**> hEt;
122			vector<TH1F**> hEtCorrEt;
123			vector<TH1F**> hAbsRspEt;
124			vector<TH1F**> hRelRspEt;
125
126	schiefer	1.1	vector<TH1F**> hEta;
127	schiefer	1.3	vector<TH1F**> hEtCorrEta;
128			vector<TH1F**> hAbsRspEta;
129			vector<TH1F**> hRelRspEta;
130
131	schiefer	1.1	vector<TH1F**> hPhi;
132	schiefer	1.3	vector<TH1F**> hEtCorrPhi;
133			vector<TH1F**> hAbsRspPhi;
134			vector<TH1F**> hRelRspPhi;
135
136			vector<TH1F***> hEtEtaEt;
137			vector<TH1F***> hEtCorrEtaEt;
138			vector<TH1F***> hAbsRspEtaEt;
139			vector<TH1F***> hRelRspEtaEt;
140
141			vector<TGraphErrors*> gRspVsGenEt;
142			vector<TGraphErrors*> gSgmVsGenEt;
143			vector<TGraphErrors*> gRspVsEt;
144			vector<TGraphErrors*> gSgmVsEt;
145	schiefer	1.1	vector<TGraphErrors*> gRspVsEta;
146			vector<TGraphErrors*> gSgmVsEta;
147			vector<TGraphErrors*> gRspVsPhi;
148			vector<TGraphErrors*> gSgmVsPhi;
149
150	schiefer	1.3	vector<TGraphErrors**> gRspVsEtaEt;
151			vector<TGraphErrors**> gSgmVsEtaEt;
152	schiefer	1.1
153	schiefer	1.3	argc=3; argv[1]="-l"; argv[2]="-b";
154	schiefer	1.1	TRint* app = new TRint(argv[0],&argc,argv);
155
156
157			//
158			// loop over all input files
159			//
160	schiefer	1.3	TFile* plotfile = new TFile(output.c_str(),"RECREATE");
161	schiefer	1.2
162	schiefer	1.1	for (unsigned int i=0;i<input.size();++i) {
163	schiefer	1.3	TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
164	schiefer	1.1	TTree* t = (TTree*)f->Get(treename.c_str());
165
166			TEventList* el = new TEventList("sel","sel");
167			t->Draw(">>sel",selection.c_str());
168
169			t->SetBranchAddress("weight", &weight);
170			t->SetBranchAddress("njt", &njt);
171	schiefer	1.3	t->SetBranchAddress("jtet", jtet);
172	schiefer	1.1	t->SetBranchAddress("jtpt", jtpt);
173			t->SetBranchAddress("jteta", jteta);
174			t->SetBranchAddress("jtphi", jtphi);
175	schiefer	1.3	t->SetBranchAddress("jtgenet",jtgenet);
176	schiefer	1.1	t->SetBranchAddress("jtgenpt",jtgenpt);
177			t->SetBranchAddress("jtgendr",jtgendr);
178
179			if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
180			else
181			for (int i1=0;i1<100;i1++)
182			for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
183
184
185	schiefer	1.3	string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
186			string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
187	schiefer	1.2
188	schiefer	1.3	hGenEt .push_back(hist::initHistos("GenEt",jetalgs[i],
189			100,"jet E_{T}^{gen} [GeV]",
190			"GenEt",etbins));
191			hEtCorrGenEt.push_back(hist::initHistos("EtCorrGenEt",jetalgs[i],
192			100,"jet E_{T}^{gen} [GeV]",
193			"GenEt",etbins));
194			hAbsRspGenEt.push_back(hist::initHistos("AbsRspGenEt",jetalgs[i],
195			100,-50,150,absrsp_xtitle,
196			"GenEt",etbins));
197			hRelRspGenEt.push_back(hist::initHistos("RelRspGenEt",jetalgs[i],
198			100,0,2,relrsp_xtitle,
199			"GenEt",etbins));
200
201			hEt .push_back(hist::initHistos("Et",jetalgs[i],
202			100,"jet E_{T} [GeV]",
203			"Et",etbins));
204			hEtCorrEt .push_back(hist::initHistos("EtCorrEt",jetalgs[i],
205			100,"jet E_{T} [GeV]",
206			"Et",etbins));
207			hAbsRspEt .push_back(hist::initHistos("AbsRspEt",jetalgs[i],
208			100,-50,150,absrsp_xtitle,
209			"Et",etbins));
210			hRelRspEt .push_back(hist::initHistos("RelRspEt",jetalgs[i],
211			100,0,2,relrsp_xtitle,
212			"Et",etbins));
213
214			hEta .push_back(hist::initHistos("Eta",jetalgs[i],
215			100,"jet #eta","Eta",etabins));
216			hEtCorrEta .push_back(hist::initHistos("EtCorrEta",jetalgs[i],
217			100,"jet E_{T} [GeV]",
218			"Eta",etabins));
219			hAbsRspEta .push_back(hist::initHistos("AbsRspEta",jetalgs[i],
220			100,-50,150,absrsp_xtitle,
221			"Eta",etabins));
222			hRelRspEta .push_back(hist::initHistos("RelRspEta",jetalgs[i],
223			100,0,2,relrsp_xtitle,
224	schiefer	1.2	"Eta",etabins));
225
226	schiefer	1.3	hPhi .push_back(hist::initHistos("Phi",jetalgs[i],
227			100,"jet #phi","Phi",phibins));
228			hEtCorrPhi .push_back(hist::initHistos("EtCorrPhi",jetalgs[i],
229			100,"jet E_{T} [Gev]",
230			"Phi",phibins));
231			hAbsRspPhi .push_back(hist::initHistos("AbsRspPhi",jetalgs[i],
232			100,-50,150,absrsp_xtitle,
233			"Phi",phibins));
234			hRelRspPhi .push_back(hist::initHistos("RelRspPhi",jetalgs[i],
235			100,0,2,relrsp_xtitle,
236	schiefer	1.2	"Phi",phibins));
237
238	schiefer	1.3	hEtEtaEt .push_back(hist::initHistos("Et",jetalgs[i],
239			100,"jet E_{T} [GeV]",
240			"Eta",etabins,"Et",etbins));
241			hEtCorrEtaEt.push_back(hist::initHistos("EtCorrEtaEt",jetalgs[i],
242			100,"jet E_{T} [GeV]",
243			"Eta",etabins,"Et",etbins));
244			hAbsRspEtaEt.push_back(hist::initHistos("AbsRspEtaEt",jetalgs[i],
245			100,-50,150,absrsp_xtitle,
246			"Eta",etabins,"Et",etbins));
247			hRelRspEtaEt.push_back(hist::initHistos("RelRspEtaEt",jetalgs[i],
248			100,0,2,relrsp_xtitle,
249			"Eta",etabins,"Et",etbins));
250	schiefer	1.1
251
252			int nevts = el->GetN();
253			cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
254	schiefer	1.3
255			// prepare intermediate trees
256			float genet,et,eta,phi,etcorr;
257			TTree** tGenEt = new TTree*[netbins];
258			TTree** tEt = new TTree*[netbins];
259			TTree** tEta = new TTree*[netabins];
260			TTree** tPhi = new TTree*[nphibins];
261			TTree* tEtEtaEt = new TTree[netabins];
262
263			for (int iet=0;iet<netbins;iet++) {
264			tGenEt[iet]=new TTree("tGenEt","tGenEt");
265			tGenEt[iet]->Branch("val", &genet, "val/F");
266			tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
267			tGenEt[iet]->Branch("weight",&weight,"weight/F");
268			tEt[iet]=new TTree("tEt","tEt");
269			tEt[iet]->Branch("val", &et, "val/F");
270			tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
271			tEt[iet]->Branch("weight",&weight,"weight/F");
272			}
273			for (int ieta=0;ieta<netabins;ieta++) {
274			tEta[ieta]=new TTree("tEta","tEta");
275			tEta[ieta]->Branch("val", &eta, "val/F");
276			tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
277			tEta[ieta]->Branch("weight",&weight,"weight/F");
278			}
279			for (int iphi=0;iphi<nphibins;iphi++) {
280			tPhi[iphi]=new TTree("tPhi","tPhi");
281			tPhi[iphi]->Branch("val", &phi, "val/F");
282			tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
283			tPhi[iphi]->Branch("weight",&weight,"weight/F");
284			}
285			for (int ieta=0;ieta<netabins;ieta++) {
286			tEtEtaEt[ieta] = new TTree*[netbins];
287			for (int iet=0;iet<netbins;iet++) {
288			tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
289			tEtEtaEt[ieta][iet]->Branch("val", &et, "val/F");
290			tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
291			tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
292			}
293			}
294	schiefer	1.1
295
296			// loop over all events
297			for (int ievt=0;ievt<nevts;ievt++) {
298
299			t->GetEntry(el->GetEntry(ievt));
300
301			// loop over all jets in the event
302			for (int ijt=0;ijt<njt;ijt++) {
303
304			if (jtgendr[ijt]<0.0\|\|jtgendr[ijt]>drmax) continue;
305
306			if (ijt>0) continue;
307
308	schiefer	1.3	et = jtet[ijt];
309			genet = jtgenet[ijt];
310			eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
311			phi = jtphi[ijt];
312			etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
313	schiefer	1.1
314	schiefer	1.3	float absrsp = genet - etcorr;
315			float relrsp = etcorr/genet;
316	schiefer	1.2
317	schiefer	1.3	// genet
318			int igenet = hist::get_ibin(genet,etbins);
319			hAbsRspGenEt[i][igenet] ->Fill(absrsp,weight);
320			hRelRspGenEt[i][igenet] ->Fill(relrsp,weight);
321			tGenEt[igenet]->Fill();
322	schiefer	1.1
323	schiefer	1.3	// et
324			int iet = hist::get_ibin(et,etbins);
325			hAbsRspEt[i][iet]->Fill(absrsp,weight);
326			hRelRspEt[i][iet]->Fill(relrsp,weight);
327			tEt[iet]->Fill();
328	schiefer	1.1
329			// eta
330	schiefer	1.2	int ieta = hist::get_ibin(eta,etabins);
331	schiefer	1.3	hAbsRspEta[i][ieta]->Fill(absrsp,weight);
332			hRelRspEta[i][ieta]->Fill(relrsp,weight);
333			tEta[ieta]->Fill();
334
335	schiefer	1.1	// phi
336	schiefer	1.2	int iphi = hist::get_ibin(phi,phibins);
337	schiefer	1.3	hAbsRspPhi[i][iphi]->Fill(absrsp,weight);
338			hRelRspPhi[i][iphi]->Fill(relrsp,weight);
339			tPhi[iphi]->Fill();
340
341			// et-eta
342			hAbsRspEtaEt[i][ieta][iet]->Fill(absrsp,weight);
343			hRelRspEtaEt[i][ieta][iet]->Fill(relrsp,weight);
344			tEtEtaEt[ieta][iet]->Fill();
345	schiefer	1.1
346			} // jets
347			} // evts
348
349			cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
350	schiefer	1.3
351	schiefer	1.2
352	schiefer	1.3	// replace histograms
353			replaceHistos(netbins,hGenEt[i],hEtCorrGenEt[i],tGenEt);
354			replaceHistos(netbins,hEt[i], hEtCorrEt[i], tEt);
355			replaceHistos(netabins,hEta[i], hEtCorrEta[i], tEta);
356			replaceHistos(nphibins,hPhi[i], hEtCorrPhi[i], tPhi);
357			for (int ieta=0;ieta<netabins;ieta++)
358			replaceHistos(netbins,hEtEtaEt[i][ieta],hEtCorrEtaEt[i][ieta],tEtEtaEt[ieta]);
359	schiefer	1.2
360			TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
361			d->cd();
362	schiefer	1.1
363	schiefer	1.3	// fit response histos
364			fitHistos(hAbsRspGenEt[i],netbins,nsigma);
365			fitHistos(hRelRspGenEt[i],netbins,nsigma);
366
367			fitHistos(hAbsRspEt[i], netbins,nsigma);
368			fitHistos(hRelRspEt[i], netbins,nsigma);
369
370			fitHistos(hAbsRspEta[i], netabins,nsigma);
371			fitHistos(hRelRspEta[i], netabins,nsigma);
372
373			fitHistos(hAbsRspPhi[i], nphibins,nsigma);
374			fitHistos(hRelRspPhi[i], nphibins,nsigma);
375
376			fitHistos(hAbsRspEtaEt[i],netabins,netbins,nsigma);
377			fitHistos(hRelRspEtaEt[i],netabins,netbins,nsigma);
378
379
380			// save histograms
381			saveHistos(hGenEt[i], netbins);
382			saveHistos(hEtCorrGenEt[i],netbins);
383			saveHistos(hAbsRspGenEt[i],netbins);
384			saveHistos(hRelRspGenEt[i],netbins);
385
386			saveHistos(hEt[i], netbins);
387			saveHistos(hEtCorrEt[i], netbins);
388			saveHistos(hAbsRspEt[i], netbins);
389			saveHistos(hRelRspEt[i], netbins);
390
391			saveHistos(hEta[i], netabins);
392			saveHistos(hEtCorrEta[i], netabins);
393			saveHistos(hAbsRspEta[i], netabins);
394			saveHistos(hRelRspEta[i], netabins);
395
396			saveHistos(hPhi[i], nphibins);
397			saveHistos(hEtCorrPhi[i], nphibins);
398			saveHistos(hAbsRspPhi[i], nphibins);
399			saveHistos(hRelRspPhi[i], nphibins);
400
401			saveHistos(hEtEtaEt[i], netabins,netbins);
402			saveHistos(hEtCorrEtaEt[i],netabins,netbins);
403			saveHistos(hAbsRspEtaEt[i],netabins,netbins);
404			saveHistos(hRelRspEtaEt[i],netabins,netbins);
405	schiefer	1.1	}
406
407	schiefer	1.2	plotfile->Write();
408			plotfile->Close();
409			delete plotfile;
410	schiefer	1.1
411	schiefer	1.3	//app->Run();
412	schiefer	1.1
413			return 0;
414			}
415
416
417
418			////////////////////////////////////////////////////////////////////////////////
419			// implementation of global functions
420			////////////////////////////////////////////////////////////////////////////////
421
422			//______________________________________________________________________________
423	schiefer	1.3	void replaceHistos(int nbins,TH1F& hVar,TH1F& hEtCorr,TTree**& tVar)
424	schiefer	1.1	{
425	schiefer	1.3	TH1::SetDefaultSumw2();
426	schiefer	1.1	for (int ibin=0;ibin<nbins;ibin++) {
427
428	schiefer	1.3	TH1F* h(0);
429	schiefer	1.1
430	schiefer	1.3	tVar[ibin]->Project("h(50)","val","weight*(1)");
431			h = (TH1F*)gROOT->FindObject("h");
432			h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
433			h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
434			h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
435			delete hVar[ibin];
436			hVar[ibin] = h;
437
438			tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
439			h = (TH1F*)gROOT->FindObject("h");
440			h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
441			h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
442			h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
443			delete hEtCorr[ibin];
444			hEtCorr[ibin] = h;
445	schiefer	1.1
446	schiefer	1.3	delete tVar[ibin];
447	schiefer	1.1	}
448
449	schiefer	1.3	delete [] tVar;
450	schiefer	1.1	}
451
452			//______________________________________________________________________________
453	schiefer	1.3	void fitHistos(TH1F** histos,int n,double nsigma)
454	schiefer	1.1	{
455	schiefer	1.3	for (int i=0;i<n;i++) {
456			float nrm = histos[i]->Integral();
457			float mean = histos[i]->GetMean();
458			float rms = histos[i]->GetRMS();
459			TF1* f = new TF1("gaus","gaus",mean-nsigmarms,mean+nsigmarms);
460			f->SetParameter(0,nrm);
461			f->SetParameter(1,mean);
462			f->SetParameter(2,rms);
463			f->SetLineWidth(1);
464			f->SetLineColor(kRed);
465			histos[i]->Fit(f,"QR");
466	schiefer	1.1	}
467			}
468
469
470			//______________________________________________________________________________
471	schiefer	1.3	void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
472	schiefer	1.1	{
473	schiefer	1.3	for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
474	schiefer	1.1	}
475
476
477			//______________________________________________________________________________
478	schiefer	1.3	void saveHistos(TH1F** histos,int n)
479	schiefer	1.1	{
480	schiefer	1.3	for (int i=0;i<n;i++) histos[i]->Write();
481	schiefer	1.1	}
482
483
484			//______________________________________________________________________________
485	schiefer	1.3	void saveHistos(TH1F*** histos,int n1,int n2)
486	schiefer	1.1	{
487	schiefer	1.3	for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
488	schiefer	1.1	}