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.13 by khahn, Sat May 5 11:04:01 2012 UTC vs.
Revision 1.21 by khahn, Sat May 12 03:00:15 2012 UTC

# Line 5 | Line 5
5   #include <iostream>                 // standard I/O
6   #include <iomanip>                  // functions to format standard I/O
7   #include <bitset>  
8 + #include <map>  
9   using namespace std;
10  
11   //
# Line 34 | Line 35 | using namespace std;
35   #include "Vertex.h"
36   #include "PFCandidate.h"
37   #include "PFCandidateCol.h"
38 + #include "PileupInfoCol.h"
39 + #include "PileupEnergyDensity.h"
40   #include "TriggerMask.h"
41   #include "TriggerTable.h"
42   #include "Names.h"
# Line 46 | Line 49 | using namespace std;
49   #include "MuonSelection.h"
50   #include "ElectronSelection.h"
51   #include "IsolationSelection.h"
52 < #include "Selection.h"
52 >
53 > //#include "Selection.h"
54 > //#include "ReferenceSelection.h"
55 > #include "ReferenceSelection2.h"
56 >
57   #include "TriggerUtils.h"
58   #include "PassHLT.h"
59   #include "Angles.h"
# Line 54 | Line 61 | using namespace std;
61   #include "InfoStruct.h"
62   //#include "GenInfoStruct.h"
63   #include "WeightStruct.h"
64 < #include "GlobalDataAndFuncs.h"
64 > //#include "GlobalDataAndFuncs.h"
65 > #include "RunLumiRangeMap.h"
66 >
67   #include "AngleTuple.h"
68 + #include "FOTuple.h"
69  
70   #include "SampleWeight.h"
71   #include "EfficiencyWeightsInterface.h"
72  
73 + #include "SimpleLepton.h"
74 +
75   #ifndef CMSSW_BASE
76   #define CMSSW_BASE "../"
77   #endif
78  
79   using namespace mithep;
80  
81 + //
82 + // globals
83 + //----------------------------------------------------------------------------
84 + TH1D * hpu;
85 + RunLumiRangeMap rlrm;
86 + vector<SimpleLepton> failingLeptons ; // for fake estimation
87 + vector<SimpleLepton> passingLeptons;      // for fake estimation
88 + vector<unsigned> cutvec;
89 + vector<vector<unsigned> > zcutvec;
90 + vector<vector<unsigned> > zzcutvec;
91 + map<unsigned,float>       evtrhoMap;
92 + vector<string> cutstrs;
93 + bool passes_HLT_MC;
94  
95 + //
96 + // prototypes
97   //----------------------------------------------------------------------------
98   bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
99 + void initPUWeights();
100 + double getPUWeight(unsigned npu);
101 + void setEffiencyWeights(EventData&, WeightStruct& );
102 + void initRunLumiRangeMap();
103 + void initEvtRhoMap(map<unsigned,float> &);
104   //----------------------------------------------------------------------------
105  
106 < //=== MAIN =================================================================================================
106 >
107 > //
108 > // MAIN
109 > //----------------------------------------------------------------------------
110   int main(int argc, char** argv)
111 + //----------------------------------------------------------------------------
112   {
113 +  cutstrs.push_back(string("skim0"));              //0
114 +  cutstrs.push_back(string("skim1"));              //1
115 +  cutstrs.push_back(string("trigger"));            //2
116 +  // -------------------------------------------------
117 +  cutstrs.push_back(string("Z candidate"));        //3
118 +  cutstrs.push_back(string("good Z1"));            //4
119 +  // -------------------------------------------------
120 +  cutstrs.push_back(string("4l"));                 //5  
121 +  cutstrs.push_back(string("ZZ candidate"));       //6
122 +  cutstrs.push_back(string("good Z2"));            //7
123 +  cutstrs.push_back(string("ZZ 20/10"));           //8
124 +  cutstrs.push_back(string("resonance"));          //9
125 +  cutstrs.push_back(string("m4l>70"));             //10
126 +  cutstrs.push_back(string("m4l>100"));            //11
127 +
128 +
129 +
130 +
131    string cmsswpath(CMSSW_BASE);
132    cmsswpath.append("/src");
133    std::bitset<1024> triggerBits;
134    vector<vector<string> > inputFiles;
135    inputFiles.push_back(vector<string>());
136 <  map<string,double> entrymap; // number of unskimmed entries in each file
136 >  map<string,unsigned> entrymap; // number of unskimmed entries in each file
137 >  for( int i=0; i<cutstrs.size(); i++ ) cutvec.push_back(0);
138 >  for( int i=0; i<2; i++ )  {
139 >    zcutvec.push_back(vector<unsigned>());
140 >    for( int j=0; j<cutstrs.size(); j++ ) zcutvec[i].push_back(0);
141 >  }
142 >  for( int i=0; i<3; i++ )  {
143 >    zzcutvec.push_back(vector<unsigned>());
144 >    for( int j=0; j<cutstrs.size(); j++ ) zzcutvec[i].push_back(0);
145 >  }
146 >
147  
148    //
149    // args
# Line 96 | Line 160 | int main(int argc, char** argv)
160        return 1;
161      }
162    ctrl.dump();
99  assert( ctrl.mH != 0 );
163  
164  
165  
# Line 113 | Line 176 | int main(int argc, char** argv)
176        if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
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        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;
189      chain->AddFile(ctrl.inputfile.c_str());
190      hltchain->AddFile(ctrl.inputfile.c_str());
191    }
# Line 149 | Line 216 | int main(int argc, char** argv)
216    nt.makeKinematicsBranch(kinematics);
217    nt.makeInfoBranch(evtinfo);
218  
219 +  FOTuple foTree( nt.getFile(), (const char*)"FOtree");
220 +  foTree.makeFailedLeptonBranch(failingLeptons);
221 +  foTree.makeZLeptonBranch(passingLeptons);
222 +
223    TH1F * h_w_hpt;
224    if(ctrl.mc) {
225      nt.makeWeightBranch(weights);
226      //    nt.makeGenInfoBranch(geninfo);
227      initEfficiencyWeights();
228 +    initPUWeights();
229      /*
230      // Higgs only, pt reweighting
231      if( ctrl.mH <= 140 ) {
# Line 170 | Line 242 | int main(int argc, char** argv)
242      initRunLumiRangeMap();
243    }
244  
245 <  initMuonIDMVA();
245 >  //  initMuonIDMVA();
246    initMuonIsoMVA();
247    initElectronIDMVA();
248    initElectronIsoMVA();
249    initTrigger();
250 <  
251 <  
252 <  //##########################################################################
253 <  // Setup tree I/O
182 <  //##########################################################################
250 >
251 >  // hack
252 >  initEvtRhoMap(evtrhoMap);
253 >
254    
255    //
256 <  // Access samples and fill histograms
256 >  // Setup tree I/O
257 >  //--------------------------------------------------------------------------------------------------------------
258    TFile *inputFile=0;
259    TTree *eventTree=0;  
260    string currentFile("");
261  
190  // Data structures to store info from TTrees
262    mithep::EventHeader *info    = new mithep::EventHeader();
263    //  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
264    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
265    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
266    mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
267    mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
268 <  mithep::Array<mithep::PileupEnergyDensity>  *puArr         = new mithep::Array<mithep::PileupEnergyDensity>();
268 >  mithep::Array<mithep::PileupInfo>           *puArr         = new mithep::Array<mithep::PileupInfo>();
269 >  mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
270    mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
271    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
272    mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
# Line 205 | Line 277 | int main(int argc, char** argv)
277    TString fInfoName(Names::gkEvtHeaderBrn);
278    TString fPrimVtxName(Names::gkPVBrn);
279    TString fPFCandidateName(Names::gkPFCandidatesBrn);
280 +  TString fPileupInfoName(Names::gkPileupInfoBrn);
281    TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
282    TString fTrackName(Names::gkTrackBrn);
283    TString fTriggerMaskName(Names::gkHltBitBrn);
# Line 215 | Line 288 | int main(int argc, char** argv)
288    chain->SetBranchAddress(fMuonName,        &muonArr);
289    chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
290    chain->SetBranchAddress(fPFCandidateName, &pfArr);
291 <  chain->SetBranchAddress(fPileupEnergyDensityName, &puArr);
291 >  if( ctrl.mc ) chain->SetBranchAddress(fPileupInfoName,  &puArr);
292 >  chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
293    chain->SetBranchAddress(fTrackName, &trkArr);
294    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
295    cout << "hlttable: " << fTriggerTableName << endl;
# Line 240 | Line 314 | int main(int argc, char** argv)
314    cout << "nEntries: " << imax << endl;
315  
316  
317 <  //##########################################################################
318 <  // Loop !!!!!!!!! should alternate events here and +1 in the training ...
319 <  //##########################################################################
317 >  //
318 >  // Loop !!!!!!!!!
319 >  //--------------------------------------------------------------------------------------------------------------
320    for(UInt_t ientry=0; ientry<imax; ientry++)
321      {
322        chain->GetEntry(ientry);
323        if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
250
324        setPV( ctrl, vtxArr, vtx);
325  
326 +      cout << "origrho: " << puDArr->At(0)->Rho()
327 +           << "\torigrholoweta: " << puDArr->At(0)->RhoLowEta()
328 +           << endl;
329 +
330 +
331 +      //
332 +      // data/MC
333 +      //
334        if(ctrl.mc) {
335 <        weights.w = getWeight(xstab,entrymap,chain);
335 >        //
336 >        // xsec & PU weights
337 >        //
338 >        //      weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
339 >        int npu = -1;
340 >        for(int i=0;i<puArr->GetEntries();i++) {
341 >          if(puArr->At(i)->GetBunchCrossing() == 0)
342 >            npu = puArr->At(i)->GetPU_NumInteractions();
343 >        }
344 >        assert(npu>=0);
345 >        evtinfo.npu;
346 >        weights.npuw = getPUWeight(evtinfo.npu);
347 >        //      cout << "weight: " << weights.w << endl;
348 >        //
349 >        // trigger
350 >        //
351 >        if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
352 >          currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
353 >          hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
354 >          hltchain->GetEntry(0);
355 >          hltTable->Clear();
356 >          fillTriggerBits(hltTable, hltTableStrings );
357 >        }
358 >        if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
359 >        fillTriggerBits( hltTable, trigMask, triggerBits );
360 >        passes_HLT_MC = true; // passed to apply as extern global, just for sync
361 >        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) {
362 >          passes_HLT_MC = false;
363 >        }
364        } else {
365 +        //
366 +        // JSON
367 +        //
368 +        /*
369 +        RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
370 +        if( !(rlrm.HasRunLumi(rl)) )  {
371 +          if( ctrl.debug ) cout << "\tfails JSON" << endl;
372 +          continue;
373 +        }
374 +        */
375 +        //
376 +        // trigger
377 +        //
378          if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
379            currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
380            hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
# Line 261 | Line 383 | int main(int argc, char** argv)
383            fillTriggerBits(hltTable, hltTableStrings );
384          }
385          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
386 +        /*
387          fillTriggerBits( hltTable, trigMask, triggerBits );
265      }
266      
267      
268      //
269      // trigger
270      //
271      if( !ctrl.mc ) {
388          if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
389            if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
390            continue;
391          }
392 +        */
393        }
394 <      
394 >
395        //
396        // lepton/kinematic selection ...
397        //
398 +      failingLeptons.clear();
399 +      passingLeptons.clear();
400        EventData ret4l =
401 <        apply_HZZ4L_selection(ctrl, info,
401 >
402 >        /*
403 >        apply_HZZ4L_reference2_selection(ctrl, info,
404 >                                         vtx,
405 >                                         pfArr,
406 >                                         puDArr,
407 >                                         electronArr,
408 >                                         &electronReferencePreSelection,
409 >                                         &electronReferenceIDMVASelection,
410 >                                         &electronReferenceIsoSelection,
411 >                                         muonArr,
412 >                                         &muonReferencePreSelection,
413 >                                         &muonIDPFSelection,
414 >                                         &muonReferenceIsoSelection);
415 >        */
416 >        apply_HZZ4L_reference2_selection(ctrl, info,
417                                vtx,
418                                pfArr,
419 <                              puArr,
419 >                              puDArr,
420                                electronArr,
421 <                              &electronPreSelection,
422 <                              //                              &electronBDTSelection,
289 <                              &electronIDMVASelection,
290 <                              //                              &electronIsoSelection,
421 >                              &electronReferencePreSelection,
422 >                              &electronReferenceIDMVASelection,
423                                &electronIsoMVASelection,
424                                muonArr,
425 <                              &muonPreSelection,
425 >                              &muonReferencePreSelection,
426                                &muonIDPFSelection,
295                              //                              &muonIDMVASelection,
296                              //                              &passMuonSelection,
427                                &muonIsoMVASelection);
428 +
429 +
430 +
431        if( ctrl.debug ) cout << endl;
432 +
433 +
434 +      int Z1type=0, Z2type=0;
435 +      if( ret4l.status.selectionBits.to_ulong() >= PASS_ZCANDIDATE ) {
436 +        if( abs(ret4l.Z1leptons[0].type) == 11  ) Z1type = 11;
437 +        if( abs(ret4l.Z1leptons[0].type) == 13  ) Z1type = 13;
438 +      }  
439 +      if( ret4l.status.selectionBits.to_ulong() >= PASS_ZZCANDIDATE ) {
440 +        if( abs(ret4l.Z2leptons[0].type) == 11  ) Z2type = 11;
441 +        if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
442 +      }  
443 +      
444 +      cout  << "bits: " << ret4l.status.selectionBits  << "\t"
445 +            << "Z1: " << Z1type << "\t"
446 +            << "Z2: " << Z2type << endl;
447 +      cout << endl;
448 +                      
449        
450        if( ret4l.status.pass() ) {
301        
451          fillAngles(ret4l,angles);
452          fillKinematics(ret4l,kinematics);
453          fillEventInfo(info,evtinfo);
454 <        //if(ctrl.mc) fillGenInfo(ginfo,geninfo);
454 >        if( ctrl.mc) {
455 >        // fillGenInfo(ginfo,geninfo);
456 >          setEffiencyWeights(ret4l, weights);
457 >           if(ctrl.debug)
458 >             cout << "w: " << weights.w << "\t"
459 >                  << "won: " << weights.won << "\t"
460 >                  << "woff: " << weights.woff << endl;
461 >        }
462          
463          /*
464          // only for Higgs < 140
# Line 323 | Line 479 | int main(int argc, char** argv)
479          pass++;
480          //      if( pass > 3 ) break;
481  
482 <      }
482 >      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
483 >        //      if(ctrl.debug) {
484 >          cout << "failingLeptons : " << failingLeptons.size() << endl;
485 >          for( int f=0; f<failingLeptons.size(); f++ ) {
486 >            SimpleLepton  l = failingLeptons[f];
487 >            cout << "f: " << f << "\t"
488 >                 << "type: " << l.type << "\t"
489 >                 << "pT: " << l.vec.Pt()  << "\t"
490 >                 << "eta: " << l.vec.Eta()
491 >                 << endl;
492 >            //    }
493 >        }
494 >        foTree.Fill();
495 >      }
496      }  
497  
498 +  
499 +  foTree.getFile()->cd();
500 +  foTree.getTree()->Write();
501    nt.WriteClose();
502 +
503 +  cout << endl;
504 +  cout << endl;
505 +  for( int i=0; i<cutvec.size(); i++ ) {
506 +    cout << "cut: " << i << "\t"
507 +         << "pass: " << cutvec[i];
508 +
509 +    if( i>PASS_TRIGGER&&i<=PASS_GOODZ1 )
510 +      cout << "\t2e: " << zcutvec[0][i]
511 +           << "\t2m: " << zcutvec[1][i];
512 +
513 +    if( i>PASS_ZCANDIDATE )
514 +      cout << "\t4e: " << zzcutvec[0][i]
515 +           << "\t4m: " << zzcutvec[1][i]
516 +           << "\t2e2m: " << zzcutvec[2][i];
517 +
518 +    cout << "\t" << cutstrs[i] << endl;
519 +
520 +    cout << endl;
521 +  }
522 +
523   }
524  
525   //----------------------------------------------------------------------------
# Line 382 | Line 575 | bool setPV(ControlFlags ctrl,
575    vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
576    return true;
577   };
578 +
579 +
580 +
581 + //----------------------------------------------------------------------------
582 + void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
583 + //----------------------------------------------------------------------------
584 + {
585 +  vector<SimpleLepton> lepvec = evtdat.Z1leptons;
586 +  lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
587 +  double w_offline=-1, werr_offline=0;
588 +  double w_online=-1, werr_online=0;
589 +  
590 +  vector< pair <double,double> > wlegs; // pair here is eff & err
591 +  vector< pair <float,float> > mvec;  // pair here is eta & pt
592 +
593 +  //      vector< pair <float,float> > evec;  // pair here is eta & pt
594 +  // now deal with medium vs loose
595 +  vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt
596 +  
597 +  for( int k=0; k<lepvec.size(); k++ ) {
598 +    if( abs(lepvec[k].type) == 13 ) {
599 +      mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
600 +      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
601 +                                                          lepvec[k].vec.Pt() ) );
602 +    } else {
603 +      
604 +      // now deal with medium vs loose
605 +      //              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
606 +      
607 +      std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
608 +      evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
609 +      
610 +      //              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
611 +      //                                                                 lepvec[k].vec.Pt() ) );
612 +      
613 +      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
614 +                                                          lepvec[k].vec.Pt(),
615 +                                                          lepvec[k].isTight ) );
616 +      
617 +    }
618 +  }
619 +  
620 +  pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
621 +  weights.woff    = offpair.first;
622 +  weights.werroff = offpair.second;
623 +  
624 +  pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
625 +  weights.won    = onpair.first;
626 +  weights.werron = onpair.second;
627 + }
628 +
629 +
630 + //----------------------------------------------------------------------------
631 + void initRunLumiRangeMap()
632 + //----------------------------------------------------------------------------
633 + {
634 +  /*
635 +  rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
636 +  //  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
637 +  rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
638 +  rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
639 +  rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
640 +  // r11b
641 +  rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));  
642 +  */
643 + };
644 +
645 +
646 + //----------------------------------------------------------------------------
647 + void initPUWeights()
648 + //----------------------------------------------------------------------------
649 + {
650 +  TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
651 +  hpu = (TH1D*)(puf->Get("puWeights"));
652 + }
653 +
654 + //----------------------------------------------------------------------------
655 + double getPUWeight(unsigned npu)
656 + //----------------------------------------------------------------------------
657 + {
658 +  return hpu->GetBinContent(hpu->FindBin(npu));
659 + }
660 +
661 + //----------------------------------------------------------------------------
662 + void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
663 + //----------------------------------------------------------------------------
664 + {
665 +  unsigned evt;
666 +  float rho;
667 +
668 +  FILE * f = fopen("evtrho.dat", "r");
669 +  while( fscanf( f, "%u:%f", &evt, &rho ) ) {
670 +    evtrhoMap[evt] = rho;
671 +  }
672 +  
673 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines