Jeng/scripts/skim_data.c

#define skim_data_cxx
#include "skim_data.h"
#include <TH2.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"
#include <iostream>
#include <TLegend.h>
#include <TPaveStats.h>
#include <TTree.h>
#include <cmath>
#include "map.h"

std::map<TString,TString> histName;
std::map<TString,TH1*> hs;
std::map<TString,TCanvas*> ts;

// Output file
// New
TString fileName = "skimmed361p3_20100621/Data_Incl_NewRelIso_wErrors_TandL_26.root";
//TString fileName = "skimmed361p3_20100621/MC_Incl_NewRelIso_wErrors_TandL_2.root";
//TString fileName = "skimmed361p3_20100621/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_26.root";
//TString fileName = "skimmed361p3_20100621/MC_Wjets_Incl_NewRelIso_wErrors_TandL_26.root";

// Old
//TString fileName = "skimmed361p3_20100621/Data_Incl_OldRelIso_wErrors_TandL_26.root";
//TString fileName = "skimmed361p3_20100621/MC_Incl_OldRelIso_wErrors_TandL_2.root";
//TString fileName = "skimmed361p3_20100621/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_26.root";
//TString fileName = "skimmed361p3_20100621/MC_Wjets_Incl_OldRelIso_wErrors_TandL_26.root";

// Select New or Old reliso:
TString isoname = "New";
// Jet bin:
Int_t jbin = -1; // 0-3, 4 for >=4. Put in number less than 5; -1 for all jets (inclusive), 5 for 1-3
double binwidth = 0.001; // 0.02,0.02,0.05,0.1 for New RelIso
// Which region (true for CD or false for AB)
bool regionCD = false;

// d0sigCut and d0Cut == -1 means not used
double d0sigCut = -1; // 3.
double d0Cut = 0.02; // 200 microns
// Histo range
double xMin = 0.0;
double xMax = 0.0;

Int_t fjbin;
Double_t fna; // Num of QCD events in signal region
Double_t fna1; // Num of QCD events in signal region

int nbin;

void skim_data::Loop()
{
  //   In a ROOT session, you can do:
  //      Modify jbin, binwidth, wxmin, wxmax
  //      Root > .L skim_data.C++
  //      Root > skim_data t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries

  if ( isoname.Contains("Old") ) {
    xMax = 1+binwidth;
    nbin = (1.0-xMin)/binwidth+1;
  }
  else if ( isoname.Contains("New") ) {
    xMax = 2.0;
    nbin = (xMax-xMin)/binwidth;
  }
  else {
    cout << "No Rel Iso defined!" << endl;
    return;
  }
  cout << "nbin for histo = " << nbin << endl;
  if ( regionCD ) cout << "==> Selecting Region CD." << endl;
  
  TTree *ftree = new TTree("myTree","myTree");
  ftree->AutoSave();
  ftree->Branch("jbin",&fjbin,"jbin/I");
  ftree->Branch("na",&fna,"na/D");
  ftree->Branch("na1",&fna1,"na1/D");
  ftree->SetBranchAddress("jbin",&fjbin);
  ftree->SetBranchAddress("na",&fna);
  ftree->SetBranchAddress("na1",&fna1);
  
  TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
  TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
  
  h_RelIso_qcd->Sumw2();
    
  // define fit region:
  double AX0 = 0.;
  double AXf = 0.;
  double LooseCut = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
    LooseCut = 0.9;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 1./19.; // Tight: 1/19
    LooseCut = 1./9.;  // Loose: 1/9
  }
  
  // Weights:
  double wqcd = 1.;

  // Na = num of qcd in signal region
  double Na[2] = {0.,0.};
  double Nqcd[2] = {0.,0.};
  
  double reliso = 0.0;
  
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntries();

  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    bool signalRegion1 = false; // Is tight
    bool signalRegion2 = false; // Is loose
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso;
      if ( reliso > AX0 ) signalRegion1 = true;
      if ( reliso > LooseCut ) signalRegion2 = true;
    }
    else if (isoname.Contains("New")) {
      /*       reliso = top_muon_new_reliso; */
      reliso = 1./top_muon_old_reliso - 1.;
      if ( reliso < AXf ) signalRegion1 = true;
      if ( reliso < LooseCut ) signalRegion2 = true;
    }
    else {
      cout<<"No RelIso defined!!!"<<endl;
      return;
    }

    double d0 = TMath::Abs(top_muon_d0);
    double d0sig = TMath::Abs(top_muon_d0)/top_muon_d0Error;
    bool correctBin = false; // jet bin condition
    bool goodD0sig = false; // d0sig condition
    
    if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}

    if ( jbin == -1 ) { correctBin = true; }
    else if ( jbin == 0 ){
      if ( top_njets == jbin ) { correctBin = true; }
    }
    else if ( jbin == 5 ){
      if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
    }
    else if ( jbin == 4 ){
      if ( top_njets >= jbin ) { correctBin = true; }
    }
    else {
      if ( top_njets == jbin ) { correctBin = true; }
    }
    
    if ( !regionCD && ((d0sigCut != -1 && d0sig < d0sigCut) || (d0Cut != -1 && d0 < d0Cut)) ) goodD0sig = true;
    if ( regionCD && ((d0sigCut != -1 && d0sig > d0sigCut) || (d0Cut != -1 && d0 > d0Cut)) ) goodD0sig = true;
    
    // Counting
    if( correctBin && goodD0sig
        && top_muon_pt > 20
        //&& top_trk_muon == 1
        //&& top_met > 5
        ) {
      h_RelIso_qcd->Fill(reliso);
      if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
      if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
    }
  }
  fjbin=jbin;
  fna=Na[0];
  fna1=Na[1];
  ftree->Fill();
  cout << "\n";
  if (isoname.Contains("New"))
    cout << " Tight Cut: " << AXf << endl;
  else if (isoname.Contains("Old"))
    cout << " Tight Cut: " << AX0 << endl;
  cout << " N events   = " << round(Nqcd[0],0) << endl;
  cout << "=========================" << endl;
  cout << "\n";
  cout << " Loose Cut: " << LooseCut << endl;
  cout << " N events   = " << round(Nqcd[1],0) << endl << endl;

  TFile *file0 = new TFile(fileName,"RECREATE");
  file0->mkdir("histos");
  file0->mkdir("hstack");
  file0->cd("histos");  
  TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);  
  THStack *hs = new THStack("hs","RelIso (stacked)");
  
  // Take care of scaling
  h_RelIso_qcd->Scale(wqcd);

  // Save for later use
  h_RelIso_qcd->Write();
  
  h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
  h_RelIso_qcd->SetLineColor(40);
  h_RelIso_qcd->SetFillStyle(3018);
  h_RelIso_qcd->SetFillColor(38);
  h_RelIso_qcd->GetXaxis()->SetTitle("RelIso");
  hs->Add(h_RelIso_qcd);
  hs->Draw();
  if ( isoname.Contains("New") )
    hs->GetXaxis()->SetTitle("RelIso'");
  else if ( isoname.Contains("Old") )
    hs->GetXaxis()->SetTitle("RelIso");
  hs->GetYaxis()->SetTitle("Events");
  hs->GetXaxis()->SetLabelFont(62);
  hs->GetXaxis()->SetTitleFont(22);
  hs->GetXaxis()->SetTitleSize(0.04);
  hs->GetYaxis()->SetLabelFont(62);
  hs->GetYaxis()->SetTitleFont(62);
  hs->GetYaxis()->SetTitleSize(0.04);

  // Label histo
  TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
  leg0->SetFillColor(0);
  leg0->SetFillStyle(0);
  leg0->AddEntry(h_RelIso_qcd,"Data");
  leg0->Draw();
  
  h_RelIso_all->Add(h_RelIso_qcd,1.);
  h_RelIso_all->Write();
  h_RelIso_all->SetLineWidth(2);
  
  if ( jbin == -1 )
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (inclusive)");
  else if ( jbin == 4 )
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
  else if ( jbin < 4 ) {
    TString title = isoname+" RelIso distribution (";
    title += jbin;
    title += " jets)";
    h_RelIso_all->SetTitle(title);
  }
  else if ( jbin == 5 ) {
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
  }
  else {
    cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
    return;
  }
  h_RelIso_all->SetMinimum(0);
  
  file0->cd("hstack");
  hs->Write();
  
  file0->cd();
  ftree->Write();
  file0->Close();
}

double skim_data::round(double num, int places)
{
    double temp, mult;
    mult = pow(10.0, places);
    temp = floor(num * mult + 0.5);
    temp = temp / mult;
    return temp;
}

Revision:	1.1
Committed:	Tue Jun 22 22:37:55 2010 UTC (14 years, 10 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Log Message:	Add file
#	User	Rev	Content
1	jengbou	1.1	#define skim_data_cxx
2			#include "skim_data.h"
3			#include <TH2.h>
4			#include <THStack.h>
5			#include <TStyle.h>
6			#include <TCanvas.h>
7			#include <TMath.h>
8			#include <TMultiGraph.h>
9			#include <TGraphErrors.h>
10			#include <TCanvas.h>
11			#include "Math/GSLIntegrator.h"
12			#include "Math/WrappedTF1.h"
13			#include <iostream>
14			#include <TLegend.h>
15			#include <TPaveStats.h>
16			#include <TTree.h>
17			#include <cmath>
18			#include "map.h"
19
20			std::map<TString,TString> histName;
21			std::map<TString,TH1*> hs;
22			std::map<TString,TCanvas*> ts;
23
24			// Output file
25			// New
26			TString fileName = "skimmed361p3_20100621/Data_Incl_NewRelIso_wErrors_TandL_26.root";
27			//TString fileName = "skimmed361p3_20100621/MC_Incl_NewRelIso_wErrors_TandL_2.root";
28			//TString fileName = "skimmed361p3_20100621/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_26.root";
29			//TString fileName = "skimmed361p3_20100621/MC_Wjets_Incl_NewRelIso_wErrors_TandL_26.root";
30
31			// Old
32			//TString fileName = "skimmed361p3_20100621/Data_Incl_OldRelIso_wErrors_TandL_26.root";
33			//TString fileName = "skimmed361p3_20100621/MC_Incl_OldRelIso_wErrors_TandL_2.root";
34			//TString fileName = "skimmed361p3_20100621/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_26.root";
35			//TString fileName = "skimmed361p3_20100621/MC_Wjets_Incl_OldRelIso_wErrors_TandL_26.root";
36
37			// Select New or Old reliso:
38			TString isoname = "New";
39			// Jet bin:
40			Int_t jbin = -1; // 0-3, 4 for >=4. Put in number less than 5; -1 for all jets (inclusive), 5 for 1-3
41			double binwidth = 0.001; // 0.02,0.02,0.05,0.1 for New RelIso
42			// Which region (true for CD or false for AB)
43			bool regionCD = false;
44
45			// d0sigCut and d0Cut == -1 means not used
46			double d0sigCut = -1; // 3.
47			double d0Cut = 0.02; // 200 microns
48			// Histo range
49			double xMin = 0.0;
50			double xMax = 0.0;
51
52			Int_t fjbin;
53			Double_t fna; // Num of QCD events in signal region
54			Double_t fna1; // Num of QCD events in signal region
55
56			int nbin;
57
58			void skim_data::Loop()
59			{
60			// In a ROOT session, you can do:
61			// Modify jbin, binwidth, wxmin, wxmax
62			// Root > .L skim_data.C++
63			// Root > skim_data t
64			// Root > t.GetEntry(12); // Fill t data members with entry number 12
65			// Root > t.Show(); // Show values of entry 12
66			// Root > t.Show(16); // Read and show values of entry 16
67			// Root > t.Loop(); // Loop on all entries
68
69			if ( isoname.Contains("Old") ) {
70			xMax = 1+binwidth;
71			nbin = (1.0-xMin)/binwidth+1;
72			}
73			else if ( isoname.Contains("New") ) {
74			xMax = 2.0;
75			nbin = (xMax-xMin)/binwidth;
76			}
77			else {
78			cout << "No Rel Iso defined!" << endl;
79			return;
80			}
81			cout << "nbin for histo = " << nbin << endl;
82			if ( regionCD ) cout << "==> Selecting Region CD." << endl;
83
84			TTree *ftree = new TTree("myTree","myTree");
85			ftree->AutoSave();
86			ftree->Branch("jbin",&fjbin,"jbin/I");
87			ftree->Branch("na",&fna,"na/D");
88			ftree->Branch("na1",&fna1,"na1/D");
89			ftree->SetBranchAddress("jbin",&fjbin);
90			ftree->SetBranchAddress("na",&fna);
91			ftree->SetBranchAddress("na1",&fna1);
92
93			TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
94			TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
95
96			h_RelIso_qcd->Sumw2();
97
98			// define fit region:
99			double AX0 = 0.;
100			double AXf = 0.;
101			double LooseCut = 0.;
102
103			cout<<"Use "<<isoname<<" RelIso"<<endl;
104
105			if ( isoname.Contains("Old") ) {
106			// old
107			AX0 = 0.95;
108			AXf = 1.0;
109			LooseCut = 0.9;
110			}
111			else if ( isoname.Contains("New") ) {
112			// new
113			AX0 = 0.0;
114			AXf = 1./19.; // Tight: 1/19
115			LooseCut = 1./9.; // Loose: 1/9
116			}
117
118			// Weights:
119			double wqcd = 1.;
120
121			// Na = num of qcd in signal region
122			double Na[2] = {0.,0.};
123			double Nqcd[2] = {0.,0.};
124
125			double reliso = 0.0;
126
127			if (fChain == 0) return;
128
129			Long64_t nentries = fChain->GetEntries();
130
131			Long64_t nbytes = 0, nb = 0;
132			for (Long64_t jentry=0; jentry<nentries;jentry++) {
133			Long64_t ientry = LoadTree(jentry);
134			if (ientry < 0) break;
135			nb = fChain->GetEntry(jentry); nbytes += nb;
136			// if (Cut(ientry) < 0) continue;
137			bool signalRegion1 = false; // Is tight
138			bool signalRegion2 = false; // Is loose
139
140			TString filename(fChain->GetCurrentFile()->GetName());
141			if (isoname.Contains("Old")) {
142			reliso = top_muon_old_reliso;
143			if ( reliso > AX0 ) signalRegion1 = true;
144			if ( reliso > LooseCut ) signalRegion2 = true;
145			}
146			else if (isoname.Contains("New")) {
147			/* reliso = top_muon_new_reliso; */
148			reliso = 1./top_muon_old_reliso - 1.;
149			if ( reliso < AXf ) signalRegion1 = true;
150			if ( reliso < LooseCut ) signalRegion2 = true;
151			}
152			else {
153			cout<<"No RelIso defined!!!"<<endl;
154			return;
155			}
156
157			double d0 = TMath::Abs(top_muon_d0);
158			double d0sig = TMath::Abs(top_muon_d0)/top_muon_d0Error;
159			bool correctBin = false; // jet bin condition
160			bool goodD0sig = false; // d0sig condition
161
162			if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
163
164			if ( jbin == -1 ) { correctBin = true; }
165			else if ( jbin == 0 ){
166			if ( top_njets == jbin ) { correctBin = true; }
167			}
168			else if ( jbin == 5 ){
169			if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
170			}
171			else if ( jbin == 4 ){
172			if ( top_njets >= jbin ) { correctBin = true; }
173			}
174			else {
175			if ( top_njets == jbin ) { correctBin = true; }
176			}
177
178			if ( !regionCD && ((d0sigCut != -1 && d0sig < d0sigCut) \|\| (d0Cut != -1 && d0 < d0Cut)) ) goodD0sig = true;
179			if ( regionCD && ((d0sigCut != -1 && d0sig > d0sigCut) \|\| (d0Cut != -1 && d0 > d0Cut)) ) goodD0sig = true;
180
181			// Counting
182			if( correctBin && goodD0sig
183			&& top_muon_pt > 20
184			//&& top_trk_muon == 1
185			//&& top_met > 5
186			) {
187			h_RelIso_qcd->Fill(reliso);
188			if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
189			if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
190			}
191			}
192			fjbin=jbin;
193			fna=Na[0];
194			fna1=Na[1];
195			ftree->Fill();
196			cout << "\n";
197			if (isoname.Contains("New"))
198			cout << " Tight Cut: " << AXf << endl;
199			else if (isoname.Contains("Old"))
200			cout << " Tight Cut: " << AX0 << endl;
201			cout << " N events = " << round(Nqcd[0],0) << endl;
202			cout << "=========================" << endl;
203			cout << "\n";
204			cout << " Loose Cut: " << LooseCut << endl;
205			cout << " N events = " << round(Nqcd[1],0) << endl << endl;
206
207			TFile *file0 = new TFile(fileName,"RECREATE");
208			file0->mkdir("histos");
209			file0->mkdir("hstack");
210			file0->cd("histos");
211			TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
212			THStack *hs = new THStack("hs","RelIso (stacked)");
213
214			// Take care of scaling
215			h_RelIso_qcd->Scale(wqcd);
216
217			// Save for later use
218			h_RelIso_qcd->Write();
219
220			h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
221			h_RelIso_qcd->SetLineColor(40);
222			h_RelIso_qcd->SetFillStyle(3018);
223			h_RelIso_qcd->SetFillColor(38);
224			h_RelIso_qcd->GetXaxis()->SetTitle("RelIso");
225			hs->Add(h_RelIso_qcd);
226			hs->Draw();
227			if ( isoname.Contains("New") )
228			hs->GetXaxis()->SetTitle("RelIso'");
229			else if ( isoname.Contains("Old") )
230			hs->GetXaxis()->SetTitle("RelIso");
231			hs->GetYaxis()->SetTitle("Events");
232			hs->GetXaxis()->SetLabelFont(62);
233			hs->GetXaxis()->SetTitleFont(22);
234			hs->GetXaxis()->SetTitleSize(0.04);
235			hs->GetYaxis()->SetLabelFont(62);
236			hs->GetYaxis()->SetTitleFont(62);
237			hs->GetYaxis()->SetTitleSize(0.04);
238
239			// Label histo
240			TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
241			leg0->SetFillColor(0);
242			leg0->SetFillStyle(0);
243			leg0->AddEntry(h_RelIso_qcd,"Data");
244			leg0->Draw();
245
246			h_RelIso_all->Add(h_RelIso_qcd,1.);
247			h_RelIso_all->Write();
248			h_RelIso_all->SetLineWidth(2);
249
250			if ( jbin == -1 )
251			h_RelIso_all->SetTitle( isoname+" RelIso distribution (inclusive)");
252			else if ( jbin == 4 )
253			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
254			else if ( jbin < 4 ) {
255			TString title = isoname+" RelIso distribution (";
256			title += jbin;
257			title += " jets)";
258			h_RelIso_all->SetTitle(title);
259			}
260			else if ( jbin == 5 ) {
261			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
262			}
263			else {
264			cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
265			return;
266			}
267			h_RelIso_all->SetMinimum(0);
268
269			file0->cd("hstack");
270			hs->Write();
271
272			file0->cd();
273			ftree->Write();
274			file0->Close();
275			}
276
277			double skim_data::round(double num, int places)
278			{
279			double temp, mult;
280			mult = pow(10.0, places);
281			temp = floor(num * mult + 0.5);
282			temp = temp / mult;
283			return temp;
284			}
285