ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/skim.C
(Generate patch)

Comparing UserCode/Jeng/scripts/skim.C (file contents):
Revision 1.1 by jengbou, Wed Apr 15 06:35:28 2009 UTC vs.
Revision 1.5 by jengbou, Wed May 6 02:38:07 2009 UTC

# Line 15 | Line 15
15   #include <TPaveStats.h>
16   #include <TTree.h>
17  
18 < TString fileName = "QCDbin4.root";
18 > // Output file
19 > // New
20 > // TString fileName = "skimmed226/QCDbin123_NewRelIso_0053.root";
21 > // TString fileName = "skimmed226/QCDbin4_NewRelIso_rebin3_TandL.root";
22 > // TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL.root";
23 > // Old
24 > TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin3_wErrors_TandL.root";
25 > // TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin1_FastSim_095.root";
26 >
27 > bool FastSim = false;
28 > // Select New or Old reliso:
29 > TString isoname = "Old";
30 > // Jet bin:
31 > Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets, 5 for 1-3
32 > double binwidth = 0.005; // 0.02,0.02,0.05,0.1 for New RelIso
33 > // Which region (true for CD or false for AB)
34 > bool regionCD = false;
35 > double d0sigCut = 3.; // 3.
36   // Histo range
37   double xMin = 0.0;
38 < double xMax = 2.0;
22 < // Jet bin:
23 < Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
24 < double binwidth = 0.1; // 0.02,0.02,0.05,0.1
25 < int nbin = (xMax-xMin)/binwidth;
26 < // Specify fitting function:
27 < TString fitfunc = "landau"; // this version is cut for Landau
28 < // Select New or Old reliso:
29 < TString isoname = "New";
30 < double d0sigCut = 3.;
38 > double xMax = 0.0;
39  
40   Int_t fjbin;
41   Double_t fna; // Num of QCD events in signal region
42 + Double_t fna1; // Num of QCD events in signal region
43 +
44 + int nbin;
45  
46   void skim::Loop()
47   {
# Line 42 | Line 53 | void skim::Loop()
53    //      Root > t.Show();       // Show values of entry 12
54    //      Root > t.Show(16);     // Read and show values of entry 16
55    //      Root > t.Loop();       // Loop on all entries
56 +
57 +  if ( isoname.Contains("Old") ) {
58 +    xMax = 1+binwidth;
59 +    nbin = (1.0-xMin)/binwidth+1;
60 +  }
61 +  else if ( isoname.Contains("New") ) {
62 +    xMax = 2.0;
63 +    nbin = (xMax-xMin)/binwidth;
64 +  }
65 +  else {
66 +    cout << "No Rel Iso defined!" << endl;
67 +    return;
68 +  }
69 +  cout << "nbin for histo = " << nbin << endl;
70 +  if ( regionCD ) cout << "==> Selecting Region CD." << endl;
71 +  
72    TTree *ftree = new TTree("myTree","myTree");
73    ftree->AutoSave();
74    ftree->Branch("jbin",&fjbin,"jbin/I");
75    ftree->Branch("na",&fna,"na/D");
76 +  ftree->Branch("na1",&fna1,"na1/D");
77    ftree->SetBranchAddress("jbin",&fjbin);
78    ftree->SetBranchAddress("na",&fna);
79 +  ftree->SetBranchAddress("na1",&fna1);
80    
81    TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
82    TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
83    TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
84    TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
85    TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
86 +  
87 +  h_RelIso_ttbar->Sumw2();
88 +  h_RelIso_qcd->Sumw2();
89 +  h_RelIso_wjets->Sumw2();
90 +  h_RelIso_zjets->Sumw2();
91      
92    // define fit region:
93    double AX0 = 0.;
94    double AXf = 0.;
95 <  double BX0 = 0.;
62 <  double BXf = 0.;
95 >  double LooseCut = 0.;
96    
97    cout<<"Use "<<isoname<<" RelIso"<<endl;
98    
# Line 67 | Line 100 | void skim::Loop()
100      // old
101      AX0 = 0.95;
102      AXf = 1.0;
103 <    BX0 = 0.2;
71 <    BXf = 0.8;
103 >    LooseCut = 0.9;
104    }
105    else if ( isoname.Contains("New") ) {
106      // new
107      AX0 = 0.0;
108 <    AXf = 0.053;
109 <    BX0 = 0.4;
78 <    BXf = 2.0;
108 >    AXf = 0.053; // Tight: 0.053
109 >    LooseCut = 0.11;  // Loose: 0.11
110    }
111    
112 <  // weights
113 <  double wttbar = 0.0101;
114 <  double wqcd   = 1.3161;
115 <  double wwjets = 0.0977;
116 <  double wzjets = 0.078;
112 >  // Weights:
113 >  double wttbar,wqcd,wwjets,wzjets;
114 >  wttbar = wqcd = wwjets = wzjets = 0.;
115 >  if ( ! FastSim ) {
116 >    cout << "Full sim sample weights used " << endl;
117 >    // Full Sim
118 >    wttbar = 0.0101;
119 >    wqcd   = 0.3907;
120 >    wwjets = 0.0883;
121 >    wzjets = 0.0744;
122 >  }
123 >  else {
124 >    // Fast Sim
125 >    cout << "Fast sim sample weights used " << endl;
126 >    wttbar = 0.0008;
127 >    wqcd   = 0.676;
128 >    wwjets = 0.0091;
129 >    wzjets = 0.0094;
130 >  }
131    
132    // Na = num of qcd in signal region
133 <  double Na;
134 <  Na = 0;
133 >  double Na[2] = {0.,0.};
134 >  double Nttbar[2] = {0.,0.};
135 >  double NWjets[2] = {0.,0.};
136 >  double NZjets[2] = {0.,0.};
137 >  double Nqcd[2] = {0.,0.};
138    
91 –  double Nttbar, NWjets, NZjets, Nqcd;
92 –  Nttbar = NWjets = NZjets = Nqcd = 0;
139    double reliso = 0.0;
140    
141    if (fChain == 0) return;
# Line 102 | Line 148 | void skim::Loop()
148      if (ientry < 0) break;
149      nb = fChain->GetEntry(jentry);   nbytes += nb;
150      // if (Cut(ientry) < 0) continue;
151 <    bool signalRegion = false;
152 <    bool fitRegion = false;
151 >    bool signalRegion1 = false; // Is tight
152 >    bool signalRegion2 = false; // Is loose
153      
154      TString filename(fChain->GetCurrentFile()->GetName());
155      if (isoname.Contains("Old")) {
156        reliso = top_muon_old_reliso[0];
157 <      if ( reliso > AX0 ) signalRegion = true;
158 <      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
157 >      if ( reliso > AX0 ) signalRegion1 = true;
158 >      if ( reliso > LooseCut ) signalRegion2 = true;
159      }
160      else if (isoname.Contains("New")) {
161        reliso = top_muon_new_reliso[0];
162 <      if ( reliso < AXf ) signalRegion = true;
163 <      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
162 >      if ( reliso < AXf ) signalRegion1 = true;
163 >      if ( reliso < LooseCut ) signalRegion2 = true;
164      }
165      else {
166        cout<<"No RelIso defined!!!"<<endl;
167        return;
168      }
169 <    
169 >
170      double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
171      bool correctBin = false; // jet bin condition
172      bool goodD0sig = false; // d0sig condition
127 –
128 –    if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
129 –    if (top_njets == jbin && jbin != 4) correctBin = true;
130 –    else if (top_njets >= jbin && jbin == 4) correctBin = true;
131 –    if ( d0sig < d0sigCut ) goodD0sig = true;
173      
174 <    // Jet bin
174 >    if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
175 >    if ( jbin == 0 ){
176 >      if ( top_njets > jbin ) { correctBin = true; }
177 >    }
178 >    else if ( jbin == 5 ){
179 >      if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
180 >    }
181 >    else if ( jbin == 4 ){
182 >      if ( top_njets >= jbin ) { correctBin = true; }
183 >    }
184 >    else {
185 >      if ( top_njets == jbin ) { correctBin = true; }
186 >    }
187 >    
188 >    if ( !regionCD && (d0sig < d0sigCut) ) goodD0sig = true;
189 >    if ( regionCD && (d0sig > d0sigCut) ) goodD0sig = true;
190 >    
191 >    // Counting
192      if( correctBin && goodD0sig ) {
193 <      if ( filename.Contains("TTJets") ) {
193 >      if ( filename.Contains("TT") ) {
194          h_RelIso_ttbar->Fill(reliso);
195 <        if ( signalRegion ) { Nttbar += wttbar; }
195 >        if ( signalRegion1 ) { Nttbar[0] += wttbar; }
196 >        if ( signalRegion2 ) { Nttbar[1] += wttbar; }
197        }
198 <      if ( filename.Contains("MuPt15") ) {
198 >      if ( filename.Contains("Mu") ) {
199          h_RelIso_qcd->Fill(reliso);
200 <        if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
200 >        if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
201 >        if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
202        }
203        if ( filename.Contains("WJets") ) {
204          h_RelIso_wjets->Fill(reliso);
205 <        if ( signalRegion ) { NWjets += wwjets; }      
205 >        if ( signalRegion1 ) { NWjets[0] += wwjets; }  
206 >        if ( signalRegion2 ) { NWjets[1] += wwjets; }  
207        }
208        if ( filename.Contains("ZJets") ) {
209          h_RelIso_zjets->Fill(reliso);
210 <        if ( signalRegion ) { NZjets += wzjets; }
210 >        if ( signalRegion1 ) { NZjets[0] += wzjets; }
211 >        if ( signalRegion2 ) { NZjets[1] += wzjets; }
212        }
213      }
214    }
215    fjbin=jbin;
216 <  fna=Na;
216 >  fna=Na[0];
217 >  fna1=Na[1];
218    ftree->Fill();
219    cout << "\n";
220 <  cout << " N ttbar = " << Nttbar << endl;
221 <  cout << " N qcd = " << Nqcd << endl;
222 <  cout << " N Wjets = " << NWjets << endl;
223 <  cout << " N Zjets = " << NZjets << endl << endl;
220 >  if (isoname.Contains("New"))
221 >    cout << " Tight Cut: " << AXf << endl;
222 >  else if (isoname.Contains("Old"))
223 >    cout << " Tight Cut: " << AX0 << endl;
224 >  cout << " N ttbar = " << Nttbar[0] << endl;
225 >  cout << " N qcd   = " << Nqcd[0] << endl;
226 >  cout << " N Wjets = " << NWjets[0] << endl;
227 >  cout << " N Zjets = " << NZjets[0] << endl << endl;
228 >  cout << "=========================" << endl;
229 >  cout << "\n";
230 >  cout << " Loose Cut: " << LooseCut << endl;
231 >  cout << " N ttbar = " << Nttbar[1] << endl;
232 >  cout << " N qcd   = " << Nqcd[1] << endl;
233 >  cout << " N Wjets = " << NWjets[1] << endl;
234 >  cout << " N Zjets = " << NZjets[1] << endl << endl;
235  
236    TFile *file0 = new TFile(fileName,"RECREATE");
237    file0->mkdir("histos");
238 +  file0->mkdir("hstack");
239    file0->cd("histos");  
240    TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);  
241    THStack *hs = new THStack("hs","RelIso (stacked)");
# Line 193 | Line 268 | void skim::Loop()
268    h_RelIso_ttbar->SetMarkerColor(2);
269    h_RelIso_ttbar->SetLineColor(2);
270    h_RelIso_ttbar->SetFillColor(2);
271 +  h_RelIso_ttbar->GetXaxis()->SetTitle("RelIso");
272    hs->Add(h_RelIso_ttbar);
273    hs->Draw();
274 +  if ( isoname.Contains("New") )
275 +    hs->GetXaxis()->SetTitle("RelIso'");
276 +  else if ( isoname.Contains("Old") )
277 +    hs->GetXaxis()->SetTitle("RelIso");
278 +  hs->GetYaxis()->SetTitle("Events");
279 +  hs->GetXaxis()->SetLabelFont(62);
280 +  hs->GetXaxis()->SetTitleFont(22);
281 +  hs->GetXaxis()->SetTitleSize(0.04);
282 +  hs->GetYaxis()->SetLabelFont(62);
283 +  hs->GetYaxis()->SetTitleFont(62);
284 +  hs->GetYaxis()->SetTitleSize(0.04);
285  
286    // Label histo
287 <  TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55);
287 >  TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
288 >  leg0->SetFillColor(0);
289 >  leg0->SetFillStyle(0);
290    leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
291    leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
292    leg0->AddEntry(h_RelIso_wjets,"W+jets");
# Line 219 | Line 308 | void skim::Loop()
308      title += " jets)";
309      h_RelIso_all->SetTitle(title);
310    }
311 +  else if ( jbin == 5 ) {
312 +    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
313 +  }
314    else {
315 <    cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;
315 >    cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
316      return;
317    }
318    h_RelIso_all->SetMinimum(0);
319 <
319 >  
320 >  file0->cd("hstack");
321 >  hs->Write();
322 >  
323    file0->cd();
324    ftree->Write();
230 – //   file0->Write();
325    file0->Close();
326   }
327  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines