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.55 by anlevin, Wed Sep 12 12:13:27 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 65 | Line 69 | using namespace std;
69   #include "RunLumiRangeMap.h"
70   #include "SampleWeight.h"
71   #include "EfficiencyWeightsInterface.h"
72 + #include "SelectionFuncs.h"
73  
74   #include "SimpleLepton.h"
75  
# Line 77 | Line 82 | using namespace mithep;
82   //
83   // globals
84   //----------------------------------------------------------------------------
85 < TH1D * hpu;
85 > TH1D *hpu_2011, *hpu_2012;
86   RunLumiRangeMap rlrm;
87   vector<SimpleLepton> failingLeptons ; // for fake estimation
88   vector<SimpleLepton> passingLeptons;      // for fake estimation
# Line 86 | Line 91 | vector<vector<unsigned> > zcutvec;
91   vector<vector<unsigned> > zzcutvec;
92   map<unsigned,float>       evtrhoMap;
93   vector<string> cutstrs;
94 < bool passes_HLT_MC;
94 > bool passes_HLT_MC,passes_HLT;
95   vector<bool>   PFnoPUflag;;
96  
97   //
98   // prototypes
99   //----------------------------------------------------------------------------
100 < 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 );
100 > void initRunLumiRangeMap(ControlFlags &ctrl);
101   //----------------------------------------------------------------------------
102  
103  
# 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 +  TrigInfo ti;
256 +  initAnalysisTriggers(ti);
257  
250  // hack
251  initEvtRhoMap(evtrhoMap);
252
253  
258    //
259    // Setup tree I/O
260    //--------------------------------------------------------------------------------------------------------------
# Line 258 | Line 262 | int main(int argc, char** argv)
262    TTree *eventTree=0;  
263    string currentFile("");
264  
265 <  mithep::EventHeader *info    = new mithep::EventHeader();
265 >  UInt_t fNMaxTriggers = 1024;
266 >  mithep::EventHeader                         *info          = new mithep::EventHeader();
267 >  mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
268    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
269    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
270    mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
# Line 269 | Line 275 | int main(int argc, char** argv)
275    mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
276    mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
277    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
278 <  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
278 >  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable(fNMaxTriggers);
279    vector<string>                              *hltTableStrings  = new vector<string>();
280 +  vector<string>                              *hltLabelStrings  = new vector<string>();
281 +  mithep::Array<mithep::TriggerObject>        *hltObjArr     = new mithep::Array<mithep::TriggerObject>();
282 +  mithep::Array<mithep::TriggerObjectRel>     *hltRelsArr    = new mithep::Array<mithep::TriggerObjectRel>();
283 +  std::vector<std::string>                    *hltTab        = new vector<string>();
284 +  std::vector<std::string>                    *hltLab        = new vector<string>();
285 +  TriggerObjectsTable                         *fTrigObjs     = new TriggerObjectsTable(hltTable,fNMaxTriggers);
286 +
287 +  RunInfo                                     *runInfo       = new RunInfo();
288  
289    TString fElectronName(Names::gkElectronBrn);
290    TString fMuonName(Names::gkMuonBrn);
291    TString fInfoName(Names::gkEvtHeaderBrn);
292 +  TString fPFMetName("PFMet");
293    TString fPrimVtxName(Names::gkPVBrn);
294    TString fPFCandidateName(Names::gkPFCandidatesBrn);
295    TString fPileupInfoName(Names::gkPileupInfoBrn);
# Line 284 | Line 299 | int main(int argc, char** argv)
299    TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
300    TString fTriggerMaskName(Names::gkHltBitBrn);
301    TString fTriggerTableName(Names::gkHltTableBrn);
302 +  TString fTriggerLabelName(Names::gkHltLabelBrn);
303 +  TString fTriggerObjectName(Names::gkHltObjBrn);
304 +  TString fTriggerObjectRelsName(Form("HLTObjectsRelation"));
305 +
306 + //   TString fHLTLabName(Names::gkHltLabelBrn);
307 + //   TString fHLTTabName(Names::gkHltTableBrn);
308 + //   TString fHLTTabNamePub(Form("%sFwk",fHLTTabName.Data()));
309 + //   TString fHLTLabNamePub(Form("%sFwk",fHLTLabName.Data()));
310 + //   TString fObjsNamePub(Form("%sFwk",fTriggerObjectName.Data()));
311 +
312 +
313    
314    chain->SetBranchAddress(fInfoName,        &info);
315 +  chain->SetBranchAddress(fPFMetName,       &metArr);
316    chain->SetBranchAddress(fElectronName,    &electronArr);
317    chain->SetBranchAddress(fMuonName,        &muonArr);
318    chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
# Line 298 | Line 325 | int main(int argc, char** argv)
325    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
326    chain->SetBranchAddress(fTrackName, &trkArr);
327    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
328 +  chain->SetBranchAddress(fTriggerObjectName,  &hltObjArr);
329 +  chain->SetBranchAddress(fTriggerObjectRelsName,  &hltRelsArr);
330 +  //  chain->SetBranchAddress(fHLTTabNamePub,  &hltTab);
331 +  //  chain->SetBranchAddress(fHLTLabNamePub,  &hltLab);
332 +
333    cout << "hlttable: " << fTriggerTableName << endl;
334    hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
335 +  hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
336 +
337 +  runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);
338  
339    mithep::Vertex              vtx;          // best primary vertex in the event
340  
# Line 314 | Line 349 | int main(int argc, char** argv)
349  
350    // only 1 HLT table / file ???
351    hltchain->GetEntry(0);
317  cerr << "here... " << endl;
352  
353 <  int imax = chain->GetEntries();
354 <  cout << "nEntries: " << imax << endl;
353 >  int imax;
354 >  //-1 means process all of the events
355 >  if(ctrl.n_events_to_process != -2 && ctrl.n_events_to_process != -1){
356 >    if(ctrl.n_events_to_process > chain->GetEntries()){
357 >      cerr << "error: n_events_to_process is larger than chain->GetEntries()" << endl;
358 >      return 1;
359 >    }
360 >    imax = ctrl.n_events_to_process;
361 >  }
362 >  else
363 >    imax = chain->GetEntries();
364  
365 +  cout << "number of events to process: " << imax << endl;
366  
367    //
368    // Loop !!!!!!!!!
369    //--------------------------------------------------------------------------------------------------------------
370 +  int lastHltEntry=-2;
371    for(UInt_t ientry=0; ientry<imax; ientry++)
372      {
373        chain->GetEntry(ientry);
374        if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
375  
376 < #ifdef HACKED_RHOS
377 <      float rho = evtrhoMap[info->EvtNum()];
378 < #endif
379 <      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
380 <      setEra( fname, ctrl );
376 >      if( ctrl.debug ) {
377 >        cout << "-----------------------------------------------------------------" << endl;
378 >        cout << "-----------------------------------------------------------------" << endl;
379 >        cout << "Run: " << info->RunNum()
380 >             << "\tEvt: " << info->EvtNum()
381 >             << "\tLumi: " << info->LumiSec()
382 >             << endl;
383 >        cerr << "Run: " << info->RunNum()
384 >             << "\tEvt: " << info->EvtNum()
385 >             << "\tLumi: " << info->LumiSec()
386 >             << endl;
387 >        cout << "-----------------------------------------------------------------" << endl;
388 >      }
389  
390 +      runchain->GetEvent(info->RunEntry());
391 +      hltchain->GetEntry(runInfo->HltEntry());
392 +      hltTable->Clear();
393 +      fillTriggerNames(hltTable, hltTableStrings );
394 +      if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
395 +      lastHltEntry = runInfo->HltEntry();
396 +      fillTriggerBits( hltTable, trigMask, triggerBits );
397 +      setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
398 +      // printTriggerObjs( hltTable, fTrigObjs);
399 +      
400 +      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
401 +      if(ctrl.debug) cout << "era is " << ctrl.era << endl;
402 +      if( ctrl.era == 0  ) {
403 +        setEra( fname, ctrl );
404 +        if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
405 +      }
406  
407        //
408        // pfNoPU
# Line 344 | Line 413 | int main(int argc, char** argv)
413        assert(pfnopu_size == pfArr->GetEntries());
414  
415        //
416 +      // trigger
417 +      //
418 +      if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
419 +        currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
420 +      }
421 +      if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
422 +
423 +      //
424        // data/MC
425        //
426        if(ctrl.mc) {
# Line 358 | Line 435 | int main(int argc, char** argv)
435          //
436          weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
437          cout << "xsec weight: " << weights.w << endl;
438 <        int npu = -1;
439 <        for(int i=0;i<puArr->GetEntries();i++) {
440 <          if(puArr->At(i)->GetBunchCrossing() == 0)
441 <            npu = puArr->At(i)->GetPU_NumInteractions();
442 <        }
443 <        assert(npu>=0);
367 <        evtinfo.npu;
368 <        weights.npuw = getPUWeight(evtinfo.npu);
369 <        
438 >
439 >        unsigned npu = getNPU(puArr, 0);
440 >        weights.npu  = npu;
441 >        weights.npuw = getPUWeight(ctrl.era,npu);
442 >        cerr << "weights.npuw: " <<  weights.npuw << endl;
443 >
444          //
445          // pt reweighting for Higgs < 140, F11
446          //
# Line 378 | Line 452 | int main(int argc, char** argv)
452            }
453          }
454  
455 <        //
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 );
455 >        // trigger
456          passes_HLT_MC = true; // passed to apply as extern global, just for sync
457          if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
458            passes_HLT_MC = false;
459          }
460 +
461        } else {
462          //
463          // JSON
464          //
465 <        /*
466 <        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
467 <        if( !(rlrm.HasRunLumi(rl)) )  {
468 <          if( ctrl.debug ) cout << "\tfails JSON" << endl;
469 <          continue;
465 >        if(!(ctrl.noJSON) ) {
466 >          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
467 >          if( !(rlrm.HasRunLumi(rl)) )  {
468 >            if( ctrl.debug ) cout << "fails JSON" << endl;
469 >            continue;
470 >          }
471          }
472 <        */
408 <        //
472 >        
473          // trigger
474 <        //
475 <        if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
476 <          currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
477 <          hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
478 <          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 <        */
474 >        passes_HLT = true;
475 >        //if( !passHLT(ctrl, triggerBits,ti,info->RunNum()) ) {
476 >        //  passes_HLT = false;
477 >        //}
478 >
479        }
480  
481        //
# Line 437 | Line 490 | int main(int argc, char** argv)
490          apply_HZZ4L_reference_selection(ctrl, info,
491                                          vtxArr,
492                                          pfArr,
440 #ifdef HACKED_RHOS
441                                        rho,
442 #else
493                                          puDArr,
444 #endif
494                                          electronArr,
495                                          &electronReferencePreSelection,
496                                          &electronReferenceIDMVASelectionV1,
# Line 464 | Line 513 | int main(int argc, char** argv)
513          if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
514        }  
515        
467 #ifdef SYNC
468      cout  << "bits: " << ret4l.status.selectionBits  << "\t"
469            << "Z1: " << Z1type << "\t"
470            << "Z2: " << Z2type << endl;
471      cout << endl;
472 #endif      
473      
516        if( ret4l.status.pass() ) {
517          fillAngles(ret4l,angles);
518          fillKinematics(ret4l,kinematics);
519          fillMassErrors(ret4l,muonArr,electronArr,kinematics);
520 <        fillEventInfo(info,evtinfo);
520 >        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
521 >        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
522          if( ctrl.mc) {
523 <          setEffiencyWeights(ret4l, weights);
524 <           if(ctrl.debug)
525 <             cout << "w: " << weights.w << "\t"
526 <                  << "won: " << weights.won << "\t"
527 <                  << "woff: " << weights.woff << endl;
523 >          //setEffiencyWeights(ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
524 >          if(ctrl.debug)
525 >            cout << "w: " << weights.w << "\t"
526 >                 << "npuw: " << weights.npuw << "\t"
527 >                 << "won: " << weights.won << "\t"
528 >                 << "woff: " << weights.woff << endl;
529          }
530          
531  
# Line 491 | Line 535 | int main(int argc, char** argv)
535                << "\trun: " << evtinfo.run
536                << "\tevt: " << evtinfo.evt
537                << "\tlumi: " << evtinfo.lumi
538 <              << "\tcostheta1: " << angles.costheta1
539 <              << "\tcostheta2: " << angles.costheta2
540 <              << "\tcostheta*: " << angles.costhetastar
541 <              << "\tbdt: " << angles.bdt
538 >              << "\tchannel: " << kinematics.channel
539 >              << "\tm4l: " << kinematics.m4l
540 >              << "\tmZ1: " << kinematics.mZ1
541 >              << "\tmZ2: " << kinematics.mZ2
542                << endl;
543          pass++;
544          //      if( pass > 3 ) break;
545  
546        } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
547 <                 ret4l.status.selectionBits.to_ulong() < PASS_4L) { // save for fakes ...
548 <        //      if(ctrl.debug) {
547 >                 ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
548 >        if(ctrl.debug) {
549            cout << "failingLeptons : " << failingLeptons.size() << endl;
550            for( int f=0; f<failingLeptons.size(); f++ ) {
551              SimpleLepton  l = failingLeptons[f];
# Line 510 | Line 554 | int main(int argc, char** argv)
554                   << "pT: " << l.vec.Pt()  << "\t"
555                   << "eta: " << l.vec.Eta()
556                   << endl;
557 <            //    }
557 >          }
558          }
559 +        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
560 +        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
561          foTree.Fill();
562 +      } else {
563 +        if(ctrl.debug)
564 +          cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
565        }
566 <    }  
566 >      if(ctrl.debug)
567 >        cout << "passes HLT: " << passes_HLT << endl;
568 >    }
569  
570    
571    foTree.getFile()->cd();
572    foTree.getTree()->Write();
573    nt.WriteClose();
574  
575 <  cout << endl;
576 <  cout << endl;
575 >  if(ctrl.debug)
576 >    cout << endl << endl;
577    for( int i=0; i<cutvec.size(); i++ ) {
578      cout << "cut: " << i << "\t"
579           << "pass: " << cutvec[i];
# Line 544 | Line 595 | int main(int argc, char** argv)
595   }
596  
597  
547
598   //----------------------------------------------------------------------------
599 < 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 <
597 < //----------------------------------------------------------------------------
598 < void initRunLumiRangeMap()
599 > void initRunLumiRangeMap(ControlFlags &ctrl)
600   //----------------------------------------------------------------------------
601   {
602    /*
# Line 607 | Line 608 | void initRunLumiRangeMap()
608    // r11b
609    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
610    */
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 }
620
621 //----------------------------------------------------------------------------
622 double getPUWeight(unsigned npu)
623 //----------------------------------------------------------------------------
624 {
625  return hpu->GetBinContent(hpu->FindBin(npu));
626 }
627
628 //----------------------------------------------------------------------------
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
611  
612 <  return pfNoPileUpflag.size();
613 < }
614 <
615 < void setEra(string fname, ControlFlags &ctrl)
616 < {
617 <  if( ctrl.mc && (fname.find("f11-") != string::npos)  ) ctrl.era=2011 ;
618 <  if( ctrl.mc && (fname.find("s12-") != string::npos)  ) ctrl.era=2012 ;
619 <  if( !ctrl.mc && (fname.find("r11-") != string::npos)  ) ctrl.era=2011 ;
725 <  if( !ctrl.mc && (fname.find("r12-") != string::npos)  ) ctrl.era=2012 ;
726 < }
612 >  // 2012 only for now ...
613 >  std::string jsonFile;
614 >  if(ctrl.jsonFile=="")
615 >    jsonFile = "./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt";
616 >  else
617 >    jsonFile = ctrl.jsonFile;
618 >  rlrm.AddJSONFile(jsonFile);  
619 > };

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines