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.3 by jengbou, Mon Apr 20 23:41:31 2009 UTC

# Line 15 | Line 15
15   #include <TPaveStats.h>
16   #include <TTree.h>
17  
18 < TString fileName = "QCDbin4.root";
18 > // Output file
19 > TString fileName = "skimmed/QCDbinall_NewRelIso.root";
20 > bool FastSim = false;
21 > TString isoname = "New";
22 > // Jet bin:
23 > Int_t jbin = 0; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets
24 > double binwidth = 0.01; // 0.02,0.02,0.05,0.1 for New RelIso
25   // Histo range
26   double xMin = 0.0;
27 < 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;
27 > double xMax = 0.0;
28   // Specify fitting function:
27 TString fitfunc = "landau"; // this version is cut for Landau
29   // Select New or Old reliso:
29 TString isoname = "New";
30   double d0sigCut = 3.;
31  
32   Int_t fjbin;
33   Double_t fna; // Num of QCD events in signal region
34  
35 + int nbin;
36 +
37   void skim::Loop()
38   {
39    //   In a ROOT session, you can do:
# Line 42 | Line 44 | void skim::Loop()
44    //      Root > t.Show();       // Show values of entry 12
45    //      Root > t.Show(16);     // Read and show values of entry 16
46    //      Root > t.Loop();       // Loop on all entries
47 +
48 +  if ( isoname.Contains("Old") ) {
49 +    xMax = 1+binwidth;
50 +    nbin = (1.0-xMin)/binwidth+1;
51 +  }
52 +  else if ( isoname.Contains("New") ) {
53 +    xMax = 2.0;
54 +    nbin = (xMax-xMin)/binwidth;
55 +  }
56 +  else {
57 +    cout << "No Rel Iso defined!" << endl;
58 +    return;
59 +  }
60 +  cout << "nbin for histo = " << nbin << endl;
61    TTree *ftree = new TTree("myTree","myTree");
62    ftree->AutoSave();
63    ftree->Branch("jbin",&fjbin,"jbin/I");
# Line 58 | Line 74 | void skim::Loop()
74    // define fit region:
75    double AX0 = 0.;
76    double AXf = 0.;
61  double BX0 = 0.;
62  double BXf = 0.;
77    
78    cout<<"Use "<<isoname<<" RelIso"<<endl;
79    
# Line 67 | Line 81 | void skim::Loop()
81      // old
82      AX0 = 0.95;
83      AXf = 1.0;
70    BX0 = 0.2;
71    BXf = 0.8;
84    }
85    else if ( isoname.Contains("New") ) {
86      // new
87      AX0 = 0.0;
88 <    AXf = 0.053;
77 <    BX0 = 0.4;
78 <    BXf = 2.0;
88 >    AXf = 0.053; // default:0.053 ; loose:0.11
89    }
90    
91 <  // weights
92 <  double wttbar = 0.0101;
93 <  double wqcd   = 1.3161;
94 <  double wwjets = 0.0977;
95 <  double wzjets = 0.078;
91 >  // Weights:
92 >  double wttbar,wqcd,wwjets,wzjets;
93 >  wttbar = wqcd = wwjets = wzjets = 0.;
94 >  if ( ! FastSim ) {
95 >    cout << "Full sim sample weights used " << endl;
96 >    // Full Sim
97 >    wttbar = 0.0101;
98 >    wqcd   = 1.3161;
99 >    wwjets = 0.0977;
100 >    wzjets = 0.078;
101 >  }
102 >  else {
103 >    // Fast Sim
104 >    cout << "Fast sim sample weights used " << endl;
105 >    wttbar = 0.0008;
106 >    wqcd   = 0.676;
107 >    wwjets = 0.0091;
108 >    wzjets = 0.0094;
109 >  }
110    
111    // Na = num of qcd in signal region
112    double Na;
# Line 103 | Line 127 | void skim::Loop()
127      nb = fChain->GetEntry(jentry);   nbytes += nb;
128      // if (Cut(ientry) < 0) continue;
129      bool signalRegion = false;
106    bool fitRegion = false;
130      
131      TString filename(fChain->GetCurrentFile()->GetName());
132      if (isoname.Contains("Old")) {
133        reliso = top_muon_old_reliso[0];
134        if ( reliso > AX0 ) signalRegion = true;
112      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
135      }
136      else if (isoname.Contains("New")) {
137        reliso = top_muon_new_reliso[0];
138        if ( reliso < AXf ) signalRegion = true;
117      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
139      }
140      else {
141        cout<<"No RelIso defined!!!"<<endl;
142        return;
143      }
144 <    
144 >
145      double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
146      bool correctBin = false; // jet bin condition
147      bool goodD0sig = false; // d0sig condition
148  
149      if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
150 <    if (top_njets == jbin && jbin != 4) correctBin = true;
151 <    else if (top_njets >= jbin && jbin == 4) correctBin = true;
150 >    if ( jbin == 0 ){
151 >      if ( top_njets > jbin ) { correctBin = true; }
152 >    }
153 >    else if ( jbin == 4 ){
154 >      if ( top_njets >= jbin ) { correctBin = true; }
155 >    }
156 >    else {
157 >      if ( top_njets == jbin ) { correctBin = true; }
158 >    }
159 >    
160      if ( d0sig < d0sigCut ) goodD0sig = true;
161      
162 <    // Jet bin
162 >    // Counting
163      if( correctBin && goodD0sig ) {
164 <      if ( filename.Contains("TTJets") ) {
164 >      if ( filename.Contains("TT") ) {
165          h_RelIso_ttbar->Fill(reliso);
166          if ( signalRegion ) { Nttbar += wttbar; }
167        }
168 <      if ( filename.Contains("MuPt15") ) {
168 >      if ( filename.Contains("Mu") ) {
169          h_RelIso_qcd->Fill(reliso);
170          if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
171        }
# Line 161 | Line 190 | void skim::Loop()
190  
191    TFile *file0 = new TFile(fileName,"RECREATE");
192    file0->mkdir("histos");
193 +  file0->mkdir("hstack");
194    file0->cd("histos");  
195    TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);  
196    THStack *hs = new THStack("hs","RelIso (stacked)");
# Line 224 | Line 254 | void skim::Loop()
254      return;
255    }
256    h_RelIso_all->SetMinimum(0);
257 <
257 >  
258 >  file0->cd("hstack");
259 >  hs->Write();
260 >  
261    file0->cd();
262    ftree->Write();
230 //   file0->Write();
263    file0->Close();
264   }
265  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines