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.44 by dkralph, Tue Jun 19 08:35:51 2012 UTC

# Line 40 | Line 40 | using namespace std;
40   #include "MCParticle.h"
41   #include "TriggerMask.h"
42   #include "TriggerTable.h"
43 + #include "TriggerObject.h"
44 + #include "TriggerObjectRel.h"
45   #include "Names.h"
46   #include "BaseMod.h"
47 + #include "TriggerObjectsTable.h"
48  
49   //
50   // our headers
# Line 50 | Line 53 | using namespace std;
53   #include "MuonSelection.h"
54   #include "ElectronSelection.h"
55   #include "IsolationSelection.h"
56 <
54 < //#include "Selection.h"
55 < //#include "ReferenceSelection.h"
56 < #include "ReferenceSelection2.h"
56 > #include "ReferenceSelection.h"
57  
58   #include "TriggerUtils.h"
59   #include "PassHLT.h"
60   #include "Angles.h"
61   #include "KinematicsStruct.h"
62   #include "InfoStruct.h"
63 < //#include "GenInfoStruct.h"
63 > #include "GenInfoStruct.h"
64   #include "WeightStruct.h"
65 //#include "GlobalDataAndFuncs.h"
66 #include "RunLumiRangeMap.h"
67
65   #include "AngleTuple.h"
66   #include "FOTuple.h"
67  
68 + #include "RunLumiRangeMap.h"
69   #include "SampleWeight.h"
70   #include "EfficiencyWeightsInterface.h"
71 + #include "SelectionFuncs.h"
72  
73   #include "SimpleLepton.h"
74  
# Line 97 | Line 96 | vector<bool>   PFnoPUflag;;
96   //
97   // prototypes
98   //----------------------------------------------------------------------------
100 bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
99   void initPUWeights();
100   double getPUWeight(unsigned npu);
103 void setEffiencyWeights(EventData&, WeightStruct& );
101   void initRunLumiRangeMap();
105 void initEvtRhoMap(map<unsigned,float> &);
106 unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
107                         vector<bool> &pfNoPileUpFlag,
108                         mithep::Array<Vertex>      * vtxArr );
102   //----------------------------------------------------------------------------
103  
104  
# Line 117 | Line 110 | int main(int argc, char** argv)
110   {
111    cutstrs.push_back(string("skim0"));              //0
112    cutstrs.push_back(string("skim1"));              //1
113 +  cutstrs.push_back(string("skim2"));              //1
114    cutstrs.push_back(string("trigger"));            //2
115    // -------------------------------------------------
116    cutstrs.push_back(string("Z candidate"));        //3
# Line 128 | Line 122 | int main(int argc, char** argv)
122    cutstrs.push_back(string("ZZ 20/10"));           //8
123    cutstrs.push_back(string("resonance"));          //9
124    cutstrs.push_back(string("m4l>70"));             //10
125 <  cutstrs.push_back(string("m4l>100"));            //11
125 >  cutstrs.push_back(string("m4l>100, mZ2 > 12"));                  //11
126  
127  
128  
# Line 175 | Line 169 | int main(int argc, char** argv)
169    TChain * hltchain = new TChain("HLT");
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 182 | 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      }
184    } else {
185      cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
186 <    unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
187 <    cout << "tmpent: " << tmpent << endl;
188 <    entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
189 <    cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
186 >    unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
187 >    entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
188 >    cout << "unskimmed entries: " << unsk_ents << endl;
189 >    total_unskimmed += unsk_ents;
190      chain->AddFile(ctrl.inputfile.c_str());
191      hltchain->AddFile(ctrl.inputfile.c_str());
192    }
193 +  // // write the total number of unskimmed events that went into making this output file to a text file
194 +  // writeEntries(ctrl,total_unskimmed);
195  
196    const char * ofname;
197    if( strcmp( ctrl.outputfile.c_str(), "") ) {
# Line 207 | Line 205 | int main(int argc, char** argv)
205    SimpleTable xstab(xspath.c_str());
206  
207  
208 +
209    //
210    // Setup
211    //--------------------------------------------------------------------------------------------------------------
# Line 214 | Line 213 | int main(int argc, char** argv)
213    KinematicsStruct kinematics;
214    InfoStruct evtinfo;
215    WeightStruct weights;
216 <  //  GenInfoStruct geninfo;
216 >  GenInfoStruct geninfo;
217  
218    AngleTuple nt( (const char*)ofname, (const char*)"zznt");
219    nt.makeAngleBranch(angles);
220    nt.makeKinematicsBranch(kinematics);
221    nt.makeInfoBranch(evtinfo);
222  
223 +  
224    FOTuple foTree( nt.getFile(), (const char*)"FOtree");
225    foTree.makeFailedLeptonBranch(failingLeptons);
226 <  foTree.makeZLeptonBranch(passingLeptons);
226 >  foTree.makePassedLeptonBranch(passingLeptons);
227  
228 <  TH1F * h_w_hpt;
228 >  TH1F * h_wf11_hpt;
229    if(ctrl.mc) {
230      nt.makeWeightBranch(weights);
231 <    //    nt.makeGenInfoBranch(geninfo);
232 <    initEfficiencyWeights();
231 >    if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
232 >    if(ctrl.efftype != string("WTF?")) initEfficiencyWeights();
233      initPUWeights();
234 <    /*
234 >
235      // Higgs only, pt reweighting
236      if( ctrl.mH <= 140 ) {
237        char wptname[256];
# Line 239 | Line 239 | int main(int argc, char** argv)
239        TFile * f_hpt = new TFile(wptname);
240        f_hpt->Print();
241        sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
242 <      h_w_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
243 <      h_w_hpt->Print();
242 >      h_wf11_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
243      }
244 <    */
244 >
245    } else {
246      initRunLumiRangeMap();
247    }
248  
249    //  initMuonIDMVA();
250 <  initMuonIsoMVA();
252 <  initElectronIDMVA();
253 <  initElectronIsoMVA();
250 >  initElectronIDMVAV1();
251    initTrigger();
252  
256  // hack
257  initEvtRhoMap(evtrhoMap);
258
259  
253    //
254    // Setup tree I/O
255    //--------------------------------------------------------------------------------------------------------------
# Line 264 | Line 257 | int main(int argc, char** argv)
257    TTree *eventTree=0;  
258    string currentFile("");
259  
260 +  UInt_t fNMaxTriggers = 1024;
261    mithep::EventHeader *info    = new mithep::EventHeader();
268  //  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
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);
# Line 287 | Line 288 | int main(int argc, char** argv)
288    TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
289    TString fTrackName(Names::gkTrackBrn);
290    TString fMCParticleName(Names::gkMCPartBrn);
291 +  TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
292    TString fTriggerMaskName(Names::gkHltBitBrn);
293    TString fTriggerTableName(Names::gkHltTableBrn);
294 +  TString fTriggerLabelName(Names::gkHltLabelBrn);
295 +  TString fTriggerObjectName(Names::gkHltObjBrn);
296 +  TString fTriggerObjectRelsName(Form("HLTObjectsRelation"));
297 +
298 + //   TString fHLTLabName(Names::gkHltLabelBrn);
299 + //   TString fHLTTabName(Names::gkHltTableBrn);
300 + //   TString fHLTTabNamePub(Form("%sFwk",fHLTTabName.Data()));
301 + //   TString fHLTLabNamePub(Form("%sFwk",fHLTLabName.Data()));
302 + //   TString fObjsNamePub(Form("%sFwk",fTriggerObjectName.Data()));
303 +
304 +
305    
306    chain->SetBranchAddress(fInfoName,        &info);
307    chain->SetBranchAddress(fElectronName,    &electronArr);
# Line 298 | Line 311 | int main(int argc, char** argv)
311    if( ctrl.mc ) {
312      chain->SetBranchAddress(fPileupInfoName,  &puArr);
313      chain->SetBranchAddress(fMCParticleName,  &mcArr);
314 +    chain->SetBranchAddress(fMCEvtInfoName,  &mcEvtInfo);
315    }
316    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
317    chain->SetBranchAddress(fTrackName, &trkArr);
318    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
319 +  chain->SetBranchAddress(fTriggerObjectName,  &hltObjArr);
320 +  chain->SetBranchAddress(fTriggerObjectRelsName,  &hltRelsArr);
321 +  //  chain->SetBranchAddress(fHLTTabNamePub,  &hltTab);
322 +  //  chain->SetBranchAddress(fHLTLabNamePub,  &hltLab);
323 +
324    cout << "hlttable: " << fTriggerTableName << endl;
325    hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
326 +  hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
327  
328    mithep::Vertex              vtx;          // best primary vertex in the event
329  
# Line 323 | Line 343 | int main(int argc, char** argv)
343    int imax = chain->GetEntries();
344    cout << "nEntries: " << imax << endl;
345  
326
346    //
347    // Loop !!!!!!!!!
348    //--------------------------------------------------------------------------------------------------------------
349    for(UInt_t ientry=0; ientry<imax; ientry++)
350      {
351        chain->GetEntry(ientry);
352 +      if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
353  
354 <      /*
355 <      bool found_tau=false;
356 <      if( ctrl.mc) {  // && ctrl.debug
357 <        cout << "nparticles: " << mcArr->GetEntries() << endl;
358 <        cout << "-----------------------------------" << endl;
359 <        for( int i=0; i<mcArr->GetEntries(); i++ ) {
360 <          const mithep::MCParticle *mc = mcArr->At(i);
361 <          cout << "i: " << i << "\t"
342 <               << "id: " << mc->AbsPdgId()
343 <               << endl;
344 <          if( mc->Is(MCParticle::kTau) ) {
345 <            found_tau=true;
346 <          }
347 <        }
354 >      if( ctrl.debug ) {
355 >        cout << "-----------------------------------------------------------------" << endl;
356 >        cout << "-----------------------------------------------------------------" << endl;
357 >        cout << "Run: " << info->RunNum()
358 >             << "\tEvt: " << info->EvtNum()
359 >             << "\tLumi: " << info->LumiSec()
360 >             << endl;
361 >        cout << "-----------------------------------------------------------------" << endl;
362        }
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      */
363  
367      if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
368      setPV( ctrl, vtxArr, vtx);
369      //if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
364  
365 <      float rho = evtrhoMap[info->EvtNum()];
365 >      string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
366 >      setEra( fname, ctrl );
367 >
368 >
369  
370        //
371        // pfNoPU
# Line 383 | Line 380 | int main(int argc, char** argv)
380        //
381        if(ctrl.mc) {
382          //
383 +        // gen info
384 +        //
385 +        if( ctrl.fillGen )
386 +          fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
387 +
388 +        //
389          // xsec & PU weights
390          //
391          weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
# Line 395 | Line 398 | int main(int argc, char** argv)
398          assert(npu>=0);
399          evtinfo.npu;
400          weights.npuw = getPUWeight(evtinfo.npu);
401 <        //      cout << "weight: " << weights.w << 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 >
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 +        printTriggerObjs( hltTable, fTrigObjs );
459 +
460          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
461          /*
462          fillTriggerBits( hltTable, trigMask, triggerBits );
# Line 452 | Line 476 | int main(int argc, char** argv)
476  
477  
478          // reference
479 <        /*
480 <        apply_HZZ4L_reference2_selection(ctrl, info,
481 <                                         vtx,
482 <                                         pfArr,
483 <                                         //rho,
484 <                                         puDArr,
485 <                                         electronArr,
486 <                                         &electronReferencePreSelection,
487 <                                         &electronReferenceIDMVASelection,
488 <                                         &electronReferenceIsoSelection,
489 <                                         muonArr,
490 <                                         &muonReferencePreSelection,
491 <                                         &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 <
479 >        apply_HZZ4L_reference_selection(ctrl, info,
480 >                                        vtxArr,
481 >                                        pfArr,
482 >                                        puDArr,
483 >                                        electronArr,
484 >                                        &electronReferencePreSelection,
485 >                                        &electronReferenceIDMVASelectionV1,
486 >                                        &electronReferenceIsoSelection,
487 >                                        muonArr,
488 >                                        &muonReferencePreSelection,
489 >                                        &muonIDPFSelection,
490 >                                        &muonReferenceIsoSelection);
491 >      
492        if( ctrl.debug ) cout << endl;
493  
494  
# Line 509 | Line 512 | int main(int argc, char** argv)
512        if( ret4l.status.pass() ) {
513          fillAngles(ret4l,angles);
514          fillKinematics(ret4l,kinematics);
515 +        fillMassErrors(ret4l,muonArr,electronArr,kinematics);
516          fillEventInfo(info,evtinfo);
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(ret4l, weights);
520 >          } else {
521 >            weights.won = weights.woff = 1.;
522 >          }
523 >          if(ctrl.debug)
524 >            cout << "w: " << weights.w << "\t"
525 >                 << "won: " << weights.won << "\t"
526 >                 << "woff: " << weights.woff << endl;
527          }
528          
529 <        /*
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 <        */
529 >
530          nt.Fill();
531          
532          cerr  << "PASS:: "
# Line 538 | Line 541 | int main(int argc, char** argv)
541          pass++;
542          //      if( pass > 3 ) break;
543  
544 <      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
545 <        //      if(ctrl.debug) {
544 >      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
545 >                 ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
546 >        if(ctrl.debug) {
547            cout << "failingLeptons : " << failingLeptons.size() << endl;
548            for( int f=0; f<failingLeptons.size(); f++ ) {
549              SimpleLepton  l = failingLeptons[f];
# Line 548 | Line 552 | int main(int argc, char** argv)
552                   << "pT: " << l.vec.Pt()  << "\t"
553                   << "eta: " << l.vec.Eta()
554                   << endl;
555 <            //    }
555 >          }
556          }
557          foTree.Fill();
558 +      } else {
559 +        cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
560        }
561      }  
562  
# Line 581 | Line 587 | int main(int argc, char** argv)
587  
588   }
589  
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
640
641 //----------------------------------------------------------------------------
642 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
590  
591   //----------------------------------------------------------------------------
592   void initRunLumiRangeMap()
# Line 700 | Line 601 | void initRunLumiRangeMap()
601    // r11b
602    rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
603    */
604 +
605 +  // 2012 only for now ...
606 +  rlrm.AddJSONFile(std::string("./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt"));  
607 +
608   };
609  
610  
# Line 717 | Line 622 | double getPUWeight(unsigned npu)
622   {
623    return hpu->GetBinContent(hpu->FindBin(npu));
624   }
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 }
812

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines