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.23 by khahn, Sun May 20 19:08:59 2012 UTC vs.
Revision 1.49 by khahn, Fri Jun 22 03:10:32 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 50 | Line 54 | using namespace std;
54   #include "MuonSelection.h"
55   #include "ElectronSelection.h"
56   #include "IsolationSelection.h"
57 <
54 < //#include "Selection.h"
55 < //#include "ReferenceSelection.h"
56 < #include "ReferenceSelection2.h"
57 > #include "ReferenceSelection.h"
58  
59   #include "TriggerUtils.h"
60   #include "PassHLT.h"
61   #include "Angles.h"
62   #include "KinematicsStruct.h"
63   #include "InfoStruct.h"
64 < //#include "GenInfoStruct.h"
64 > #include "GenInfoStruct.h"
65   #include "WeightStruct.h"
65 //#include "GlobalDataAndFuncs.h"
66 #include "RunLumiRangeMap.h"
67
66   #include "AngleTuple.h"
67   #include "FOTuple.h"
68  
69 + #include "RunLumiRangeMap.h"
70   #include "SampleWeight.h"
71   #include "EfficiencyWeightsInterface.h"
72 + #include "SelectionFuncs.h"
73  
74   #include "SimpleLepton.h"
75  
# Line 82 | 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 97 | Line 97 | vector<bool>   PFnoPUflag;;
97   //
98   // prototypes
99   //----------------------------------------------------------------------------
100 < bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
101 < void initPUWeights();
102 < double getPUWeight(unsigned npu);
103 < void setEffiencyWeights(EventData&, WeightStruct& );
104 < void initRunLumiRangeMap();
105 < void initEvtRhoMap(map<unsigned,float> &);
106 < unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
107 <                         vector<bool> &pfNoPileUpFlag,
108 <                         mithep::Array<Vertex>      * vtxArr );
100 > void initRunLumiRangeMap(ControlFlags &ctrl);
101   //----------------------------------------------------------------------------
102  
103  
# Line 117 | Line 109 | int main(int argc, char** argv)
109   {
110    cutstrs.push_back(string("skim0"));              //0
111    cutstrs.push_back(string("skim1"));              //1
112 +  cutstrs.push_back(string("skim2"));              //1
113    cutstrs.push_back(string("trigger"));            //2
114    // -------------------------------------------------
115    cutstrs.push_back(string("Z candidate"));        //3
# Line 128 | Line 121 | int main(int argc, char** argv)
121    cutstrs.push_back(string("ZZ 20/10"));           //8
122    cutstrs.push_back(string("resonance"));          //9
123    cutstrs.push_back(string("m4l>70"));             //10
124 <  cutstrs.push_back(string("m4l>100"));            //11
124 >  cutstrs.push_back(string("m4l>100, mZ2 > 12"));                  //11
125  
126  
127  
# Line 175 | Line 168 | int main(int argc, char** argv)
168    TChain * hltchain = new TChain("HLT");
169  
170    string fname;
171 +  unsigned total_unskimmed=0;
172    if( !ctrl.inputfiles.empty() ) {
173      ifstream f(ctrl.inputfiles.c_str());
174      while (f >> fname) {
# Line 182 | Line 176 | int main(int argc, char** argv)
176        cout << "adding inputfile : " << fname.c_str() << endl;
177        entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
178        cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
179 +      total_unskimmed += entrymap[string(fname.c_str())];
180        chain->AddFile(fname.c_str());
181        hltchain->AddFile(fname.c_str());
182      }
183    } else {
184      cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
185 <    unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
186 <    cout << "tmpent: " << tmpent << endl;
187 <    entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
188 <    cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
185 >    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
186 >    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
187 >    cout << "unskimmed entries: " << unsk_ents << endl;
188 >    total_unskimmed += unsk_ents;
189      chain->AddFile(ctrl.inputfile.c_str());
190      hltchain->AddFile(ctrl.inputfile.c_str());
191    }
192 +  // // write the total number of unskimmed events that went into making this output file to a text file
193 +  // writeEntries(ctrl,total_unskimmed);
194  
195    const char * ofname;
196    if( strcmp( ctrl.outputfile.c_str(), "") ) {
# Line 207 | Line 204 | int main(int argc, char** argv)
204    SimpleTable xstab(xspath.c_str());
205  
206  
207 +
208    //
209    // Setup
210    //--------------------------------------------------------------------------------------------------------------
# Line 214 | Line 212 | int main(int argc, char** argv)
212    KinematicsStruct kinematics;
213    InfoStruct evtinfo;
214    WeightStruct weights;
215 <  //  GenInfoStruct geninfo;
215 >  GenInfoStruct geninfo;
216  
217    AngleTuple nt( (const char*)ofname, (const char*)"zznt");
218    nt.makeAngleBranch(angles);
219    nt.makeKinematicsBranch(kinematics);
220    nt.makeInfoBranch(evtinfo);
221  
222 +  
223    FOTuple foTree( nt.getFile(), (const char*)"FOtree");
224    foTree.makeFailedLeptonBranch(failingLeptons);
225 <  foTree.makeZLeptonBranch(passingLeptons);
225 >  foTree.makePassedLeptonBranch(passingLeptons);
226  
227 <  TH1F * h_w_hpt;
227 >  TH1F * h_wf11_hpt;
228    if(ctrl.mc) {
229      nt.makeWeightBranch(weights);
230 <    //    nt.makeGenInfoBranch(geninfo);
230 >    if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
231      initEfficiencyWeights();
232      initPUWeights();
233 <    /*
233 >
234      // Higgs only, pt reweighting
235      if( ctrl.mH <= 140 ) {
236        char wptname[256];
# Line 239 | Line 238 | int main(int argc, char** argv)
238        TFile * f_hpt = new TFile(wptname);
239        f_hpt->Print();
240        sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
241 <      h_w_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
243 <      h_w_hpt->Print();
241 >      h_wf11_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
242      }
243 <    */
243 >
244    } else {
245      initRunLumiRangeMap();
246    }
247  
248    //  initMuonIDMVA();
249 <  initMuonIsoMVA();
252 <  initElectronIDMVA();
253 <  initElectronIsoMVA();
249 >  initElectronIDMVAV1();
250    initTrigger();
251  
256  // hack
257  initEvtRhoMap(evtrhoMap);
258
259  
252    //
253    // Setup tree I/O
254    //--------------------------------------------------------------------------------------------------------------
# Line 264 | Line 256 | int main(int argc, char** argv)
256    TTree *eventTree=0;  
257    string currentFile("");
258  
259 +  UInt_t fNMaxTriggers = 1024;
260    mithep::EventHeader *info    = new mithep::EventHeader();
261 <  //  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
261 >  mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
262    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
263    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
264    mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
# Line 274 | Line 267 | int main(int argc, char** argv)
267    mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
268    mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
269    mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
270 +  mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
271    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
272 <  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
272 >  mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable(fNMaxTriggers);
273    vector<string>                              *hltTableStrings  = new vector<string>();
274 +  vector<string>                              *hltLabelStrings  = new vector<string>();
275 +  mithep::Array<mithep::TriggerObject>        *hltObjArr     = new mithep::Array<mithep::TriggerObject>();
276 +  mithep::Array<mithep::TriggerObjectRel>     *hltRelsArr    = new mithep::Array<mithep::TriggerObjectRel>();
277 +  std::vector<std::string>                    *hltTab        = new vector<string>();
278 +  std::vector<std::string>                    *hltLab        = new vector<string>();
279 +  TriggerObjectsTable                         *fTrigObjs     = new TriggerObjectsTable(hltTable,fNMaxTriggers);
280 +
281  
282    TString fElectronName(Names::gkElectronBrn);
283    TString fMuonName(Names::gkMuonBrn);
284    TString fInfoName(Names::gkEvtHeaderBrn);
285 +  TString fPFMetName("PFMet");
286    TString fPrimVtxName(Names::gkPVBrn);
287    TString fPFCandidateName(Names::gkPFCandidatesBrn);
288    TString fPileupInfoName(Names::gkPileupInfoBrn);
289    TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
290    TString fTrackName(Names::gkTrackBrn);
291    TString fMCParticleName(Names::gkMCPartBrn);
292 +  TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
293    TString fTriggerMaskName(Names::gkHltBitBrn);
294    TString fTriggerTableName(Names::gkHltTableBrn);
295 +  TString fTriggerLabelName(Names::gkHltLabelBrn);
296 +  TString fTriggerObjectName(Names::gkHltObjBrn);
297 +  TString fTriggerObjectRelsName(Form("HLTObjectsRelation"));
298 +
299 + //   TString fHLTLabName(Names::gkHltLabelBrn);
300 + //   TString fHLTTabName(Names::gkHltTableBrn);
301 + //   TString fHLTTabNamePub(Form("%sFwk",fHLTTabName.Data()));
302 + //   TString fHLTLabNamePub(Form("%sFwk",fHLTLabName.Data()));
303 + //   TString fObjsNamePub(Form("%sFwk",fTriggerObjectName.Data()));
304 +
305 +
306    
307    chain->SetBranchAddress(fInfoName,        &info);
308 +  chain->SetBranchAddress(fPFMetName,       &metArr);
309    chain->SetBranchAddress(fElectronName,    &electronArr);
310    chain->SetBranchAddress(fMuonName,        &muonArr);
311    chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
# Line 298 | Line 313 | int main(int argc, char** argv)
313    if( ctrl.mc ) {
314      chain->SetBranchAddress(fPileupInfoName,  &puArr);
315      chain->SetBranchAddress(fMCParticleName,  &mcArr);
316 +    chain->SetBranchAddress(fMCEvtInfoName,  &mcEvtInfo);
317    }
318    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
319    chain->SetBranchAddress(fTrackName, &trkArr);
320    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
321 +  chain->SetBranchAddress(fTriggerObjectName,  &hltObjArr);
322 +  chain->SetBranchAddress(fTriggerObjectRelsName,  &hltRelsArr);
323 +  //  chain->SetBranchAddress(fHLTTabNamePub,  &hltTab);
324 +  //  chain->SetBranchAddress(fHLTLabNamePub,  &hltLab);
325 +
326    cout << "hlttable: " << fTriggerTableName << endl;
327    hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
328 +  hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
329  
330    mithep::Vertex              vtx;          // best primary vertex in the event
331  
# Line 318 | Line 340 | int main(int argc, char** argv)
340  
341    // only 1 HLT table / file ???
342    hltchain->GetEntry(0);
321  cerr << "here... " << endl;
343  
344    int imax = chain->GetEntries();
345    cout << "nEntries: " << imax << endl;
346  
326
347    //
348    // Loop !!!!!!!!!
349    //--------------------------------------------------------------------------------------------------------------
350    for(UInt_t ientry=0; ientry<imax; ientry++)
351      {
352        chain->GetEntry(ientry);
353 +      if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
354  
355 <      /*
356 <      bool found_tau=false;
357 <      if( ctrl.mc) {  // && ctrl.debug
358 <        cout << "nparticles: " << mcArr->GetEntries() << endl;
359 <        cout << "-----------------------------------" << endl;
360 <        for( int i=0; i<mcArr->GetEntries(); i++ ) {
361 <          const mithep::MCParticle *mc = mcArr->At(i);
362 <          cout << "i: " << i << "\t"
342 <               << "id: " << mc->AbsPdgId()
343 <               << endl;
344 <          if( mc->Is(MCParticle::kTau) ) {
345 <            found_tau=true;
346 <          }
347 <        }
355 >      if( ctrl.debug ) {
356 >        cout << "-----------------------------------------------------------------" << endl;
357 >        cout << "-----------------------------------------------------------------" << endl;
358 >        cout << "Run: " << info->RunNum()
359 >             << "\tEvt: " << info->EvtNum()
360 >             << "\tLumi: " << info->LumiSec()
361 >             << endl;
362 >        cout << "-----------------------------------------------------------------" << endl;
363        }
349      cout << "-----------------------------------" << endl;
350      cout << "found tau: " << found_tau << endl;
351      */
352
353      // events that fail bc ele overlaps w/ reco mu
354      /*
355      if( info->EvtNum() != 67315 &&
356          info->EvtNum() != 288693 ) continue;
357      */
358
359      // events where they apparently don't veto electrons ...
360      /*
361      if( info->EvtNum() != 141713 &&
362          info->EvtNum() != 178019 &&
363          info->EvtNum() != 60823 &&
364          info->EvtNum() != 192795 ) continue;
365      */
364  
367      if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
368      setPV( ctrl, vtxArr, vtx);
369      //if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
365  
366 <      float rho = evtrhoMap[info->EvtNum()];
366 >      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
367 >      if(ctrl.debug) cout << "era is " << ctrl.era << endl;
368 >      if( ctrl.era == 0  ) {
369 >        setEra( fname, ctrl );
370 >        if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
371 >      }
372  
373        //
374        // pfNoPU
# Line 383 | Line 383 | int main(int argc, char** argv)
383        //
384        if(ctrl.mc) {
385          //
386 +        // gen info
387 +        //
388 +        if( ctrl.fillGen )
389 +          fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
390 +
391 +        //
392          // xsec & PU weights
393          //
394          weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
395          cout << "xsec weight: " << weights.w << endl;
396 <        int npu = -1;
397 <        for(int i=0;i<puArr->GetEntries();i++) {
398 <          if(puArr->At(i)->GetBunchCrossing() == 0)
399 <            npu = puArr->At(i)->GetPU_NumInteractions();
396 >
397 >        unsigned npu = getNPU(puArr, 0);
398 >        weights.npu  = npu;
399 >        weights.npuw = getPUWeight(ctrl.era,npu);
400 >        cerr << "weights.npuw: " <<  weights.npuw << endl;
401 >
402 >        //
403 >        // pt reweighting for Higgs < 140, F11
404 >        //
405 >        if( ctrl.fillGen ) {
406 >          if( (fname.find("f11-h115zz4l") != string::npos) ||
407 >              (fname.find("f11-h120zz4l") != string::npos) ||
408 >              (fname.find("f11-h130zz4l") != string::npos)  ) {
409 >            weights.w *= h_wf11_hpt->GetBinContent(h_wf11_hpt->FindBin(geninfo.pt_zz));
410 >          }
411          }
412 <        assert(npu>=0);
396 <        evtinfo.npu;
397 <        weights.npuw = getPUWeight(evtinfo.npu);
398 <        //      cout << "weight: " << weights.w << endl;
412 >
413          //
414          // trigger
415          //
416          if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
417            currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
418            hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
419 +          hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
420            hltchain->GetEntry(0);
421            hltTable->Clear();
422 <          fillTriggerBits(hltTable, hltTableStrings );
422 >          fillTriggerNames(hltTable, hltTableStrings );
423          }
424          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
425          fillTriggerBits( hltTable, trigMask, triggerBits );
426          passes_HLT_MC = true; // passed to apply as extern global, just for sync
427 <        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) {
427 >        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
428            passes_HLT_MC = false;
429          }
430 +
431 +        setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
432 +
433        } else {
434          //
435          // JSON
436          //
437 <        /*
438 <        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
439 <        if( !(rlrm.HasRunLumi(rl)) )  {
440 <          if( ctrl.debug ) cout << "\tfails JSON" << endl;
441 <          continue;
437 >        if(!(ctrl.noJSON) ) {
438 >          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
439 >          if( !(rlrm.HasRunLumi(rl)) )  {
440 >            if( ctrl.debug ) cout << "fails JSON" << endl;
441 >            continue;
442 >          }
443          }
444 <        */
444 >        
445          //
446          // trigger
447          //
448          if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
449            currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
450            hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
451 +          hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
452            hltchain->GetEntry(0);
453            hltTable->Clear();
454 <          fillTriggerBits(hltTable, hltTableStrings );
454 >          fillTriggerNames(hltTable, hltTableStrings );
455          }
456 +
457 +        setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
458 +
459          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
460          /*
461          fillTriggerBits( hltTable, trigMask, triggerBits );
# Line 452 | Line 475 | int main(int argc, char** argv)
475  
476  
477          // reference
478 <        /*
479 <        apply_HZZ4L_reference2_selection(ctrl, info,
480 <                                         vtx,
481 <                                         pfArr,
482 <                                         //rho,
483 <                                         puDArr,
484 <                                         electronArr,
485 <                                         &electronReferencePreSelection,
486 <                                         &electronReferenceIDMVASelection,
487 <                                         &electronReferenceIsoSelection,
488 <                                         muonArr,
489 <                                         &muonReferencePreSelection,
490 <                                         &muonIDPFSelection,
468 <                                         &muonReferenceIsoSelection);
469 <        */
470 <
471 <        // evolved
472 <        apply_HZZ4L_reference2_selection(ctrl, info,
473 <                                         vtx,
474 <                                         pfArr,
475 <                                         //rho,
476 <                                         puDArr,
477 <                                         electronArr,
478 <                                         &electronReferencePreSelection,
479 <                                         //&electronPreSelection,
480 <                                         &electronReferenceIDMVASelection,
481 <                                         &electronIsoMVASelection,
482 <                                         muonArr,
483 <                                         &muonReferencePreSelection,
484 <                                         //&muonPreSelection,
485 <                                         &muonIDPFSelection,
486 <                                         &muonIsoMVASelection);
487 <
488 <
478 >        apply_HZZ4L_reference_selection(ctrl, info,
479 >                                        vtxArr,
480 >                                        pfArr,
481 >                                        puDArr,
482 >                                        electronArr,
483 >                                        &electronReferencePreSelection,
484 >                                        &electronReferenceIDMVASelectionV1,
485 >                                        &electronReferenceIsoSelection,
486 >                                        muonArr,
487 >                                        &muonReferencePreSelection,
488 >                                        &muonIDPFSelection,
489 >                                        &muonReferenceIsoSelection);
490 >      
491        if( ctrl.debug ) cout << endl;
492  
493  
# Line 500 | Line 502 | int main(int argc, char** argv)
502        }  
503        
504   #ifdef SYNC
505 <      cout  << "bits: " << ret4l.status.selectionBits  << "\t"
506 <            << "Z1: " << Z1type << "\t"
507 <            << "Z2: " << Z2type << endl;
508 <      cout << endl;
505 >      if(ctrl.debug)
506 >        cout  << "bits: " << ret4l.status.selectionBits  << "\t"
507 >              << "Z1: " << Z1type << "\t"
508 >              << "Z2: " << Z2type << endl << endl;
509   #endif      
510        
511        if( ret4l.status.pass() ) {
512          fillAngles(ret4l,angles);
513          fillKinematics(ret4l,kinematics);
514 <        fillEventInfo(info,evtinfo);
514 >        fillMassErrors(ret4l,muonArr,electronArr,kinematics);
515 >        TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
516 >        fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
517          if( ctrl.mc) {
518 <        // fillGenInfo(ginfo,geninfo);
519 <          setEffiencyWeights(ret4l, weights);
520 <           if(ctrl.debug)
521 <             cout << "w: " << weights.w << "\t"
522 <                  << "won: " << weights.won << "\t"
523 <                  << "woff: " << weights.woff << endl;
518 >          //      if( std::string(ctrl.efftype) != std::string("WTF?")) {
519 >            setEffiencyWeights(ctrl.era, ret4l, weights);
520 >            //    } else {
521 >            //      weights.won = weights.woff = 1.;
522 >            //    }
523 >          if(ctrl.debug)
524 >            cout << "w: " << weights.w << "\t"
525 >                 << "npuw: " << weights.npuw << "\t"
526 >                 << "won: " << weights.won << "\t"
527 >                 << "woff: " << weights.woff << endl;
528          }
529          
530 <        /*
523 <        // only for Higgs < 140
524 <        geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz));
525 <        angles.bdt = reader->EvaluateMVA("BDTG");
526 <        */
530 >
531          nt.Fill();
532          
533          cerr  << "PASS:: "
534                << "\trun: " << evtinfo.run
535                << "\tevt: " << evtinfo.evt
536                << "\tlumi: " << evtinfo.lumi
537 <              << "\tcostheta1: " << angles.costheta1
538 <              << "\tcostheta2: " << angles.costheta2
539 <              << "\tcostheta*: " << angles.costhetastar
540 <              << "\tbdt: " << angles.bdt
537 >              << "\tchannel: " << kinematics.channel
538 >              << "\tm4l: " << kinematics.m4l
539 >              << "\tmZ1: " << kinematics.mZ1
540 >              << "\tmZ2: " << kinematics.mZ2
541                << endl;
542          pass++;
543          //      if( pass > 3 ) break;
544  
545 <      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
546 <        //      if(ctrl.debug) {
545 >      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
546 >                 ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
547 >        if(ctrl.debug) {
548            cout << "failingLeptons : " << failingLeptons.size() << endl;
549            for( int f=0; f<failingLeptons.size(); f++ ) {
550              SimpleLepton  l = failingLeptons[f];
# Line 548 | Line 553 | int main(int argc, char** argv)
553                   << "pT: " << l.vec.Pt()  << "\t"
554                   << "eta: " << l.vec.Eta()
555                   << endl;
556 <            //    }
556 >          }
557          }
558          foTree.Fill();
559 +      } else {
560 +        if(ctrl.debug)
561 +          cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
562        }
563      }  
564  
# Line 559 | Line 567 | int main(int argc, char** argv)
567    foTree.getTree()->Write();
568    nt.WriteClose();
569  
570 <  cout << endl;
571 <  cout << endl;
570 >  if(ctrl.debug)
571 >    cout << endl << endl;
572    for( int i=0; i<cutvec.size(); i++ ) {
573      cout << "cut: " << i << "\t"
574           << "pass: " << cutvec[i];
# Line 581 | Line 589 | int main(int argc, char** argv)
589  
590   }
591  
584 //----------------------------------------------------------------------------
585 //
586 // Get primary vertices
587 // Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
588 // NOTE: if no PV is found from fitting tracks, the beamspot is used
589 //
590 //----------------------------------------------------------------------------
591 bool setPV(ControlFlags ctrl,
592           const mithep::Array<mithep::Vertex> * vtxArr,
593           mithep::Vertex & vtx)
594 //----------------------------------------------------------------------------
595 {
596  const Vertex *bestPV = 0;
597  float best_sumpt=-1;
598
599  // good PV requirements
600  const UInt_t   fMinNTracksFit = 0;
601  const Double_t fMinNdof       = 4;
602  const Double_t fMaxAbsZ       = 24;
603  const Double_t fMaxRho        = 2;
604  
605  for(int i=0; i<vtxArr->GetEntries(); ++i) {
606    const Vertex *pv = (Vertex*)(vtxArr->At(i));
607    if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
608    
609    // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
610    if(!pv->IsValid())                                continue;
611    if(pv->NTracksFit()       < fMinNTracksFit)       continue;
612    if(pv->Ndof()                 < fMinNdof)         continue;
613    if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
614    if(pv->Position().Rho()   > fMaxRho)              continue;    
615    
616    // take the first ...
617    bestPV = pv;
618    break;
619
620    // this never reached ...    
621    float tmp_sumpt=0;
622    for( int t=0; t<pv->NTracks(); t++ )
623      tmp_sumpt += pv->Trk(t)->Pt();
624    
625    if( tmp_sumpt > best_sumpt  ) {
626      bestPV = pv;
627      best_sumpt = tmp_sumpt;
628      if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
629    }
630  }
631  // sync
632  if(!bestPV) bestPV = vtxArr->At(0);
633  //if(!bestPV) return false;
634  vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
635  vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
636  return true;
637 };
638
639
592  
593   //----------------------------------------------------------------------------
594 < void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
643 < //----------------------------------------------------------------------------
644 < {
645 <  vector<SimpleLepton> lepvec = evtdat.Z1leptons;
646 <  lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
647 <  double w_offline=-1, werr_offline=0;
648 <  double w_online=-1, werr_online=0;
649 <  
650 <  vector< pair <double,double> > wlegs; // pair here is eff & err
651 <  vector< pair <float,float> > mvec;  // pair here is eta & pt
652 <
653 <  //      vector< pair <float,float> > evec;  // pair here is eta & pt
654 <  // now deal with medium vs loose
655 <  vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt
656 <  
657 <  for( int k=0; k<lepvec.size(); k++ ) {
658 <    if( abs(lepvec[k].type) == 13 ) {
659 <      mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
660 <      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
661 <                                                          lepvec[k].vec.Pt() ) );
662 <    } else {
663 <      
664 <      // now deal with medium vs loose
665 <      //              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
666 <      
667 <      std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
668 <      evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
669 <      
670 <      //              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
671 <      //                                                                 lepvec[k].vec.Pt() ) );
672 <      
673 <      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
674 <                                                          lepvec[k].vec.Pt(),
675 <                                                          lepvec[k].isTight ) );
676 <      
677 <    }
678 <  }
679 <  
680 <  pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
681 <  weights.woff    = offpair.first;
682 <  weights.werroff = offpair.second;
683 <  
684 <  pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
685 <  weights.won    = onpair.first;
686 <  weights.werron = onpair.second;
687 < }
688 <
689 <
690 < //----------------------------------------------------------------------------
691 < void initRunLumiRangeMap()
594 > void initRunLumiRangeMap(ControlFlags &ctrl)
595   //----------------------------------------------------------------------------
596   {
597    /*
# Line 700 | Line 603 | void initRunLumiRangeMap()
603    // r11b
604    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
605    */
703 };
704
705
706 //----------------------------------------------------------------------------
707 void initPUWeights()
708 //----------------------------------------------------------------------------
709 {
710  TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
711  hpu = (TH1D*)(puf->Get("puWeights"));
712 }
713
714 //----------------------------------------------------------------------------
715 double getPUWeight(unsigned npu)
716 //----------------------------------------------------------------------------
717 {
718  return hpu->GetBinContent(hpu->FindBin(npu));
719 }
720
721 //----------------------------------------------------------------------------
722 void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
723 //----------------------------------------------------------------------------
724 {
725  unsigned evt;
726  float rho;
727
728  cout << "initialzing EvtRhoMap ";
729  //FILE * f = fopen("evtrho.dat", "r");
730  //  FILE * f = fopen("allrho.cali.uniq.dat", "r");
731  FILE * f = fopen("allrhoAA.cali.uniq.dat", "r");
732  int lcount=0;
733  while( fscanf( f, "%u:%f", &evt, &rho ) ) {
734    if( feof(f) ) break;
735    evtrhoMap[evt] = rho;
736    lcount++;
737    if( !(lcount%1000) ) cout << "."; flush(cout);
738  }
739  cout << " done" << endl;
740 }
741
742 //----------------------------------------------------------------------------
743 unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
744                     vector<bool> &pfNoPileUpflag,
745                     mithep::Array<Vertex>      * vtxArr )
746 //----------------------------------------------------------------------------
747 {
748  for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
749    const PFCandidate *pf = fPFCandidates->At(i);
750    assert(pf);
751
752    if(pf->PFType() == PFCandidate::eHadron) {
753      if(pf->HasTrackerTrk() &&
754         vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
755         vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
756
757        pfNoPileUpflag.push_back(1);
758
759      } else {
760
761        Bool_t vertexFound = kFALSE;
762        const Vertex *closestVtx = 0;
763        Double_t dzmin = 10000;
764        
765        // loop over vertices
766        for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
767          const Vertex *vtx = vtxArr->At(j);
768          assert(vtx);
769          
770          if(pf->HasTrackerTrk() &&
771             vtx->HasTrack(pf->TrackerTrk()) &&
772             vtx->TrackWeight(pf->TrackerTrk()) > 0) {
773            vertexFound = kTRUE;
774            closestVtx = vtx;
775            break;
776          }
777          Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
778          if(dz < dzmin) {
779            closestVtx = vtx;
780            dzmin = dz;
781          }
782        }
783
784        //      if(fCheckClosestZVertex) {
785        if(1) {
786          // Fallback: if track is not associated with any vertex,
787          // associate it with the vertex closest in z
788          if(vertexFound || closestVtx != vtxArr->At(0)) {
789            //      pfPileUp->Add(pf);
790            pfNoPileUpflag.push_back(0);
791          } else {
792            pfNoPileUpflag.push_back(1);
793          }
794        } else {
795          if(vertexFound && closestVtx != vtxArr->At(0)) {
796            //      pfPileUp->Add(pf);
797            pfNoPileUpflag.push_back(0);
798          } else {
799            //      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
800            pfNoPileUpflag.push_back(1);
801          }
802        }
803      } //hadron & trk stuff
804    } else { // hadron
805      //      PFCandidate * ptr = pfNoPileUp->AddNew();
806      pfNoPileUpflag.push_back(1);
807    }
808  } // Loop over PF candidates
809
810  return pfNoPileUpflag.size();
811 }
606  
607 +  // 2012 only for now ...
608 +  std::string jsonFile;
609 +  if(ctrl.jsonFile=="")
610 +    jsonFile = "./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt";
611 +  else
612 +    jsonFile = ctrl.jsonFile;
613 +  rlrm.AddJSONFile(jsonFile);  
614 + };

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines