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

Comparing UserCode/MitHzz4l/Selection/src/applySelection.cc (file contents):
Revision 1.33 by khahn, Sun Jun 3 15:59:56 2012 UTC vs.
Revision 1.56 by anlevin, Fri Oct 12 12:26:12 2012 UTC

# Line 33 | Line 33 | using namespace std;
33   #include "Electron.h"
34   #include "Muon.h"
35   #include "Vertex.h"
36 + #include "PFMet.h"
37   #include "PFCandidate.h"
38   #include "PFCandidateCol.h"
39   #include "PileupInfoCol.h"
# Line 40 | Line 41 | using namespace std;
41   #include "MCParticle.h"
42   #include "TriggerMask.h"
43   #include "TriggerTable.h"
44 + #include "TriggerObject.h"
45 + #include "TriggerObjectRel.h"
46   #include "Names.h"
47   #include "BaseMod.h"
48 + #include "TriggerObjectsTable.h"
49  
50   //
51   // our headers
# Line 49 | Line 53 | using namespace std;
53   #include "ParseArgs.h"
54   #include "MuonSelection.h"
55   #include "ElectronSelection.h"
56 + #include "ElectronMomentumCorrection.h"
57   #include "IsolationSelection.h"
58   #include "ReferenceSelection.h"
59  
# Line 65 | Line 70 | using namespace std;
70   #include "RunLumiRangeMap.h"
71   #include "SampleWeight.h"
72   #include "EfficiencyWeightsInterface.h"
73 + #include "SelectionFuncs.h"
74  
75   #include "SimpleLepton.h"
76  
# Line 77 | Line 83 | using namespace mithep;
83   //
84   // globals
85   //----------------------------------------------------------------------------
86 < TH1D * hpu;
86 > TH1D *hpu_2011, *hpu_2012;
87   RunLumiRangeMap rlrm;
88   vector<SimpleLepton> failingLeptons ; // for fake estimation
89   vector<SimpleLepton> passingLeptons;      // for fake estimation
# Line 86 | Line 92 | vector<vector<unsigned> > zcutvec;
92   vector<vector<unsigned> > zzcutvec;
93   map<unsigned,float>       evtrhoMap;
94   vector<string> cutstrs;
95 < bool passes_HLT_MC;
95 > bool passes_HLT_MC,passes_HLT;
96   vector<bool>   PFnoPUflag;;
97  
98   //
99   // prototypes
100   //----------------------------------------------------------------------------
101 < void initPUWeights();
96 < double getPUWeight(unsigned npu);
97 < void setEra(string, ControlFlags&);
98 < void setEffiencyWeights(EventData&, WeightStruct& );
99 < void initRunLumiRangeMap();
100 < void initEvtRhoMap(map<unsigned,float> &);
101 < unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
102 <                         vector<bool> &pfNoPileUpFlag,
103 <                         const mithep::Array<mithep::Vertex>     * vtxArr );
101 > void initRunLumiRangeMap(ControlFlags &ctrl);
102   //----------------------------------------------------------------------------
103  
106
104   //
105   // MAIN
106   //----------------------------------------------------------------------------
# Line 169 | Line 166 | int main(int argc, char** argv)
166    //--------------------------------------------------------------------------------------------------------------
167    TChain * chain = new TChain("Events");
168    TChain * hltchain = new TChain("HLT");
169 +  TChain * runchain = new TChain(Names::gkRunTreeName);
170  
171    string fname;
172 +  unsigned total_unskimmed=0;
173    if( !ctrl.inputfiles.empty() ) {
174      ifstream f(ctrl.inputfiles.c_str());
175      while (f >> fname) {
# Line 178 | Line 177 | int main(int argc, char** argv)
177        cout << "adding inputfile : " << fname.c_str() << endl;
178        entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
179        cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
180 +      total_unskimmed += entrymap[string(fname.c_str())];
181        chain->AddFile(fname.c_str());
182        hltchain->AddFile(fname.c_str());
183 +      runchain->AddFile(fname.c_str());
184      }
185    } else {
186      cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
187 <    unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
188 <    cout << "tmpent: " << tmpent << endl;
189 <    entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
190 <    cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
187 >    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
188 >    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
189 >    cout << "unskimmed entries: " << unsk_ents << endl;
190 >    total_unskimmed += unsk_ents;
191      chain->AddFile(ctrl.inputfile.c_str());
192      hltchain->AddFile(ctrl.inputfile.c_str());
193 +    runchain->AddFile(ctrl.inputfile.c_str());
194    }
195 +  // // write the total number of unskimmed events that went into making this output file to a text file
196 +  writeEntries(ctrl,total_unskimmed);
197  
198    const char * ofname;
199    if( strcmp( ctrl.outputfile.c_str(), "") ) {
# Line 203 | Line 207 | int main(int argc, char** argv)
207    SimpleTable xstab(xspath.c_str());
208  
209  
210 +
211    //
212    // Setup
213    //--------------------------------------------------------------------------------------------------------------
# Line 219 | Line 224 | int main(int argc, char** argv)
224  
225    
226    FOTuple foTree( nt.getFile(), (const char*)"FOtree");
227 +  foTree.makeInfoBranch(evtinfo);
228    foTree.makeFailedLeptonBranch(failingLeptons);
229    foTree.makePassedLeptonBranch(passingLeptons);
230  
231    TH1F * h_wf11_hpt;
232    if(ctrl.mc) {
233      nt.makeWeightBranch(weights);
234 <    //    nt.makeGenInfoBranch(geninfo);
235 <    initEfficiencyWeights();
234 >    if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
235 >    //initEfficiencyWeights();
236      initPUWeights();
237  
238      // Higgs only, pt reweighting
# Line 240 | Line 246 | int main(int argc, char** argv)
246      }
247  
248    } else {
249 <    initRunLumiRangeMap();
249 >    initRunLumiRangeMap(ctrl);
250    }
251  
252    //  initMuonIDMVA();
253    initElectronIDMVAV1();
254    initTrigger();
255 +  ElectronMomentumCorrection electron_mom_corr;
256 +  electron_mom_corr.initialize_regression("../EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
257 +  TrigInfo ti;
258 +  initAnalysisTriggers(ti);
259  
250  // hack
251  initEvtRhoMap(evtrhoMap);
252
253  
260    //
261    // Setup tree I/O
262    //--------------------------------------------------------------------------------------------------------------
# Line 258 | Line 264 | int main(int argc, char** argv)
264    TTree *eventTree=0;  
265    string currentFile("");
266  
267 <  mithep::EventHeader *info    = new mithep::EventHeader();
267 >  UInt_t fNMaxTriggers = 1024;
268 >  mithep::EventHeader                         *info          = new mithep::EventHeader();
269 >  mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
270    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
271    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
272    mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
# Line 269 | Line 277 | int main(int argc, char** argv)
277    mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
278    mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
279    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
280 <  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
280 >  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable(fNMaxTriggers);
281    vector<string>                              *hltTableStrings  = new vector<string>();
282 +  vector<string>                              *hltLabelStrings  = new vector<string>();
283 +  mithep::Array<mithep::TriggerObject>        *hltObjArr     = new mithep::Array<mithep::TriggerObject>();
284 +  mithep::Array<mithep::TriggerObjectRel>     *hltRelsArr    = new mithep::Array<mithep::TriggerObjectRel>();
285 +  std::vector<std::string>                    *hltTab        = new vector<string>();
286 +  std::vector<std::string>                    *hltLab        = new vector<string>();
287 +  TriggerObjectsTable                         *fTrigObjs     = new TriggerObjectsTable(hltTable,fNMaxTriggers);
288 +
289 +  RunInfo                                     *runInfo       = new RunInfo();
290  
291    TString fElectronName(Names::gkElectronBrn);
292    TString fMuonName(Names::gkMuonBrn);
293    TString fInfoName(Names::gkEvtHeaderBrn);
294 +  TString fPFMetName("PFMet");
295    TString fPrimVtxName(Names::gkPVBrn);
296    TString fPFCandidateName(Names::gkPFCandidatesBrn);
297    TString fPileupInfoName(Names::gkPileupInfoBrn);
# Line 284 | Line 301 | int main(int argc, char** argv)
301    TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
302    TString fTriggerMaskName(Names::gkHltBitBrn);
303    TString fTriggerTableName(Names::gkHltTableBrn);
304 +  TString fTriggerLabelName(Names::gkHltLabelBrn);
305 +  TString fTriggerObjectName(Names::gkHltObjBrn);
306 +  TString fTriggerObjectRelsName(Form("HLTObjectsRelation"));
307 +
308 + //   TString fHLTLabName(Names::gkHltLabelBrn);
309 + //   TString fHLTTabName(Names::gkHltTableBrn);
310 + //   TString fHLTTabNamePub(Form("%sFwk",fHLTTabName.Data()));
311 + //   TString fHLTLabNamePub(Form("%sFwk",fHLTLabName.Data()));
312 + //   TString fObjsNamePub(Form("%sFwk",fTriggerObjectName.Data()));
313 +
314 +
315    
316    chain->SetBranchAddress(fInfoName,        &info);
317 +  chain->SetBranchAddress(fPFMetName,       &metArr);
318    chain->SetBranchAddress(fElectronName,    &electronArr);
319    chain->SetBranchAddress(fMuonName,        &muonArr);
320    chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
# Line 298 | Line 327 | int main(int argc, char** argv)
327    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
328    chain->SetBranchAddress(fTrackName, &trkArr);
329    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
330 +  chain->SetBranchAddress(fTriggerObjectName,  &hltObjArr);
331 +  chain->SetBranchAddress(fTriggerObjectRelsName,  &hltRelsArr);
332 +  //  chain->SetBranchAddress(fHLTTabNamePub,  &hltTab);
333 +  //  chain->SetBranchAddress(fHLTLabNamePub,  &hltLab);
334 +
335    cout << "hlttable: " << fTriggerTableName << endl;
336    hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
337 +  hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
338 +
339 +  runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);
340  
341    mithep::Vertex              vtx;          // best primary vertex in the event
342  
# Line 314 | Line 351 | int main(int argc, char** argv)
351  
352    // only 1 HLT table / file ???
353    hltchain->GetEntry(0);
317  cerr << "here... " << endl;
354  
355 <  int imax = chain->GetEntries();
356 <  cout << "nEntries: " << imax << endl;
355 >  int imax;
356 >  //-1 means process all of the events
357 >  if(ctrl.n_events_to_process != -2 && ctrl.n_events_to_process != -1){
358 >    if(ctrl.n_events_to_process > chain->GetEntries()){
359 >      cerr << "error: n_events_to_process is larger than chain->GetEntries()" << endl;
360 >      return 1;
361 >    }
362 >    imax = ctrl.n_events_to_process;
363 >  }
364 >  else
365 >    imax = chain->GetEntries();
366  
367 +  cout << "number of events to process: " << imax << endl;
368  
369    //
370    // Loop !!!!!!!!!
371    //--------------------------------------------------------------------------------------------------------------
372 +  int lastHltEntry=-2;
373    for(UInt_t ientry=0; ientry<imax; ientry++)
374      {
375        chain->GetEntry(ientry);
376        if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
377  
378 < #ifdef HACKED_RHOS
379 <      float rho = evtrhoMap[info->EvtNum()];
380 < #endif
381 <      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
382 <      setEra( fname, ctrl );
378 >      if( ctrl.debug ) {
379 >        cout << "-----------------------------------------------------------------" << endl;
380 >        cout << "-----------------------------------------------------------------" << endl;
381 >        cout << "Run: " << info->RunNum()
382 >             << "\tEvt: " << info->EvtNum()
383 >             << "\tLumi: " << info->LumiSec()
384 >             << endl;
385 >        cerr << "Run: " << info->RunNum()
386 >             << "\tEvt: " << info->EvtNum()
387 >             << "\tLumi: " << info->LumiSec()
388 >             << endl;
389 >        cout << "-----------------------------------------------------------------" << endl;
390 >      }
391  
392 +      runchain->GetEvent(info->RunEntry());
393 +      hltchain->GetEntry(runInfo->HltEntry());
394 +      hltTable->Clear();
395 +      fillTriggerNames(hltTable, hltTableStrings );
396 +      if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
397 +      lastHltEntry = runInfo->HltEntry();
398 +      fillTriggerBits( hltTable, trigMask, triggerBits );
399 +      setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
400 +      // printTriggerObjs( hltTable, fTrigObjs);
401 +      
402 +      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
403 +      if(ctrl.debug) cout << "era is " << ctrl.era << endl;
404 +      if( ctrl.era == 0  ) {
405 +        setEra( fname, ctrl );
406 +        if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
407 +      }
408  
409        //
410        // pfNoPU
# Line 344 | Line 415 | int main(int argc, char** argv)
415        assert(pfnopu_size == pfArr->GetEntries());
416  
417        //
418 +      // trigger
419 +      //
420 +      if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
421 +        currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
422 +      }
423 +      if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
424 +
425 +      //
426        // data/MC
427        //
428        if(ctrl.mc) {
# Line 358 | Line 437 | int main(int argc, char** argv)
437          //
438          weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
439          cout << "xsec weight: " << weights.w << endl;
440 <        int npu = -1;
441 <        for(int i=0;i<puArr->GetEntries();i++) {
442 <          if(puArr->At(i)->GetBunchCrossing() == 0)
443 <            npu = puArr->At(i)->GetPU_NumInteractions();
444 <        }
445 <        assert(npu>=0);
367 <        evtinfo.npu;
368 <        weights.npuw = getPUWeight(evtinfo.npu);
369 <        
440 >
441 >        unsigned npu = getNPU(puArr, 0);
442 >        weights.npu  = npu;
443 >        weights.npuw = getPUWeight(ctrl.era,npu);
444 >        cerr << "weights.npuw: " <<  weights.npuw << endl;
445 >
446          //
447          // pt reweighting for Higgs < 140, F11
448          //
# Line 378 | Line 454 | int main(int argc, char** argv)
454            }
455          }
456  
457 <        //
382 <        // trigger
383 <        //
384 <        if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
385 <          currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
386 <          hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
387 <          hltchain->GetEntry(0);
388 <          hltTable->Clear();
389 <          fillTriggerNames(hltTable, hltTableStrings );
390 <        }
391 <        if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
392 <        fillTriggerNames( hltTable, trigMask, triggerBits );
457 >        // trigger
458          passes_HLT_MC = true; // passed to apply as extern global, just for sync
459          if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
460            passes_HLT_MC = false;
461          }
462 +
463        } else {
464          //
465          // JSON
466          //
467 <        /*
468 <        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
469 <        if( !(rlrm.HasRunLumi(rl)) )  {
470 <          if( ctrl.debug ) cout << "\tfails JSON" << endl;
471 <          continue;
467 >        if(!(ctrl.noJSON) ) {
468 >          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
469 >          if( !(rlrm.HasRunLumi(rl)) )  {
470 >            if( ctrl.debug ) cout << "fails JSON" << endl;
471 >            continue;
472 >          }
473          }
474 <        */
408 <        //
474 >        
475          // trigger
476 <        //
477 <        if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
478 <          currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
479 <          hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
480 <          hltchain->GetEntry(0);
415 <          hltTable->Clear();
416 <          fillTriggerNames(hltTable, hltTableStrings );
417 <        }
418 <        if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
419 <        /*
420 <        fillTriggerNames( hltTable, trigMask, triggerBits );
421 <        if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
422 <          if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
423 <          continue;
424 <        }
425 <        */
476 >        passes_HLT = true;
477 >        //if( !passHLT(ctrl, triggerBits,ti,info->RunNum()) ) {
478 >        //  passes_HLT = false;
479 >        //}
480 >
481        }
482  
483        //
# Line 437 | Line 492 | int main(int argc, char** argv)
492          apply_HZZ4L_reference_selection(ctrl, info,
493                                          vtxArr,
494                                          pfArr,
440 #ifdef HACKED_RHOS
441                                        rho,
442 #else
495                                          puDArr,
444 #endif
496                                          electronArr,
497                                          &electronReferencePreSelection,
498                                          &electronReferenceIDMVASelectionV1,
# Line 464 | Line 515 | int main(int argc, char** argv)
515          if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
516        }  
517        
467 #ifdef SYNC
468      cout  << "bits: " << ret4l.status.selectionBits  << "\t"
469            << "Z1: " << Z1type << "\t"
470            << "Z2: " << Z2type << endl;
471      cout << endl;
472 #endif      
473      
518        if( ret4l.status.pass() ) {
519          fillAngles(ret4l,angles);
520          fillKinematics(ret4l,kinematics);
521          fillMassErrors(ret4l,muonArr,electronArr,kinematics);
522 <        fillEventInfo(info,evtinfo);
522 >        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
523 >        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
524          if( ctrl.mc) {
525 <          setEffiencyWeights(ret4l, weights);
526 <           if(ctrl.debug)
527 <             cout << "w: " << weights.w << "\t"
528 <                  << "won: " << weights.won << "\t"
529 <                  << "woff: " << weights.woff << endl;
525 >          //setEffiencyWeights(ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
526 >          if(ctrl.debug)
527 >            cout << "w: " << weights.w << "\t"
528 >                 << "npuw: " << weights.npuw << "\t"
529 >                 << "won: " << weights.won << "\t"
530 >                 << "woff: " << weights.woff << endl;
531          }
532          
533  
# Line 491 | Line 537 | int main(int argc, char** argv)
537                << "\trun: " << evtinfo.run
538                << "\tevt: " << evtinfo.evt
539                << "\tlumi: " << evtinfo.lumi
540 <              << "\tcostheta1: " << angles.costheta1
541 <              << "\tcostheta2: " << angles.costheta2
542 <              << "\tcostheta*: " << angles.costhetastar
543 <              << "\tbdt: " << angles.bdt
540 >              << "\tchannel: " << kinematics.channel
541 >              << "\tm4l: " << kinematics.m4l
542 >              << "\tmZ1: " << kinematics.mZ1
543 >              << "\tmZ2: " << kinematics.mZ2
544                << endl;
545          pass++;
546          //      if( pass > 3 ) break;
547  
548        } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
549 <                 ret4l.status.selectionBits.to_ulong() < PASS_4L) { // save for fakes ...
550 <        //      if(ctrl.debug) {
549 >                 ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
550 >        if(ctrl.debug) {
551            cout << "failingLeptons : " << failingLeptons.size() << endl;
552            for( int f=0; f<failingLeptons.size(); f++ ) {
553              SimpleLepton  l = failingLeptons[f];
# Line 510 | Line 556 | int main(int argc, char** argv)
556                   << "pT: " << l.vec.Pt()  << "\t"
557                   << "eta: " << l.vec.Eta()
558                   << endl;
559 <            //    }
559 >          }
560          }
561 +        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
562 +        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
563          foTree.Fill();
564 +      } else {
565 +        if(ctrl.debug)
566 +          cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
567        }
568 <    }  
568 >      if(ctrl.debug)
569 >        cout << "passes HLT: " << passes_HLT << endl;
570 >    }
571  
572    
573    foTree.getFile()->cd();
574    foTree.getTree()->Write();
575    nt.WriteClose();
576  
577 <  cout << endl;
578 <  cout << endl;
577 >  if(ctrl.debug)
578 >    cout << endl << endl;
579    for( int i=0; i<cutvec.size(); i++ ) {
580      cout << "cut: " << i << "\t"
581           << "pass: " << cutvec[i];
# Line 544 | Line 597 | int main(int argc, char** argv)
597   }
598  
599  
547
548 //----------------------------------------------------------------------------
549 void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
550 //----------------------------------------------------------------------------
551 {
552  vector<SimpleLepton> lepvec = evtdat.Z1leptons;
553  lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
554  double w_offline=-1, werr_offline=0;
555  double w_online=-1, werr_online=0;
556  
557  vector< pair <double,double> > wlegs; // pair here is eff & err
558  vector< pair <float,float> > mvec;  // pair here is eta & pt
559
560  //      vector< pair <float,float> > evec;  // pair here is eta & pt
561  // now deal with medium vs loose
562  vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt
563  
564  for( int k=0; k<lepvec.size(); k++ ) {
565    if( abs(lepvec[k].type) == 13 ) {
566      mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
567      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
568                                                          lepvec[k].vec.Pt() ) );
569    } else {
570      
571      // now deal with medium vs loose
572      //              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
573      
574      std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
575      evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
576      
577      //              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
578      //                                                                 lepvec[k].vec.Pt() ) );
579      
580      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
581                                                          lepvec[k].vec.Pt(),
582                                                          lepvec[k].isTight ) );
583      
584    }
585  }
586  
587  pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
588  weights.woff    = offpair.first;
589  weights.werroff = offpair.second;
590  
591  pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
592  weights.won    = onpair.first;
593  weights.werron = onpair.second;
594 }
595
596
600   //----------------------------------------------------------------------------
601 < void initRunLumiRangeMap()
601 > void initRunLumiRangeMap(ControlFlags &ctrl)
602   //----------------------------------------------------------------------------
603   {
604    /*
# Line 607 | Line 610 | void initRunLumiRangeMap()
610    // r11b
611    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
612    */
610 };
611
612
613 //----------------------------------------------------------------------------
614 void initPUWeights()
615 //----------------------------------------------------------------------------
616 {
617  TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
618  hpu = (TH1D*)(puf->Get("puWeights"));
619 }
613  
614 < //----------------------------------------------------------------------------
615 < double getPUWeight(unsigned npu)
616 < //----------------------------------------------------------------------------
617 < {
618 <  return hpu->GetBinContent(hpu->FindBin(npu));
619 < }
620 <
621 < //----------------------------------------------------------------------------
629 < void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
630 < //----------------------------------------------------------------------------
631 < {
632 <  unsigned evt;
633 <  float rho;
634 <
635 <  cout << "initialzing EvtRhoMap ";
636 <  //FILE * f = fopen("evtrho.dat", "r");
637 <  //  FILE * f = fopen("allrho.cali.uniq.dat", "r");
638 <  FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
639 <  int lcount=0;
640 <  while( fscanf( f, "%u:%f", &evt, &rho ) ) {
641 <    if( feof(f) ) break;
642 <    evtrhoMap[evt] = rho;
643 <    lcount++;
644 <    if( !(lcount%1000) ) cout << "."; flush(cout);
645 <  }
646 <  cout << " done" << endl;
647 < }
648 <
649 < //----------------------------------------------------------------------------
650 < unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
651 <                         vector<bool> &pfNoPileUpflag,
652 <                         const mithep::Array<mithep::Vertex>      * vtxArr )
653 < //----------------------------------------------------------------------------
654 < {
655 <  for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
656 <    const PFCandidate *pf = fPFCandidates->At(i);
657 <    assert(pf);
658 <
659 <    if(pf->PFType() == PFCandidate::eHadron) {
660 <      if(pf->HasTrackerTrk() &&
661 <         vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
662 <         vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
663 <
664 <        pfNoPileUpflag.push_back(1);
665 <
666 <      } else {
667 <
668 <        Bool_t vertexFound = kFALSE;
669 <        const Vertex *closestVtx = 0;
670 <        Double_t dzmin = 10000;
671 <        
672 <        // loop over vertices
673 <        for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
674 <          const Vertex *vtx = vtxArr->At(j);
675 <          assert(vtx);
676 <          
677 <          if(pf->HasTrackerTrk() &&
678 <             vtx->HasTrack(pf->TrackerTrk()) &&
679 <             vtx->TrackWeight(pf->TrackerTrk()) > 0) {
680 <            vertexFound = kTRUE;
681 <            closestVtx = vtx;
682 <            break;
683 <          }
684 <          Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
685 <          if(dz < dzmin) {
686 <            closestVtx = vtx;
687 <            dzmin = dz;
688 <          }
689 <        }
690 <
691 <        //      if(fCheckClosestZVertex) {
692 <        if(1) {
693 <          // Fallback: if track is not associated with any vertex,
694 <          // associate it with the vertex closest in z
695 <          if(vertexFound || closestVtx != vtxArr->At(0)) {
696 <            //      pfPileUp->Add(pf);
697 <            pfNoPileUpflag.push_back(0);
698 <          } else {
699 <            pfNoPileUpflag.push_back(1);
700 <          }
701 <        } else {
702 <          if(vertexFound && closestVtx != vtxArr->At(0)) {
703 <            //      pfPileUp->Add(pf);
704 <            pfNoPileUpflag.push_back(0);
705 <          } else {
706 <            //      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
707 <            pfNoPileUpflag.push_back(1);
708 <          }
709 <        }
710 <      } //hadron & trk stuff
711 <    } else { // hadron
712 <      //      PFCandidate * ptr = pfNoPileUp->AddNew();
713 <      pfNoPileUpflag.push_back(1);
714 <    }
715 <  } // Loop over PF candidates
716 <
717 <  return pfNoPileUpflag.size();
718 < }
719 <
720 < void setEra(string fname, ControlFlags &ctrl)
721 < {
722 <  if( ctrl.mc && (fname.find("f11-") != string::npos)  ) ctrl.era=2011 ;
723 <  if( ctrl.mc && (fname.find("s12-") != string::npos)  ) ctrl.era=2012 ;
724 <  if( !ctrl.mc && (fname.find("r11-") != string::npos)  ) ctrl.era=2011 ;
725 <  if( !ctrl.mc && (fname.find("r12-") != string::npos)  ) ctrl.era=2012 ;
726 < }
614 >  // 2012 only for now ...
615 >  std::string jsonFile;
616 >  if(ctrl.jsonFile=="")
617 >    jsonFile = "./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt";
618 >  else
619 >    jsonFile = ctrl.jsonFile;
620 >  rlrm.AddJSONFile(jsonFile);  
621 > };

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines