ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/plotSignalFakes.cc
(Generate patch)

Comparing UserCode/MitHzz4l/NonMCBackground/src/plotSignalFakes.cc (file contents):
Revision 1.4 by dkralph, Tue Jul 31 16:46:22 2012 UTC vs.
Revision 1.6 by dkralph, Tue Oct 23 11:22:41 2012 UTC

# Line 9 | Line 9
9   #include "TStyle.h"
10   #include "TH1D.h"
11   #include "TNtuple.h"
12 + #include "TROOT.h"
13  
14   #include "CPlot.h"
15   #include "FOArgs.h"
# Line 20 | Line 21
21   #include "TriggerUtils.h"
22   #include "JetDefs.h"
23   #include "JetInfoStruct.h"
24 + #include "FR_struct.h"
25 + #include "SelectionFuncs.h"
26 + #include "SampleWeight.h"
27 + #include "MitStyleRemix.h"
28  
29   #include "TMVA/Reader.h"
30   #include "TMVA/Tools.h"
# Line 34 | Line 39 | using namespace std;
39   using namespace RooFit;
40   using namespace mithep;
41  
42 < TH1D* hpu_2011;
43 < TH1D* hpu_2012;
44 < //----------------------------------------------------------------------------------------
45 < class EvtsOfInterest
41 < {
42 < public:
43 <  EvtsOfInterest(TString config="");
44 <  bool has(unsigned run, unsigned evt);
45 <  vector<unsigned> runs,evts;
46 < };
47 < EvtsOfInterest::EvtsOfInterest(TString config)
48 < {
49 <  if(config!="") {
50 <    ifstream ifs(config);
51 <    assert(ifs.is_open());
52 <    string line;
53 <    while(getline(ifs,line)) {
54 <      stringstream ss(line);
55 <      unsigned run,evt;
56 <      ss >> run >> evt;
57 <      runs.push_back(run);
58 <      evts.push_back(evt);
59 <    }
60 <    ifs.close();
61 <  }
62 < }
63 < bool EvtsOfInterest::has(unsigned run, unsigned evt)
64 < {
65 <  for(unsigned i=0; i<runs.size(); i++) {
66 <    if(runs[i]==run && evts[i]==evt)
67 <      return true;
68 <  }
69 <  return false;
70 < }
42 > // TH1D* hpu_2011;
43 > // TH1D* hpu_2012;
44 > // TH1D* hpu_2011_me;
45 > // TH1D* hpu_2012_me;
46   //----------------------------------------------------------------------------------------
47   TCanvas *can;
73 //----------------------------------------------------------------------------------------
74 class lepFrs
75 {
76 public:
77  double fwgt_3,fwgt_4;
78  double fwgt_lo_3,fwgt_lo_4;
79  double fwgt_hi_3,fwgt_hi_4;
80 };
81 //----------------------------------------------------------------------------------------
82 class dr_struct
83 {
84 public:
85  dr_struct();
86  void check();
87  float l1j1,l2j1,l3j1,l4j1;
88  float l1j2,l2j2,l3j2,l4j2;
89  float l1j3,l2j3,l3j3,l4j3;
90  float l1j4,l2j4,l3j4,l4j4;
91 };
92 dr_struct::dr_struct()
93 {  
94  l1j1 = l2j1 = l3j1 = l4j1 = l1j2 = l2j2 = l3j2 = l4j2 = l1j3 = l2j3 = l3j3 = l4j3 = l1j4 = l2j4 = l3j4 = l4j4 = -1;
95 }
96 void dr_struct::check()
97 {
98  if(fabs(l1j1) < 0.2) { cout << "l1j1 too small: " << l1j1 << endl; }
99  if(fabs(l2j1) < 0.2) { cout << "l2j1 too small: " << l2j1 << endl; }
100  if(fabs(l3j1) < 0.2) { cout << "l3j1 too small: " << l3j1 << endl; }
101  if(fabs(l4j1) < 0.2) { cout << "l4j1 too small: " << l4j1 << endl; }
102  if(fabs(l1j2) < 0.2) { cout << "l1j2 too small: " << l1j2 << endl; }
103  if(fabs(l2j2) < 0.2) { cout << "l2j2 too small: " << l2j2 << endl; }
104  if(fabs(l3j2) < 0.2) { cout << "l3j2 too small: " << l3j2 << endl; }
105  if(fabs(l4j2) < 0.2) { cout << "l4j2 too small: " << l4j2 << endl; }
106  if(fabs(l1j3) < 0.2) { cout << "l1j3 too small: " << l1j3 << endl; }
107  if(fabs(l2j3) < 0.2) { cout << "l2j3 too small: " << l2j3 << endl; }
108  if(fabs(l3j3) < 0.2) { cout << "l3j3 too small: " << l3j3 << endl; }
109  if(fabs(l4j3) < 0.2) { cout << "l4j3 too small: " << l4j3 << endl; }
110  if(fabs(l1j4) < 0.2) { cout << "l1j4 too small: " << l1j4 << endl; }
111  if(fabs(l2j4) < 0.2) { cout << "l2j4 too small: " << l2j4 << endl; }
112  if(fabs(l3j4) < 0.2) { cout << "l3j4 too small: " << l3j4 << endl; }
113  if(fabs(l4j4) < 0.2) { cout << "l4j4 too small: " << l4j4 << endl; }
114 }  
48  
49 + TString dummy_integral_str(TH1D *hist, int nDecPlaces) { return TString(""); };
50   bool findGoodJets(vector<SimpleLepton> &goodJets, filestuff *fs,
51                    SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4);
52   void fillHist(CSample *cs, TString channel, TString type, TString var, double val, double wgt, double wgt_lo=0, double wgt_hi=0);
53 < void fillAllHists(FOFlags ctrl, CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine,
53 > void fillAllHists(FOFlags ctrl, CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine, Angles angles,
54                    SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4,
55                    double wgt, double wgt_lo=0, double wgt_hi=0);
56   void fillAllJetHists(FOFlags ctrl,  CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine,
57                        SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4,
58                        FusionMva &fusion, vector<SimpleLepton> &goodJets, JetInfoStruct ji,
59                        double wgt, double wgt_lo=0, double wgt_hi=0);
60 < void makeHTML(FOFlags &ctrl, TString type, TString plotLabel);
60 > void makeHTML(FOFlags &ctrl, TString type, TString plotLabel, TString fullOutDir);
61   map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
62   bool passHlt(FOFlags &ctrl, TrigInfo ti, InfoStruct *info, unsigned lep1matchBits,
63               unsigned lep2matchBits, unsigned lep3matchBits,
64               unsigned lep4matchBits);
131 dr_struct fill_dr_struct(vector<SimpleLepton> jets, SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4);
65   void init_cuts(vector<TString> &cutstrs, map<TString,int> &cutvec);
133 void fillLeptonFakeWeights(lepFrs &fwgts, FR_struct *fr, SimpleLepton lep3, SimpleLepton lep4);
66   void resetCutVect(map<TString,int> &cutvec) {
67    for(map<TString,int>::iterator it=cutvec.begin(); it!=cutvec.end(); it++)
68      cutvec[(*it).first] = 0;
69   }
70   map<TString,TH1D*> setHistSet(CSample *cs, TString type, TString var);
139 void shiftOverflows(map<TString,TH1D*> &hset);
140 void scaleHists(map<TString,TH1D*> &hset);
71   void fillFakeTuple(CSample *cs, filestuff *fs, EventData evtdata, unsigned ientry, double fillweight);
72 + void incrementSsofCounters(FOFlags &ctrl, filestuff *fs, double minMz2,
73 +                           pair<int,int> *best_z_indices, double wgt);
74 + map<TString,double> init_counter(TString fname="");
75 + map<TString,double> initSsofRatios(vector<CSample*> &samplev);
76 + double getSsofWgt(TString channel, map<TString,double> &ssofRatios);
77   //----------------------------------------------------------------------------------------
78   int main(int argc, char** argv)
79   {
80 +  SetStyle();
81    initPUWeights();
82 +  ControlFlags dummyCtrl;
83 +  // TFile * puf_me;
84 +  // puf_me = new TFile("data/pileup/PUWeights_S12To2012_190456-199011.root");
85 +  // hpu_2012_me = (TH1D*)(puf_me->Get("puWeights"));
86 +  // hpu_2012_me->SetDirectory(0);
87 +  // puf_me->Close();
88 +  // puf_me = new TFile("data/pileup/PUWeights_F11To2011.root");
89 +  // hpu_2011_me = (TH1D*)(puf_me->Get("puWeights"));
90 +  // hpu_2011_me->SetDirectory(0);
91 +  // puf_me->Close();
92 +
93    vector<TString> cutstrs;
94    map<TString,int> cutvec;
95    init_cuts(cutstrs, cutvec);
96  
150  EvtsOfInterest eoiAdish("/tmp/adish-uniq");
151  EvtsOfInterest eoiMe("/tmp/me-uniq");
152
97    // arguments...
98    FOFlags ctrl;
99    parse_foargs( argc, argv, ctrl );
100    ctrl.dump();
157  bool makeFakeTuples = (ctrl.extraArgs.Contains("makeFakeTuples"));
101    FusionMva fusion(ctrl.uncert);//"./weights/againstZZ-fusion_BDTG.weights.xml");
102    bool makeJetTuple=false; // jet tuple for vbf mva training
103    assert(ctrl.faketype=="SR"   || // plot events in signal region, if requested including ntuples of the extrapolation from 2P2F and 3P1F regions
# Line 163 | Line 106 | int main(int argc, char** argv)
106    if(ctrl.heavyFlavor) assert(ctrl.faketype=="2P2F");
107    if(ctrl.ssof) assert(ctrl.faketype=="2P2F");
108  
109 <  can = new TCanvas("can","can");
109 >  can = new TCanvas("can","can",700,500);
110  
111 <  FR_struct fr2011 = initFRs(ctrl.mufakefile2011,ctrl.elefakefile2011);
112 <  FR_struct fr2012 = initFRs(ctrl.mufakefile2012,ctrl.elefakefile2012);
111 >  FR_struct fr2011(ctrl.mufakefile2011,ctrl.elefakefile2011);
112 >  FR_struct fr2012(ctrl.mufakefile2012,ctrl.elefakefile2012);
113  
114    TString cmsswpath(CMSSW_BASE + TString("/src"));
115    string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
# Line 179 | Line 122 | int main(int argc, char** argv)
122    TString ntupledir(""),label(""),plotLabel(""),jsonFile("");
123    int puTarget;
124    readConfigFile(ctrl.config, ntupledir, label, plotLabel, jsonFile, puTarget, samplev, &ctrl, init_hists);
125 +  map<TString,double> ssofRatios;
126 +  if(ctrl.ssof && !ctrl.writessofratio && !ctrl.plotWholeSample) ssofRatios = initSsofRatios(samplev);
127  
128 +  Angles angles; // NOTE: this is not the same as fs->angles -- the one in fs is filled with whichever choice of leptons was made in applyZPlusX
129    KinematicsStruct kine; // NOTE: this is not the same as fs->kine -- the one in fs is filled with whichever choice of leptons was made in applyZPlusX
130    JetInfoStruct ji;
131  
# Line 193 | Line 139 | int main(int argc, char** argv)
139      AngleTuple *jettuple;
140      if(makeJetTuple) {
141        jettuple = new AngleTuple( (const char*)(TString(cs->name)+"_pass.root").Data(), (const char*)"zznt");
142 +      jettuple->makeAngleBranch(angles);
143        jettuple->makeKinematicsBranch(kine);
144        jettuple->makeJetInfoBranch(ji);
145        const char *str1 = "evt";
# Line 204 | Line 151 | int main(int argc, char** argv)
151      resetCutVect(cutvec);
152      for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
153        filestuff *fs = (cs->fsv)[ifs];
154 <      if(makeFakeTuples && fs->dataset_!="fakes") { // write a copy of the input trees for events that will be used for the SR fake prediction
154 >      if(ctrl.makeFakeTuples && fs->dataset_!="fakes") { // write a copy of the input trees for events that will be used for the SR fake prediction
155          TString fakefilename(fs->fname_);
156          fakefilename.ReplaceAll(".root","-fakes.root");
157 +        if(ctrl.hiStatFakes) fakefilename.ReplaceAll("-fakes.root","-fakes-histat.root");
158          fs->makeOutputFile(fakefilename);
159        }
160        cout << "\t" << fs->fname_; cout.flush();
213      FR_struct *fr = (fs->era_==2011) ? &fr2011 : &fr2012;
161        unsigned nDuplSkipped=0;
162        for(unsigned ientry=0; ientry<fs->getentries("FOtree"); ientry++) {
163          fs->getentry(ientry,"FOs","FOtree");
164          fs->getentry(ientry,"Zleptons","FOtree");
165          fs->getentry(ientry,"info","zznt");
166  
220        // if(eoiAdish.has(fs->info->run, fs->info->evt) &&
221        //    !eoiMe.has(fs->info->run, fs->info->evt)) {
222        //   cout << "found an event: " << fs->info->run << " " << fs->info->evt << endl;
223        //   ctrl.debug = true;
224        // } else
225        //   ctrl.debug = false;
226
167          cutvec["start"] = fs->total_entries_;
168          if(fs->isdata_) {// && !(fs->dataset_=="fakes")) {
169            // WARNING: when I write the fakeTuple from here, I fill an event *twice* if it has two fake z2 candidates, so event numbers are duplicated.
170            // in plotH4l I *don't* want to remove these, but here I *do*, because in this file I again loop over all z2 candidates
171            setMinMaxRun(fs->info->run, minRun, maxRun);
172            bool dupl = takeCareOfDuplicateEvents(fs->info->run, fs->info->evt, runEvtv, nDuplSkipped);
173 <          if(dupl) { if(ctrl.debug) cout << " failing dupl" << endl; continue; }
174 <          if(!fs->rlrm_.HasRunLumi(pair<unsigned,unsigned> (fs->info->run, fs->info->lumi))) { if(ctrl.debug) cout << " failing rlrm" << endl; continue; }
173 >          if(dupl) continue;
174 >          if(!fs->rlrm_.HasRunLumi(pair<unsigned,unsigned> (fs->info->run, fs->info->lumi))) continue;
175          }
176          cutvec["rlrmAndDupl"] += 1;
177  
178          double wgt=1;
179          if(!fs->isdata_) {
180            double xsWgt = fs->lumi_*xstab.Get(fs->dataset_)/fs->total_entries_;
181 <          double puWgt = getPUWeight(fs->era_,fs->info->npu);
182 <          if(fs->era_==2012)
183 <            puWgt = weightTrue2012(fs->info->npu);
181 >          cout << "WARNING: check the pu weihgs funciton here" << endl;
182 >          double puWgt = getPUWeight(fs->era_, "/029/", fs->info->npu);
183 >          // if(fs->era_==2012)
184 >          //   puWgt = hpu_2012_me->GetBinContent(hpu_2012_me->FindBin(fs->info->npu));
185 >          // else if(fs->era_==2011)
186 >          //   puWgt = hpu_2011_me->GetBinContent(hpu_2011_me->FindBin(fs->info->npu));
187 >          // else assert(0);
188            wgt = xsWgt*puWgt;
189 +          if(wgt!=wgt) cout << "xsWgt: " << xsWgt << " puWgt: " << puWgt << endl;
190          }
191          if(fs->dataset_=="fakes") { // if we're reading in the fake ntuple that we wrote out earlier, use the event weight that was stored (this should be the same as what we'd get if we calculated it below)
192            fs->getentry(ientry,"weights","zznt");
# Line 257 | Line 202 | int main(int argc, char** argv)
202          if(best_mz1<=0) {
203            cout << "WARNING: best z1 mass: " << best_mz1 << endl;
204          }
205 <        if(ctrl.heavyFlavor && best_mz1<60) { if(ctrl.debug) cout << " failing mz1" << endl; continue; }
205 >        if(ctrl.heavyFlavor && best_mz1<60) continue;
206          EventData evtdata;
207          evtdata.Z1leptons.push_back((*fs->passingL)[best_z_indices.first]);
208          evtdata.Z1leptons.push_back((*fs->passingL)[best_z_indices.second]);
# Line 266 | Line 211 | int main(int argc, char** argv)
211          vector<pair<SimpleLepton,SimpleLepton> > z2cands;
212          vector<TString> types;
213          TString signFlavor;
214 <        if(ctrl.heavyFlavor) signFlavor = "";
215 <        else if(ctrl.ssof)   signFlavor = "SS-SF";
216 <        else                 signFlavor = "OS-SF";
214 >        if(ctrl.heavyFlavor)            signFlavor = "";
215 >        else if(ctrl.ssof)              signFlavor = "SS-SF";
216 >        else if(ctrl.hiStatFakes)       signFlavor = "signalExclude";
217 >        else                            signFlavor = "OS-SF";
218          double minMz2 = ctrl.heavyFlavor ? 0 : 12;
219          int hiPtPassingZ2 = findZ2CandidatesPassFail(fs->passingL, z2cands, types, signFlavor, minMz2, "pass",&best_z_indices); // signal region
220          int hiPtFailingZ2 = findZ2CandidatesPassFail(fs->failingL, z2cands, types, signFlavor, minMz2, "fail"); // 2P+2F
221          int hiPtPfZ2      = findZ2CandidatesPassFail(fs->passingL, z2cands, types, signFlavor, minMz2, "PF", &best_z_indices, fs->failingL); // 3P+1F
222  
223 <        if(ctrl.debug) {
278 <          cout << "z2 cands: " << endl;
279 <          for(unsigned iz2=0; iz2<z2cands.size(); iz2++)
280 <            cout << " ---" << types[iz2] << "'ing z2" << endl;
281 <        }
223 >        if(ctrl.writessofratio) incrementSsofCounters(ctrl, fs, minMz2, &best_z_indices, wgt);
224  
225          for(unsigned iz2=0; iz2<z2cands.size(); iz2++) {
226            SimpleLepton lep1 = evtdata.Z1leptons[0];
# Line 288 | Line 230 | int main(int argc, char** argv)
230            evtdata.Z2leptons.clear();
231            evtdata.Z2leptons.push_back(lep3);
232            evtdata.Z2leptons.push_back(lep4);
233 +          fillAngles(evtdata,angles);
234            fillKinematics(evtdata,kine);
235            lepFrs fwgts;
236 <          fillLeptonFakeWeights(fwgts, fr, lep3, lep4);
236 >          fillLeptonFakeWeights(fwgts, (fs->era_==2011) ? &fr2011 : &fr2012, lep3, lep4);
237            TString channel = getChannel(lep1,lep2,lep3,lep4);
238            // bool hltPass = fs->isdata_ ? passHlt(ctrl,ti,fs->info,lep1.scID,lep2.scID,lep3.scID,lep4.scID) : true;
239  
240            if(ctrl.heavyFlavor) {
241 <            if(fabs(lep1.dz) > 0.01 || fabs(lep2.dz) > 0.01)            { cout << " failing XXX" << endl; continue; }
241 >            if(fabs(lep1.dz) > 0.01 || fabs(lep2.dz) > 0.01)            continue;
242              if( !(fabs(lep3.ip3dSig) > 2 && fabs(lep4.ip3dSig) > 8) &&
243 <                !(fabs(lep4.ip3dSig) > 2 && fabs(lep3.ip3dSig) > 8)  )  { cout << " failing XXX" << endl; continue; }
243 >                !(fabs(lep4.ip3dSig) > 2 && fabs(lep3.ip3dSig) > 8)  )  continue;
244              if( !(fabs(lep3.d0) > .002 && fabs(lep4.d0) > .01) &&
245 <                !(fabs(lep4.d0) > .002 && fabs(lep3.d0) > .01)  )       { cout << " failing XXX" << endl; continue; }
246 <            if(fabs(lep3.d0) > 0.1 || fabs(lep4.d0) > 0.1)              { cout << " failing XXX" << endl; continue; }
247 <            if(fabs(lep3.dz) > 0.2 || fabs(lep4.dz) > 0.2)              { cout << " failing XXX" << endl; continue; }
245 >                !(fabs(lep4.d0) > .002 && fabs(lep3.d0) > .01)  )       continue;
246 >            if(fabs(lep3.d0) > 0.1 || fabs(lep4.d0) > 0.1)              continue;
247 >            if(fabs(lep3.dz) > 0.2 || fabs(lep4.dz) > 0.2)              continue;
248            } else {
249 <            if(fabs(lep3.ip3dSig) >= 4 || fabs(lep4.ip3dSig) >= 4)      { if(ctrl.debug) cout << " failing sip" << endl; continue; } // |sip| < 4 is applied to lep1, lep2 in ZPlusX
250 <            if(!leptonDrReqs(lep1,lep2,lep3,lep4))                      { if(ctrl.debug) cout << " failing dr" << endl; continue; }
251 <            if(kine.m4l < 100)                                          { if(ctrl.debug) cout << " failing m4l 100" << endl; continue; }
252 <            if(!resonanceKilling(lep1,lep2,lep3,lep4))                  { if(ctrl.debug) cout << " failing reskill" << endl; continue; }
249 >            if(fabs(lep3.ip3dSig) >= 4 || fabs(lep4.ip3dSig) >= 4)      continue; // |sip| < 4 is applied to lep1, lep2 in ZPlusX
250 >            if(!leptonDrReqs(lep1,lep2,lep3,lep4))                      continue;
251 >            if(kine.m4l < 100)                                          continue;
252 >            if(!resonanceKilling(lep1,lep2,lep3,lep4))                  continue;
253            }
254 <          if(!finalLeptonPtReqs(lep1,lep2,lep3,lep4))                   { if(ctrl.debug) cout << " failing lepton pt " << endl; continue; }
254 >          if(!finalLeptonPtReqs(lep1,lep2,lep3,lep4))                   continue;
255  
256            cutvec["4lsele"] += 1;
257  
# Line 322 | Line 265 | int main(int argc, char** argv)
265              if(fs->jets->size() > 1) cutvec["twoJets"] += 1;
266              bool hazJetz = findGoodJets(goodJets, fs, lep1, lep2, lep3, lep4);
267              if(goodJets.size() > 1) cutvec["twoJetsAfter"] += 1;
268 <            if(!hazJetz) { cout << " failing XXX" << endl; continue; }
268 >            if(!hazJetz) continue;
269              ControlFlags ctrlTmp;
270              fillJetInfo(goodJets,ji,ctrlTmp);
271              if(makeJetTuple) {
# Line 331 | Line 274 | int main(int argc, char** argv)
274              }
275              fusion.setValues(ji,kine);
276            }
277 <
277 >          
278            cutvec["filling"] += 1;
279  
280            if(types[iz2]=="pass" || ((fs->dataset_=="fakes") && (types[iz2]=="fail"))) { // for the fakes ntuple, we plot the "failing" z2s as passing ones, i.e. put them in hObs, but with the proper weight
281              // llll events
282 <            if((fs->dataset_!="fakes") && (iz2!=hiPtPassingZ2)) { cout << " failing XXX" << endl; continue; } // for non-fakes, we only want the  highest-pt "passing" z2 (for the fakes these are by design failing z2s, and we want all of 'em)
282 >            if((fs->dataset_!="fakes") && (iz2!=hiPtPassingZ2)) continue; // for non-fakes, we only want the  highest-pt "passing" z2 (for the fakes these are by design failing z2s, and we want all of 'em)
283              cutvec["fillPass"] += 1;
284 <            fillAllHists( ctrl, cs, channel, "obs", fs, kine, lep1, lep2, lep3, lep4, wgt);
284 >            fillAllHists( ctrl, cs, channel, "obs", fs, kine, angles, lep1, lep2, lep3, lep4, wgt);
285              if(doJets) fillAllJetHists( ctrl, cs, channel, "obs", fs, kine, lep1, lep2, lep3, lep4, fusion, goodJets, ji, wgt);
286  
287            } else if(types[iz2]=="fail") { // lljj events
288 <            if(iz2!=hiPtFailingZ2 && ctrl.plotWholeSample) { cout << " failing XXX" << endl; continue; }
288 >            if(iz2!=hiPtFailingZ2 && ctrl.plotWholeSample) continue;
289              double centr=0,lo=0,hi=0;
290              if(ctrl.faketype=="2P2F") {
291                if(ctrl.plotWholeSample) { // make plots without weighting with the FR
# Line 351 | Line 294 | int main(int argc, char** argv)
294                  centr = wgt*fwgts.fwgt_3*fwgts.fwgt_4;
295                  lo    = wgt*fwgts.fwgt_lo_3*fwgts.fwgt_lo_4;
296                  hi    = wgt*fwgts.fwgt_hi_3*fwgts.fwgt_hi_4;
297 +                if(ctrl.ssof && !ctrl.writessofratio && !ctrl.plotWholeSample) {
298 +                  double ssof_wgt = getSsofWgt(channel,ssofRatios);
299 +                  centr *= ssof_wgt;
300 +                  lo    *= ssof_wgt;
301 +                  hi    *= ssof_wgt;
302 +                }
303                }
304              } else if(ctrl.faketype=="3P1F") { // apply the weights for extrapolation from 2P2F region to 3P1F region
305                centr = wgt*(fwgts.fwgt_3    + fwgts.fwgt_4);
# Line 358 | Line 307 | int main(int argc, char** argv)
307                hi    = wgt*(fwgts.fwgt_hi_3 + fwgts.fwgt_hi_4);
308              }
309  
310 <            // runfile << fixed << setprecision(0) << fs->info->run << " " << fixed << setprecision(0) << fs->info->evt << endl;
362 <
363 <            if(makeFakeTuples && fs->dataset_!="fakes") {
310 >            if(ctrl.makeFakeTuples && fs->dataset_!="fakes") {
311                double fillwgt;
312                bool dofill=true;
313                if(ctrl.faketype=="2P2F") {
# Line 373 | Line 320 | int main(int argc, char** argv)
320                fillFakeTuple(cs,fs,evtdata,ientry,fillwgt);
321              }
322              cutvec["fillFail"] += 1;
323 <            fillAllHists( ctrl, cs, channel, "pred", fs, kine, lep1, lep2, lep3, lep4, centr, lo, hi);
323 >            fillAllHists( ctrl, cs, channel, "pred", fs, kine, angles, lep1, lep2, lep3, lep4, centr, lo, hi);
324              if(doJets) fillAllJetHists( ctrl, cs, channel, "pred", fs, kine, lep1, lep2, lep3, lep4, fusion, goodJets, ji, centr, lo, hi);
325  
326            } else if(types[iz2]=="PF") { // lllj events
327              if(fs->dataset_!="fakes" && iz2!=hiPtPfZ2) continue; // shouldn't have this when I extrapolate to signal region
328  
329 <            if(makeFakeTuples && fs->dataset_!="fakes" && ctrl.faketype=="3P1F") {
329 >            if(ctrl.makeFakeTuples && fs->dataset_!="fakes" && ctrl.faketype=="3P1F") {
330                assert(lep3.isLoose || lep4.isLoose);
331                assert(!(lep3.isLoose && lep4.isLoose));
332                double fwgt = lep3.isLoose ? fwgts.fwgt_4 : fwgts.fwgt_3;
# Line 391 | Line 338 | int main(int argc, char** argv)
338              }
339              
340              cutvec["fillPF"] += 1;
341 <            fillAllHists( ctrl, cs, channel, "PF", fs, kine, lep1, lep2, lep3, lep4, wgt);
341 >            fillAllHists( ctrl, cs, channel, "PF", fs, kine, angles, lep1, lep2, lep3, lep4, wgt);
342              if(doJets) fillAllJetHists( ctrl, cs, channel, "PF", fs, kine, lep1, lep2, lep3, lep4, fusion, goodJets, ji, wgt);
343            } else assert(0);
344          }
345        }
346 +      if(ctrl.writessofratio) fs->writeSsofRatios();
347        cout << "\t\t" << setw(7) << nDuplSkipped << " duplicate events skipped" << endl;
348 <      if(makeFakeTuples) {
348 >      if(ctrl.makeFakeTuples) {
349          fs->outFotuple->getFile()->cd();
350          fs->outFotuple->getTree()->Write();
351          fs->outtuple->WriteClose();
352        }
353 +      fs->del();
354      }
355      for(unsigned icut=0; icut<cutstrs.size(); icut++) {
356        cout << setw(22) << cutstrs[icut] << setw(22) << cutvec[cutstrs[icut]] << endl;
# Line 416 | Line 365 | int main(int argc, char** argv)
365    vector<TString> typev;
366    typev.push_back("all");
367  
368 <  bool plotAllCats=false;
420 <  if(plotAllCats) {
368 >  if(ctrl.plotAllCats) {
369      typev.push_back("4e");
370      typev.push_back("4m");
371      typev.push_back("2e2m");
372    }
373  
374 +  ofstream txtOutFile("fakes.txt",ios_base::app);
375 +  txtOutFile << "---> " << ctrl.faketype << " " << plotLabel << endl;
376    // loop over channels
377    for(unsigned itype=0; itype<typev.size(); itype++) {
378      TString type(typev[itype]);
379 <    gSystem->mkdir(ctrl.outdir+"/"+plotLabel+"/"+type,true);
380 <    TFile runHistFile(ctrl.outdir+"/"+plotLabel+"/"+type+"/runs.root","recreate");
379 >    TString fullOutDir(ctrl.outdir+"/"+plotLabel+"/"+type);
380 >    gSystem->mkdir(fullOutDir,true);
381 >    TFile runHistFile(fullOutDir+"/runs.root","recreate");
382      // loop over variables
383      map<TString,TH1D*>::iterator it_v;
384      for(it_v=(*(samplev[0]->hists)[type+"_obs"]).begin(); it_v!=(*(samplev[0]->hists)[type+"_obs"]).end(); it_v++) {
385        TString var((*it_v).first);
386 <      CPlot cplot(var,"",(*(samplev[0]->hists)[type+"_obs"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+plotLabel+"/"+type+"/plots");
386 >      CPlot cplot(var,TString(""),(*(samplev[0]->hists)[type+"_obs"])[var]->GetXaxis()->GetTitle(),TString("events"),TString(fullOutDir+"/plots"));
387        double fakeCounter=0;
388        double ymax=0;
389        // loop over csamples
# Line 446 | Line 397 | int main(int argc, char** argv)
397            scaleHists(hs);
398          
399          if(cs->isdata && !(cs->name=="fakes")) {
400 <          if(ctrl.faketype=="SR" || (ctrl.heavyFlavor && !ctrl.plotWholeSample)) {
400 >          if(ctrl.faketype=="SR" || (ctrl.heavyFlavor && !ctrl.plotWholeSample) || (ctrl.ssof && !ctrl.plotWholeSample)) {
401              cplot.AddHist1D(hs["obs"],cs->legend+": "+integral_str(hs["obs"],0),"E");
402              ymax = max(ymax,hs["obs"]->GetMaximum());
403            }
# Line 455 | Line 406 | int main(int argc, char** argv)
406                cplot.AddHist1D(hs["pred"],cs->legend+": "+integral_str(hs["pred"],0),"E");
407                ymax = max(ymax,hs["pred"]->GetMaximum());
408              } else {
409 <              cplot.AddHist1D(hs["pred"],   "FR predic: "+integral_str(hs["pred"],3),   "hist",kRed);
410 <              cplot.AddHist1D(hs["pred_lo"],  "stat lo: "+integral_str(hs["pred_lo"],2),"hist",kRed,kDashed);
411 <              cplot.AddHist1D(hs["pred_hi"],  "stat hi: "+integral_str(hs["pred_hi"],2),"hist",kRed,kDashed);
409 >              cplot.AddHist1D(hs["pred"],    TString("FR predic: "+integral_str(hs["pred"],3)),   "hist",kRed);
410 >              cplot.AddHist1D(hs["pred_lo"], TString(  "stat lo: "+integral_str(hs["pred_lo"],2)),"hist",kRed,kDashed);
411 >              cplot.AddHist1D(hs["pred_hi"], TString(  "stat hi: "+integral_str(hs["pred_hi"],2)),"hist",kRed,kDashed);
412 >              if(var=="m4l") {
413 >                txtOutFile << "  " << type << ": " << integral_str(hs["pred"],3) << " (" << integral_str(hs["pred_lo"],2) << "," << integral_str(hs["pred_hi"],2) << ")" << endl;
414 >              }
415                ymax = max(ymax,histMax(hs["pred"],hs["pred_lo"],hs["pred_hi"]));
416              }
417            }
418            if(ctrl.faketype=="3P1F") {
419              double ratio = 0;//integrateHist(hs["PF"]) / integrateHist(hs["pred"]);
420 <            cplot.AddHist1D(hs["PF"],cs->legend+": "+integral_str(hs["PF"],0)+", ratio: "+f_to_a(ratio),"E",cs->color);
420 >            cplot.AddHist1D(hs["PF"],TString(cs->legend+": "+integral_str(hs["PF"],0)+", ratio: "+f_to_a(ratio)),"E",cs->color);
421              if(!ctrl.plotWholeSample) {
422 <              cplot.AddToStack(hs["pred"],cs->legend+"(2P2F extrap.): "+integral_str(hs["pred"],1), kRed);
422 >              cplot.AddToStack(hs["pred"],TString(cs->legend+"(2P2F extrap.): "+integral_str(hs["pred"],1)), kRed, -1);
423                ymax = max(ymax,histMax(hs["PF"],cplot.GetStack(),hs["pred"]));
424              } else {
425                ymax = max(ymax,hs["PF"]->GetMaximum());
# Line 482 | Line 436 | int main(int argc, char** argv)
436            if(ctrl.faketype=="3P1F")   mchist = hs["PF"];
437            assert(mchist);
438            bool stackMc = true;
439 <          if(stackMc) cplot.AddToStack(mchist,cs->legend+": "+integral_str(mchist,2),cs->color);
440 <          else        cplot.AddHist1D(mchist,cs->legend +": "+integral_str(mchist,3),"Ehist",cs->color);//cs->color);
439 >          if(stackMc) cplot.AddToStack(mchist,TString(cs->legend+": "+integral_str(mchist,2)),cs->color,-1);
440 >          else        cplot.AddHist1D(mchist,TString(cs->legend +": "+integral_str(mchist,3)),"Ehist",cs->color);//cs->color);
441            if(cs->name=="zz") fakeCounter -= integrateHist(hs["PF"]);
442            ymax = max(ymax,histMax(mchist,cplot.GetStack()));
443          }
# Line 493 | Line 447 | int main(int argc, char** argv)
447          }
448          if(!var.Contains("fusionMVA"))
449            cplot.SetYRange(0,1.2*ymax);
450 <        if(ctrl.faketype=="3P1F") cplot.AddTextBox("diff: "+f_to_a(fakeCounter),.7,.5,.8,.6);
450 >        if(ctrl.faketype=="3P1F") cplot.AddTextBox(TString("diff: "+f_to_a(fakeCounter)),.7,.5,.8,.6);
451        }
452        cplot.Draw(can,true,"png");
453        if(var.Contains("fusionMVA")) {
# Line 503 | Line 457 | int main(int argc, char** argv)
457        }
458      }
459      runHistFile.Close();
460 <    makeHTML(ctrl,type,plotLabel);
460 >    makeHTML(ctrl,type,plotLabel,fullOutDir);
461    }
462 +  txtOutFile.close();
463 +      
464 +  fr2011.FR_struct::~FR_struct();
465 +  fr2012.FR_struct::~FR_struct();
466   }
467   //----------------------------------------------------------------------------------------
468   map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
# Line 542 | Line 500 | map<TString,map<TString,TH1D*>* > init_h
500  
501    for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
502      (*((*it_h).second))["npv"]          = new TH1D(TString("npv") +"_"+(*it_h).first+str,";#bf{N PV};",         50,-0.5,49.5);      (*((*it_h).second))["npv"]->Sumw2();
503 <    (*((*it_h).second))["run"]          = new TH1D(TString("run") +"_"+(*it_h).first+str,";#bf{run};",          50,160800,196535);   (*((*it_h).second))["run"]->Sumw2();
503 >    (*((*it_h).second))["run"]          = new TH1D(TString("run") +"_"+(*it_h).first+str,";#bf{run};",          50,160800,203744/*196535*/);   (*((*it_h).second))["run"]->Sumw2();
504      (*((*it_h).second))["mZ1"]          = new TH1D(TString("mZ1") +"_"+(*it_h).first+str,";#bf{mZ1 [GeV]};",    20,40,120);          (*((*it_h).second))["mZ1"]->Sumw2();
505      (*((*it_h).second))["mZ2"]          = new TH1D(TString("mZ2") +"_"+(*it_h).first+str,";#bf{mZ2 [GeV]};",    30,0,115);           (*((*it_h).second))["mZ2"]->Sumw2();
506      (*((*it_h).second))["mZ2_lo"]       = new TH1D(TString("mZ2_lo") +"_"+(*it_h).first+str,";#bf{mZ2 [GeV]};", 30,0,10);           (*((*it_h).second))["mZ2_lo"]->Sumw2();
507      (*((*it_h).second))["m4l_lo"]       = new TH1D(TString("m4l_lo") +"_"+(*it_h).first+str,";#bf{m4l [GeV]};", 35,70,120);          (*((*it_h).second))["m4l_lo"]->Sumw2();
508      (*((*it_h).second))["m4l_med"]      = new TH1D(TString("m4l_med") +"_"+(*it_h).first+str,";#bf{m4l [GeV]};",35,100,180);          (*((*it_h).second))["m4l_med"]->Sumw2();
509      (*((*it_h).second))["m4l"]          = new TH1D(TString("m4l") +"_"+(*it_h).first+str,";#bf{m4l [GeV]};",    35,100,600);          (*((*it_h).second))["m4l"]->Sumw2();
510 <    (*((*it_h).second))["Z1pt"]         = new TH1D(TString("Z1pt") +"_"+(*it_h).first+str,";#bf{Z1pt [GeV]};",  20,0,200);           (*((*it_h).second))["Z1pt"]->Sumw2();
511 <    (*((*it_h).second))["ZZpt"]         = new TH1D(TString("ZZpt") +"_"+(*it_h).first+str,";#bf{ZZpt [GeV]};",  30,0,200);           (*((*it_h).second))["ZZpt"]->Sumw2();
510 >    // (*((*it_h).second))["Z1pt"]              = new TH1D(TString("Z1pt") +"_"+(*it_h).first+str,";#bf{Z1pt [GeV]};",  20,0,200);           (*((*it_h).second))["Z1pt"]->Sumw2();
511 >    // (*((*it_h).second))["ZZpt"]              = new TH1D(TString("ZZpt") +"_"+(*it_h).first+str,";#bf{ZZpt [GeV]};",  30,0,200);           (*((*it_h).second))["ZZpt"]->Sumw2();
512      (*((*it_h).second))["met"]          = new TH1D(TString("met") +"_"+(*it_h).first+str,";#bf{met [GeV]};",    30,0,100);           (*((*it_h).second))["met"]->Sumw2();
513  
514      (*((*it_h).second))["pt_l1"]        = new TH1D(TString("pt_l1") +"_"+(*it_h).first+str,";#bf{l1 p_{T} [GeV]};",     37,0,120);  (*((*it_h).second))["pt_l1"]->Sumw2();
# Line 626 | Line 584 | map<TString,map<TString,TH1D*>* > init_h
584      (*((*it_h).second))["dR_l3j4"]      = new TH1D(TString("dR_l3j4") +"_"+(*it_h).first+str,";#bf{#Delta R l3,jet 4};",75,0,6);        (*((*it_h).second))["dR_l3j4"]->Sumw2();
585      (*((*it_h).second))["dR_l4j4"]      = new TH1D(TString("dR_l4j4") +"_"+(*it_h).first+str,";#bf{#Delta R l4,jet 4};",75,0,6);        (*((*it_h).second))["dR_l4j4"]->Sumw2();
586  
587 +    int nbins=35;
588 +    (*((*it_h).second))["costheta1"]    = new TH1D(TString("costheta1") +"_"+(*it_h).first+str,";#bf{costheta1};",                nbins,-1,1);          (*((*it_h).second))["costheta1"]->Sumw2();
589 +    (*((*it_h).second))["costheta2"]    = new TH1D(TString("costheta2") +"_"+(*it_h).first+str,";#bf{costheta2};",                nbins,-1,1);          (*((*it_h).second))["costheta2"]->Sumw2();
590 +    (*((*it_h).second))["costhetastar"] = new TH1D(TString("costhetastar") +"_"+(*it_h).first+str,";#bf{costhetastar};",          nbins,-1,1);          (*((*it_h).second))["costhetastar"]->Sumw2();
591 +    (*((*it_h).second))["Phi"]          = new TH1D(TString("Phi") +"_"+(*it_h).first+str,";#bf{Phi};",                            nbins,-3.15,3.15);          (*((*it_h).second))["Phi"]->Sumw2();
592 +    (*((*it_h).second))["Phi1"]         = new TH1D(TString("Phi1") +"_"+(*it_h).first+str,";#bf{Phi1};",                          nbins,-3.15,3.15);          (*((*it_h).second))["Phi1"]->Sumw2();
593 +
594 +    (*((*it_h).second))["pt4l"]         = new TH1D(TString("pt4l") +"_"+(*it_h).first+str,";#bf{pt4l/m4l};",                      nbins-7,0,0.5);          (*((*it_h).second))["pt4l"]->Sumw2();
595 +    (*((*it_h).second))["pt4l_log"]     = new TH1D(TString("pt4l_log") +"_"+(*it_h).first+str,";#bf{pt4l/m4l};",                  nbins,0,1.5);          (*((*it_h).second))["pt4l_log"]->Sumw2();
596 +    (*((*it_h).second))["y4l"]          = new TH1D(TString("y4l") +"_"+(*it_h).first+str,";#bf{y4l};",                            nbins,-2.4,2.4);      (*((*it_h).second))["y4l"]->Sumw2();
597 +    (*((*it_h).second))["Z1pt"]         = new TH1D(TString("Z1pt") +"_"+(*it_h).first+str,";#bf{Z1pt/m4l};",                      nbins,0,0.5);          (*((*it_h).second))["Z1pt"]->Sumw2();
598 +    (*((*it_h).second))["Z1pt_log"]     = new TH1D(TString("Z1pt_log") +"_"+(*it_h).first+str,";#bf{Z1pt/m4l};",                  nbins,0,1.5);          (*((*it_h).second))["Z1pt_log"]->Sumw2();
599 +    (*((*it_h).second))["Z2pt"]         = new TH1D(TString("Z2pt") +"_"+(*it_h).first+str,";#bf{Z2pt/m4l};",                      nbins,0,0.5);          (*((*it_h).second))["Z2pt"]->Sumw2();
600 +    (*((*it_h).second))["Z2pt_log"]     = new TH1D(TString("Z2pt_log") +"_"+(*it_h).first+str,";#bf{Z2pt/m4l};",                  nbins,0,1.5);          (*((*it_h).second))["Z2pt_log"]->Sumw2();
601 +
602 +    (*((*it_h).second))["ZZdotZ1"]      = new TH1D(TString("ZZdotZ1") +"_"+(*it_h).first+str,";#bf{ZZ #bullet Z1/(m4l*mZ1)};",    nbins,-.2,2.5);          (*((*it_h).second))["ZZdotZ1"]->Sumw2();
603 +    (*((*it_h).second))["ZZdotZ2"]      = new TH1D(TString("ZZdotZ2") +"_"+(*it_h).first+str,";#bf{ZZ #bullet Z2/(m4l*mZ2)};",    nbins,-1,3);          (*((*it_h).second))["ZZdotZ2"]->Sumw2();
604 +    (*((*it_h).second))["ZZdotZ1_log"]  = new TH1D(TString("ZZdotZ1_log") +"_"+(*it_h).first+str,";#bf{ZZ #bullet Z1/(m4l*mZ1)};",nbins,-1.5,15);          (*((*it_h).second))["ZZdotZ1_log"]->Sumw2();
605 +    (*((*it_h).second))["ZZdotZ2_log"]  = new TH1D(TString("ZZdotZ2_log") +"_"+(*it_h).first+str,";#bf{ZZ #bullet Z2/(m4l*mZ2)};",nbins,-1.5,15);          (*((*it_h).second))["ZZdotZ2_log"]->Sumw2();
606 +    (*((*it_h).second))["dphi1"]        = new TH1D(TString("dphi1") +"_"+(*it_h).first+str,";#bf{dphi1};",                        nbins,-1,1);          (*((*it_h).second))["dphi1"]->Sumw2();
607 +    (*((*it_h).second))["dphi2"]        = new TH1D(TString("dphi2") +"_"+(*it_h).first+str,";#bf{dphi2};",                        nbins,-1,1);          (*((*it_h).second))["dphi2"]->Sumw2();
608 +    (*((*it_h).second))["dphi1_log"]    = new TH1D(TString("dphi1_log") +"_"+(*it_h).first+str,";#bf{dphi1};",                    nbins,-1,1);          (*((*it_h).second))["dphi1_log"]->Sumw2();
609 +    (*((*it_h).second))["dphi2_log"]    = new TH1D(TString("dphi2_log") +"_"+(*it_h).first+str,";#bf{dphi2};",                    nbins,-1,1);          (*((*it_h).second))["dphi2_log"]->Sumw2();
610 +
611      (*((*it_h).second))["fusionMVA"]            = new TH1D(TString("fusionMVA")+"_"+(*it_h).first+str,";#bf{MVA output};",        100,-1,1);      (*((*it_h).second))["fusionMVA"]->Sumw2();
612      (*((*it_h).second))["fusionMVA_lo"]         = new TH1D(TString("fusionMVA_lo")+"_"+(*it_h).first+str,";#bf{MVA output};",        100,-1,-.83);      (*((*it_h).second))["fusionMVA_lo"]->Sumw2();
613      (*((*it_h).second))["fusionMVA_med"]                = new TH1D(TString("fusionMVA_med")+"_"+(*it_h).first+str,";#bf{MVA output};",        100,-.5,1);      (*((*it_h).second))["fusionMVA_med"]->Sumw2();
# Line 636 | Line 618 | map<TString,map<TString,TH1D*>* > init_h
618    return hists;
619   }
620   //--------------------------------------------------------------------------------------------------
621 < void makeHTML(FOFlags &ctrl, TString type, TString plotLabel)
621 > void makeHTML(FOFlags &ctrl, TString type, TString plotLabel, TString fullOutDir)
622   {
623    TString tmpString(plotLabel);
624    tmpString.ReplaceAll("/"," ");
625    TString title(ctrl.faketype+", "+tmpString+", "+type);
626 <  TString htmlfname(ctrl.outdir+"/"+plotLabel+"/"+type+"/plots.html");
626 >  TString htmlfname(fullOutDir+"/plots.html");
627    ofstream htmlfile(htmlfname);
628  
629    htmlfile << "<!DOCTYPE html" << endl;
# Line 695 | Line 677 | void makeHTML(FOFlags &ctrl, TString typ
677    htmlfile << "</html>" << endl;
678    htmlfile.close();
679  
680 +  TString anglefname(htmlfname);
681 +  anglefname.ReplaceAll("plots.html","angleplots.html");
682 +  htmlfile.open(anglefname);
683 +
684 +  htmlfile << "<!DOCTYPE html" << endl;
685 +  htmlfile << "    PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
686 +  htmlfile << "<html>" << endl;
687 +
688 +  htmlfile << "<head><title>"+title+"</title></head>" << endl;
689 +  htmlfile << "<body bgcolor=\"000000\">" << endl;
690 +  htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
691 +
692 +  htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">Angles & Co.</h3>" << endl;
693 +  htmlfile << "<hr />" << endl;
694 +  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
695 +  htmlfile << "</table>" << endl;
696 +
697 +  htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">Angles:</h3>" << endl;
698 +  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
699 +
700 +  htmlfile << "<tr>" << endl;
701 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/costheta1.png\"><img src=\"plots/costheta1.png\" alt=\"plots/costheta1.png\" width=\"100%\"></a></td>" << endl;
702 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/costheta2.png\"><img src=\"plots/costheta2.png\" alt=\"plots/costheta2.png\" width=\"100%\"></a></td>" << endl;
703 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/costhetastar.png\"><img src=\"plots/costhetastar.png\" alt=\"plots/costhetastar.png\" width=\"100%\"></a></td>" << endl;
704 +  htmlfile << "<td width=\"25%\"></td>" << endl;
705 +  htmlfile << "</tr>" << endl;
706 +
707 +  htmlfile << "<tr>" << endl;
708 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Phi.png\"><img src=\"plots/Phi.png\" alt=\"plots/Phi.png\" width=\"100%\"></a></td>" << endl;
709 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Phi1.png\"><img src=\"plots/Phi1.png\" alt=\"plots/Phi1.png\" width=\"100%\"></a></td>" << endl;
710 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
711 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
712 +  htmlfile << "</tr>" << endl;
713 +
714 +  htmlfile << "</table>" << endl;
715 +  htmlfile << "<hr />" << endl;
716 +  htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">Boson pt variables: " << endl;
717 +  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
718 +
719 +  htmlfile << "<tr>" << endl;
720 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/y4l.png\"><img src=\"plots/y4l.png\" alt=\"plots/y4l.png\" width=\"100%\"></a></td>" << endl;
721 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/pt4l.png\"><img src=\"plots/pt4l.png\" alt=\"plots/pt4l.png\" width=\"100%\"></a></td>" << endl;
722 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/pt4l_log.png\"><img src=\"plots/pt4l_log.png\" alt=\"plots/pt4l_log.png\" width=\"100%\"></a></td>" << endl;
723 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
724 +  htmlfile << "<tr>" << endl;
725 +
726 +  htmlfile << "</tr>" << endl;
727 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Z1pt.png\"><img src=\"plots/Z1pt.png\" alt=\"plots/Z1pt.png\" width=\"100%\"></a></td>" << endl;
728 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Z2pt.png\"><img src=\"plots/Z2pt.png\" alt=\"plots/Z2pt.png\" width=\"100%\"></a></td>" << endl;
729 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Z1pt_log.png\"><img src=\"plots/Z1pt_log.png\" alt=\"plots/Z1pt_log.png\" width=\"100%\"></a></td>" << endl;
730 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Z2pt_log.png\"><img src=\"plots/Z2pt_log.png\" alt=\"plots/Z2pt_log.png\" width=\"100%\"></a></td>" << endl;
731 +  htmlfile << "</tr>" << endl;
732 +
733 +  htmlfile << "<tr>" << endl;
734 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ZZdotZ1.png\"><img src=\"plots/ZZdotZ1.png\" alt=\"plots/ZZdotZ1.png\" width=\"100%\"></a></td>" << endl;
735 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ZZdotZ2.png\"><img src=\"plots/ZZdotZ2.png\" alt=\"plots/ZZdotZ2.png\" width=\"100%\"></a></td>" << endl;
736 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ZZdotZ1_log.png\"><img src=\"plots/ZZdotZ1_log.png\" alt=\"plots/ZZdotZ1_log.png\" width=\"100%\"></a></td>" << endl;
737 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ZZdotZ2_log.png\"><img src=\"plots/ZZdotZ2_log.png\" alt=\"plots/ZZdotZ2_log.png\" width=\"100%\"></a></td>" << endl;
738 +  htmlfile << "</tr>" << endl;
739 +
740 +  htmlfile << "<tr>" << endl;
741 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dphi1.png\"><img src=\"plots/dphi1.png\" alt=\"plots/dphi1.png\" width=\"100%\"></a></td>" << endl;
742 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dphi2.png\"><img src=\"plots/dphi2.png\" alt=\"plots/dphi2.png\" width=\"100%\"></a></td>" << endl;
743 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dphi1_log.png\"><img src=\"plots/dphi1_log.png\" alt=\"plots/dphi1_log.png\" width=\"100%\"></a></td>" << endl;
744 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dphi2_log.png\"><img src=\"plots/dphi2_log.png\" alt=\"plots/dphi2_log.png\" width=\"100%\"></a></td>" << endl;
745 +  htmlfile << "</tr>" << endl;
746 +
747 +
748 +  htmlfile << "</table>" << endl;
749 +
750 +  htmlfile << "</body>" << endl;
751 +  htmlfile << "</html>" << endl;
752 +  htmlfile.close();
753 +
754    TString leptonfname(htmlfname);
755    leptonfname.ReplaceAll("plots.html","leptonplots.html");
756    htmlfile.open(leptonfname);
701  cout << "opening: " << leptonfname << endl;
757  
758    htmlfile << "<!DOCTYPE html" << endl;
759    htmlfile << "    PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
# Line 724 | Line 779 | void makeHTML(FOFlags &ctrl, TString typ
779    htmlfile << "</tr>" << endl;
780  
781    htmlfile << "<tr>" << endl;
782 <  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dz_l1.png\"><img src=\"plots/dz_l1.png\" alt=\"plots/dz_l1.png\" width=\"100%\"></a></td>" << endl;
783 <  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dz_l2.png\"><img src=\"plots/dz_l2.png\" alt=\"plots/dz_l2.png\" width=\"100%\"></a></td>" << endl;
782 >  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l1.png\"><img src=\"plots/ip3ds_l1.png\" alt=\"plots/ip3ds_l1.png\" width=\"100%\"></a></td>" << endl;
783 >  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l2.png\"><img src=\"plots/ip3ds_l2.png\" alt=\"plots/ip3ds_l2.png\" width=\"100%\"></a></td>" << endl;
784    htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
785    htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
786    htmlfile << "</tr>" << endl;
787  
788    htmlfile << "<tr>" << endl;
734  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l1.png\"><img src=\"plots/ip3ds_l1.png\" alt=\"plots/ip3ds_l1.png\" width=\"100%\"></a></td>" << endl;
735  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l2.png\"><img src=\"plots/ip3ds_l2.png\" alt=\"plots/ip3ds_l2.png\" width=\"100%\"></a></td>" << endl;
789    htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/d0_l1.png\"><img src=\"plots/d0_l1.png\" alt=\"plots/d0_l1.png\" width=\"100%\"></a></td>" << endl;
790    htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/d0_l2.png\"><img src=\"plots/d0_l2.png\" alt=\"plots/d0_l2.png\" width=\"100%\"></a></td>" << endl;
791 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
792 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
793 +  htmlfile << "</tr>" << endl;
794 +
795 +  htmlfile << "<tr>" << endl;
796 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dz_l1.png\"><img src=\"plots/dz_l1.png\" alt=\"plots/dz_l1.png\" width=\"100%\"></a></td>" << endl;
797 +  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/dz_l2.png\"><img src=\"plots/dz_l2.png\" alt=\"plots/dz_l2.png\" width=\"100%\"></a></td>" << endl;
798 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
799 +  htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
800    htmlfile << "</tr>" << endl;
801  
802    // htmlfile << "<tr>" << endl;
# Line 896 | Line 958 | void makeHTML(FOFlags &ctrl, TString typ
958   //----------------------------------------------------------------------------------------
959   void fillHist(CSample *cs, TString channel, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
960   {
961 +  assert(wgt==wgt);
962 +  assert(wgt_lo==wgt_lo);
963 +  assert(wgt_hi==wgt_hi);
964    (*(cs->hists)["all_"+type])[var]->Fill(           val,   wgt);
965    if(channel!="")
966       (*(cs->hists)[channel+"_"+type])[var]->Fill(   val,   wgt);
# Line 909 | Line 974 | void fillHist(CSample *cs, TString chann
974    }
975   }
976   //----------------------------------------------------------------------------------------
977 < void fillAllHists(FOFlags ctrl, CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine,
977 > void fillAllHists(FOFlags ctrl, CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine, Angles angles,
978                    SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4,
979                    double wgt, double wgt_lo, double wgt_hi)
980   {
# Line 921 | Line 986 | void fillAllHists(FOFlags ctrl, CSample
986    fillHist( cs, channel, type, "m4l_lo"      , kine.m4l       , wgt, wgt_lo, wgt_hi);                
987    fillHist( cs, channel, type, "m4l_med"     , kine.m4l       , wgt, wgt_lo, wgt_hi);                
988    fillHist( cs, channel, type, "m4l"         , kine.m4l       , wgt, wgt_lo, wgt_hi);                
924  fillHist( cs, channel, type, "Z1pt"        , kine.Z1pt      , wgt, wgt_lo, wgt_hi);                
925  fillHist( cs, channel, type, "ZZpt"        , kine.ZZpt      , wgt, wgt_lo, wgt_hi);                
989    fillHist( cs, channel, type, "met"         , fs->info->met  , wgt, wgt_lo, wgt_hi);                
990    fillHist( cs, channel, type, "ip3ds_l1"    , lep1.ip3dSig   , wgt, wgt_lo, wgt_hi);
991    fillHist( cs, channel, type, "ip3ds_l2"    , lep2.ip3dSig   , wgt, wgt_lo, wgt_hi);
# Line 950 | Line 1013 | void fillAllHists(FOFlags ctrl, CSample
1013    fillHist( cs, channel, type, "eta_l4"      , lep4.vec.Eta() , wgt, wgt_lo, wgt_hi);
1014    fillHist( cs, channel, type, "dR_l3l4_lo"  , dr(lep3,lep4)  , wgt, wgt_lo, wgt_hi);
1015    fillHist( cs, channel, type, "dR_l3l4"     , dr(lep3,lep4)  , wgt, wgt_lo, wgt_hi);
1016 +
1017 +  fillHist( cs, channel, type, "costheta1"      , angles.costheta1                  ,   wgt, wgt_lo, wgt_hi);
1018 +  fillHist( cs, channel, type, "costheta2"      , angles.costheta2                  ,   wgt, wgt_lo, wgt_hi);
1019 +  fillHist( cs, channel, type, "costhetastar"   , angles.costhetastar               ,   wgt, wgt_lo, wgt_hi);
1020 +  fillHist( cs, channel, type, "Phi"            , angles.Phi                        ,   wgt, wgt_lo, wgt_hi);
1021 +  fillHist( cs, channel, type, "Phi1"           , angles.Phi1                       ,   wgt, wgt_lo, wgt_hi);
1022 +  fillHist( cs, channel, type, "pt4l"           , kine.ZZpt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1023 +  fillHist( cs, channel, type, "pt4l_log"       , kine.ZZpt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1024 +  fillHist( cs, channel, type, "y4l"            , kine.ZZy                          ,   wgt, wgt_lo, wgt_hi);
1025 +  fillHist( cs, channel, type, "Z1pt"           , kine.Z1pt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1026 +  fillHist( cs, channel, type, "Z1pt_log"       , kine.Z1pt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1027 +  fillHist( cs, channel, type, "Z2pt"           , kine.Z2pt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1028 +  fillHist( cs, channel, type, "Z2pt_log"       , kine.Z2pt/kine.m4l                ,   wgt, wgt_lo, wgt_hi);
1029 +  fillHist( cs, channel, type, "ZZdotZ1"        , kine.ZZdotZ1/(kine.m4l*kine.mZ1)  ,   wgt, wgt_lo, wgt_hi);
1030 +  fillHist( cs, channel, type, "ZZdotZ2"        , kine.ZZdotZ2/(kine.m4l*kine.mZ2)  ,   wgt, wgt_lo, wgt_hi);
1031 +  fillHist( cs, channel, type, "ZZdotZ1_log"    , kine.ZZdotZ1/(kine.m4l*kine.mZ1)  ,   wgt, wgt_lo, wgt_hi);
1032 +  fillHist( cs, channel, type, "ZZdotZ2_log"    , kine.ZZdotZ2/(kine.m4l*kine.mZ2)  ,   wgt, wgt_lo, wgt_hi);
1033 +  fillHist( cs, channel, type, "dphi1"          , kine.ZZptZ1ptCosDphi              ,   wgt, wgt_lo, wgt_hi);
1034 +  fillHist( cs, channel, type, "dphi2"          , kine.ZZptZ2ptCosDphi              ,   wgt, wgt_lo, wgt_hi);
1035 +  fillHist( cs, channel, type, "dphi1_log"      , kine.ZZptZ1ptCosDphi              ,   wgt, wgt_lo, wgt_hi);
1036 +  fillHist( cs, channel, type, "dphi2_log"      , kine.ZZptZ2ptCosDphi              ,   wgt, wgt_lo, wgt_hi);
1037 +
1038   }
1039   //----------------------------------------------------------------------------------------
1040   void fillAllJetHists( FOFlags ctrl, CSample *cs, TString channel, TString type, filestuff *fs, KinematicsStruct kine,
# Line 959 | Line 1044 | void fillAllJetHists( FOFlags ctrl, CSam
1044   {
1045    double fMVAval = (ctrl.uncert=="") ? 0 : fusion.reader->EvaluateMVA("BDTG");
1046    dr_struct drs = fill_dr_struct(goodJets,lep1,lep2,lep3,lep4);
962  drs.check();
1047    fillHist( cs, channel, type, "pt_j1"  ,     ji.ptJet1,        wgt, wgt_lo, wgt_hi);
1048    fillHist( cs, channel, type, "pt_j2"  ,     ji.ptJet2,        wgt, wgt_lo, wgt_hi);
1049    fillHist( cs, channel, type, "pt_j3"  ,     ji.ptJet3,        wgt, wgt_lo, wgt_hi);
# Line 1033 | Line 1117 | bool passHlt(FOFlags &ctrl, TrigInfo ti,
1117    return pass;
1118   }
1119   //----------------------------------------------------------------------------------------
1036 void fillLeptonFakeWeights(lepFrs &fwgts, FR_struct *fr, SimpleLepton lep3, SimpleLepton lep4)
1037 {
1038  fwgts.fwgt_3    = get_fake_weight("",lep3,*fr);
1039  fwgts.fwgt_lo_3 = get_fake_weight("lo",lep3,*fr);
1040  fwgts.fwgt_hi_3 = get_fake_weight("hi",lep3,*fr);
1041  fwgts.fwgt_4    = get_fake_weight("",lep4,*fr);
1042  fwgts.fwgt_lo_4 = get_fake_weight("lo",lep4,*fr);
1043  fwgts.fwgt_hi_4 = get_fake_weight("hi",lep4,*fr);
1044 }
1045 //----------------------------------------------------------------------------------------
1046 dr_struct fill_dr_struct(vector<SimpleLepton> jets, SimpleLepton lep1, SimpleLepton lep2, SimpleLepton lep3, SimpleLepton lep4)
1047 {
1048  dr_struct drs;
1049  if(jets.size() > 0) {
1050    drs.l1j1 = dr(lep1,jets[0]);
1051    drs.l2j1 = dr(lep2,jets[0]);
1052    drs.l3j1 = dr(lep3,jets[0]);
1053    drs.l4j1 = dr(lep4,jets[0]);
1054  }
1055  if(jets.size() > 1) {
1056    drs.l1j2 = dr(lep1,jets[1]);
1057    drs.l2j2 = dr(lep2,jets[1]);
1058    drs.l3j2 = dr(lep3,jets[1]);
1059    drs.l4j2 = dr(lep4,jets[1]);
1060  }
1061  if(jets.size() > 2) {
1062    drs.l1j3 = dr(lep1,jets[2]);
1063    drs.l2j3 = dr(lep2,jets[2]);
1064    drs.l3j3 = dr(lep3,jets[2]);
1065    drs.l4j3 = dr(lep4,jets[2]);
1066  }
1067  if(jets.size() > 3) {
1068    drs.l1j4 = dr(lep1,jets[3]);
1069    drs.l2j4 = dr(lep2,jets[3]);
1070    drs.l3j4 = dr(lep3,jets[3]);
1071    drs.l4j4 = dr(lep4,jets[3]);
1072  }
1073  return drs;
1074 }
1075 //----------------------------------------------------------------------------------------
1120   void init_cuts( vector<TString> &cutstrs, map<TString,int> &cutvec)
1121   {
1122    cutstrs.push_back("start"  );    cutvec["start"]  = 0;
# Line 1125 | Line 1169 | map<TString,TH1D*> setHistSet(CSample *c
1169    return hset;
1170   }
1171   //----------------------------------------------------------------------------------------
1128 void shiftOverflows(map<TString,TH1D*> &hset)
1129 {
1130  for(map<TString,TH1D*>::iterator it=hset.begin(); it!=hset.end(); it++)
1131    shiftOverflows((*it).second);
1132 }
1133 //----------------------------------------------------------------------------------------
1134 void scaleHists(map<TString,TH1D*> &hset)
1135 {
1136  for(map<TString,TH1D*>::iterator it=hset.begin(); it!=hset.end(); it++)
1137    (*it).second->Scale(1./integrateHist((*it).second));
1138 }
1139 //----------------------------------------------------------------------------------------
1172   void fillFakeTuple(CSample *cs, filestuff *fs, EventData evtdata, unsigned ientry, double fillweight)
1173   {
1174    fs->getentry(ientry,"","zznt"); // make sure we've got this entry for all the branches
1175 +  fillAngles(evtdata,*fs->angles);
1176    fillKinematics(evtdata,*fs->kine); // reminder: kine has the fake information, while fs->kine *had*, until this line, whatever ZPlusX put in it
1177    fs->weights->w = fillweight;
1178    fs->outtuple->Fill();
1179    fs->outFotuple->Fill();
1180   }
1181 + //----------------------------------------------------------------------------------------
1182 + void incrementSsofCounters(FOFlags &ctrl, filestuff *fs, double minMz2,
1183 +                           pair<int,int> *best_z_indices,
1184 +                           double wgt)
1185 + {
1186 +  assert(ctrl.ssof); // only makes sense to write this out when we're thinking about ssof method
1187 +  vector<pair<SimpleLepton,SimpleLepton> > z2cands;
1188 +  vector<TString> types;
1189 +  
1190 +  int hiPtFailZ2 = findZ2CandidatesPassFail(fs->failingL, z2cands, types, "XX-SF", minMz2, "fail");
1191 +  if(hiPtFailZ2 == -1) return; // no z2 found
1192 +  SimpleLepton lep1 = (*fs->passingL)[best_z_indices->first];
1193 +  SimpleLepton lep2 = (*fs->passingL)[best_z_indices->second];
1194 +  SimpleLepton lep3 = z2cands[hiPtFailZ2].first;
1195 +  SimpleLepton lep4 = z2cands[hiPtFailZ2].second;
1196 +  TString channel = getChannel(lep1,lep2,lep3,lep4);
1197 +  TString sign = (lep3.charge == lep4.charge) ? "SS" : "OS";
1198  
1199 <          // double mjj=0,dEta=0,etaProd=0,nCentralJets=0;
1200 <          // if(jet1 && jet2) {
1201 <          //   cutvec["twoJets"] += 1;
1202 <          //   if(ctrl.debug) cout << "looping through " << goodJets.size() << " good jets" << endl;
1203 <          //   if(ctrl.debug) {cout << "jet 1 (" << jet1 << "): "; jet1->print();}
1204 <          //   if(ctrl.debug) {cout << "jet 2 (" << jet2 << "): "; jet2->print();}
1205 <          //   for(unsigned ijet=0; ijet<goodJets.size(); ijet++) {
1206 <          //     SimpleLepton *jet = goodJets[ijet];
1207 <          //     if(jet == jet1) {
1208 <          //    if(ctrl.debug) cout << "  " << jet << " this is jet 1" << endl;
1209 <          //    continue;
1210 <          //     }
1211 <          //     if(jet == jet2) {
1212 <          //    if(ctrl.debug) cout << "  " << jet << " this is jet 2" << endl;
1213 <          //    continue;
1214 <          //     }
1215 <          //     double eta = jet->vec.Eta();
1216 <          //     double eta1 = jet1->vec.Eta();
1217 <          //     double eta2 = jet2->vec.Eta();
1218 <          //     if(eta1 > eta2) {
1219 <          //    if(eta > eta2 && eta < eta1)
1220 <          //      nCentralJets++;
1221 <          //     } else if(eta2 > eta1) {
1222 <          //    if(eta > eta1 && eta < eta2)
1223 <          //      nCentralJets++;
1224 <          //     }
1225 <          //     if(ctrl.debug) cout << "  nCentralJets: " << nCentralJets << endl;
1226 <          //   }
1227 <          //   double mjjMin = 400;
1228 <          //   double dEtaMin = 4;
1229 <          //   TLorentzVector dijet(jet1->vec + jet2->vec);
1230 <          //   mjj = dijet.M();
1231 <          //   if(ctrl.debug) cout << "setting mjj: " << mjj << "(using " << jet1->vec.Pt() << " " << jet2->vec.Pt() << endl;
1232 <          //   dEta = jet1->vec.Eta() - jet2->vec.Eta();
1233 <          //   etaProd = jet1->vec.Eta()*jet2->vec.Eta();
1234 <          //   bool isVbf = (mjj > mjjMin) && (fabs(dEta) > dEtaMin) && (etaProd < 0) && (nCentralJets==0);
1235 <          // }
1236 <
1199 >  fs->counters_[sign+"_"+channel] += wgt;
1200 >  fs->counters_[sign] += wgt;
1201 > }
1202 > //----------------------------------------------------------------------------------------
1203 > double getSsofWgt(TString channel, map<TString,double> &ssofRatios)
1204 > {
1205 >  double returnval = ssofRatios[channel];
1206 >  if(returnval < 0.1 || returnval > 3) {
1207 >    cout << "\nERROR! ssof ratio out of bounds for " << channel << ": " << returnval << endl;
1208 >    assert(0);
1209 >  }
1210 >  return returnval;
1211 > }
1212 > //----------------------------------------------------------------------------------------
1213 > map<TString,double> initSsofRatios(vector<CSample*> &samplev)
1214 > {
1215 >  double OS,SS;
1216 >  double OS_4e,SS_4e;
1217 >  double OS_4m,SS_4m;
1218 >  double OS_2e2m,SS_2e2m;
1219 >  cout << "Initializing ssof ratios: " << endl;
1220 >  for(unsigned isam=0; isam<samplev.size(); isam++) {
1221 >    CSample *cs = samplev[isam];
1222 >    for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
1223 >      filestuff *fs = (cs->fsv)[ifs];
1224 >      if(fs->useSsofRatio_) {
1225 >        cout << "  adding from: " << fs->fname_ << endl;
1226 >        cout << "        " << fs->counters_["ratio"] << setw(8) << fixed << setprecision(2) << fs->counters_["ratio_4e"] << setw(8) << fixed << setprecision(2) << fs->counters_["ratio_4m"] << setw(8) << fixed << setprecision(2) << fs->counters_["ratio_2e2m"] << endl;
1227 >        OS      += fs->counters_["OS"];
1228 >        SS      += fs->counters_["SS"];
1229 >        OS_4e   += fs->counters_["OS_4e"];
1230 >        SS_4e   += fs->counters_["SS_4e"];
1231 >        OS_4m   += fs->counters_["OS_4m"];
1232 >        SS_4m   += fs->counters_["SS_4m"];
1233 >        OS_2e2m += fs->counters_["OS_2e2m"];
1234 >        SS_2e2m += fs->counters_["SS_2e2m"];
1235 >      }
1236 >    }
1237 >  }
1238 >  assert(OS!=0 && SS!=0 && OS_4e!=0 && SS_4e!=0 && OS_4m!=0 && SS_4m!=0 && OS_2e2m && SS_2e2m!=0);
1239 >  map<TString,double> ssofRatios;
1240 >  ssofRatios[""]     = OS/SS;
1241 >  ssofRatios["4e"]   = OS_4e/SS_4e;
1242 >  ssofRatios["4m"]   = OS_4m/SS_4m;
1243 >  ssofRatios["2e2m"] = OS_2e2m/SS_2e2m;
1244 >  cout << "returning: " << ssofRatios[""] << setw(8) << fixed << setprecision(2) << ssofRatios["4e"] << setw(8) << fixed << setprecision(2) << ssofRatios["4m"] << setw(8) << fixed << setprecision(2) << ssofRatios["2e2m"] << endl;
1245 >  return ssofRatios;
1246 > }
1247 >        

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines