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.16 by khahn, Wed May 9 14:57:20 2012 UTC vs.
Revision 1.27 by anlevin, Wed May 23 22:06:40 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 36 | Line 37 | using namespace std;
37   #include "PFCandidateCol.h"
38   #include "PileupInfoCol.h"
39   #include "PileupEnergyDensity.h"
40 + #include "MCParticle.h"
41   #include "TriggerMask.h"
42   #include "TriggerTable.h"
43   #include "Names.h"
# Line 50 | Line 52 | using namespace std;
52   #include "IsolationSelection.h"
53  
54   //#include "Selection.h"
55 < #include "ReferenceSelection.h"
55 > //#include "ReferenceSelection.h"
56 > #include "ReferenceSelection2.h"
57  
58   #include "TriggerUtils.h"
59   #include "PassHLT.h"
# Line 83 | Line 86 | TH1D * hpu;
86   RunLumiRangeMap rlrm;
87   vector<SimpleLepton> failingLeptons ; // for fake estimation
88   vector<SimpleLepton> passingLeptons;      // for fake estimation
89 + vector<unsigned> cutvec;
90 + vector<vector<unsigned> > zcutvec;
91 + vector<vector<unsigned> > zzcutvec;
92 + map<unsigned,float>       evtrhoMap;
93 + vector<string> cutstrs;
94 + bool passes_HLT_MC;
95 + vector<bool>   PFnoPUflag;;
96  
97   //
98   // prototypes
99   //----------------------------------------------------------------------------
90 bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
100   void initPUWeights();
101   double getPUWeight(unsigned npu);
102   void setEffiencyWeights(EventData&, WeightStruct& );
103   void initRunLumiRangeMap();
104 + void initEvtRhoMap(map<unsigned,float> &);
105 + unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
106 +                         vector<bool> &pfNoPileUpFlag,
107 +                         const mithep::VertexCol     * vtxArr );
108   //----------------------------------------------------------------------------
109  
110  
# Line 101 | Line 114 | void initRunLumiRangeMap();
114   int main(int argc, char** argv)
115   //----------------------------------------------------------------------------
116   {
117 +  cutstrs.push_back(string("skim0"));              //0
118 +  cutstrs.push_back(string("skim1"));              //1
119 +  cutstrs.push_back(string("skim2"));              //1
120 +  cutstrs.push_back(string("trigger"));            //2
121 +  // -------------------------------------------------
122 +  cutstrs.push_back(string("Z candidate"));        //3
123 +  cutstrs.push_back(string("good Z1"));            //4
124 +  // -------------------------------------------------
125 +  cutstrs.push_back(string("4l"));                 //5  
126 +  cutstrs.push_back(string("ZZ candidate"));       //6
127 +  cutstrs.push_back(string("good Z2"));            //7
128 +  cutstrs.push_back(string("ZZ 20/10"));           //8
129 +  cutstrs.push_back(string("resonance"));          //9
130 +  cutstrs.push_back(string("m4l>70"));             //10
131 +  cutstrs.push_back(string("m4l>100, mZ2 > 12"));                  //11
132 +
133 +
134 +
135 +
136    string cmsswpath(CMSSW_BASE);
137    cmsswpath.append("/src");
138    std::bitset<1024> triggerBits;
139    vector<vector<string> > inputFiles;
140    inputFiles.push_back(vector<string>());
141    map<string,unsigned> entrymap; // number of unskimmed entries in each file
142 +  for( int i=0; i<cutstrs.size(); i++ ) cutvec.push_back(0);
143 +  for( int i=0; i<2; i++ )  {
144 +    zcutvec.push_back(vector<unsigned>());
145 +    for( int j=0; j<cutstrs.size(); j++ ) zcutvec[i].push_back(0);
146 +  }
147 +  for( int i=0; i<3; i++ )  {
148 +    zzcutvec.push_back(vector<unsigned>());
149 +    for( int j=0; j<cutstrs.size(); j++ ) zzcutvec[i].push_back(0);
150 +  }
151 +
152  
153    //
154    // args
# Line 178 | Line 220 | int main(int argc, char** argv)
220    nt.makeAngleBranch(angles);
221    nt.makeKinematicsBranch(kinematics);
222    nt.makeInfoBranch(evtinfo);
223 <  FOTuple foTree( nt.getFile(), (const char*)"fo");
223 >
224 >  FOTuple foTree( nt.getFile(), (const char*)"FOtree");
225    foTree.makeFailedLeptonBranch(failingLeptons);
226    foTree.makeZLeptonBranch(passingLeptons);
227  
# Line 209 | Line 252 | int main(int argc, char** argv)
252    initElectronIDMVA();
253    initElectronIsoMVA();
254    initTrigger();
255 <  
255 >
256 >  // hack
257 >  initEvtRhoMap(evtrhoMap);
258  
259    
260    //
# Line 223 | Line 268 | int main(int argc, char** argv)
268    //  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
269    mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
270    mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
271 <  mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
271 >  const mithep::VertexCol              *vtxArr        = 0;//new mithep::Array<mithep::Vertex>();
272    mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
273    mithep::Array<mithep::PileupInfo>           *puArr         = new mithep::Array<mithep::PileupInfo>();
274    mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
275    mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
276 +  mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
277    mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
278    mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
279    vector<string>                              *hltTableStrings  = new vector<string>();
# Line 240 | Line 286 | int main(int argc, char** argv)
286    TString fPileupInfoName(Names::gkPileupInfoBrn);
287    TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
288    TString fTrackName(Names::gkTrackBrn);
289 +  TString fMCParticleName(Names::gkMCPartBrn);
290    TString fTriggerMaskName(Names::gkHltBitBrn);
291    TString fTriggerTableName(Names::gkHltTableBrn);
292    
# Line 248 | Line 295 | int main(int argc, char** argv)
295    chain->SetBranchAddress(fMuonName,        &muonArr);
296    chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
297    chain->SetBranchAddress(fPFCandidateName, &pfArr);
298 <  chain->SetBranchAddress(fPileupInfoName,  &puArr);
298 >  if( ctrl.mc ) {
299 >    chain->SetBranchAddress(fPileupInfoName,  &puArr);
300 >    chain->SetBranchAddress(fMCParticleName,  &mcArr);
301 >  }
302    chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
303    chain->SetBranchAddress(fTrackName, &trkArr);
304    chain->SetBranchAddress(fTriggerMaskName, &trigMask);
# Line 280 | Line 330 | int main(int argc, char** argv)
330    for(UInt_t ientry=0; ientry<imax; ientry++)
331      {
332        chain->GetEntry(ientry);
333 +
334 +      /*
335 +      bool found_tau=false;
336 +      if( ctrl.mc) {  // && ctrl.debug
337 +        cout << "nparticles: " << mcArr->GetEntries() << endl;
338 +        cout << "-----------------------------------" << endl;
339 +        for( int i=0; i<mcArr->GetEntries(); i++ ) {
340 +          const mithep::MCParticle *mc = mcArr->At(i);
341 +          cout << "i: " << i << "\t"
342 +               << "id: " << mc->AbsPdgId()
343 +               << endl;
344 +          if( mc->Is(MCParticle::kTau) ) {
345 +            found_tau=true;
346 +          }
347 +        }
348 +      }
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 +      */
366 +
367        if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
368 <      if( ientry>100 )break;
285 <      setPV( ctrl, vtxArr, vtx);
368 >      //if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
369  
370 +      float rho = evtrhoMap[info->EvtNum()];
371 +
372 +      //
373 +      // pfNoPU
374 +      //
375 +      //mithep::Array<PFCandidate> pfNoPileUp;
376 +      PFnoPUflag.clear();
377 +      int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
378 +      assert(pfnopu_size == pfArr->GetEntries());
379  
380        //
381        // data/MC
# Line 292 | Line 384 | int main(int argc, char** argv)
384          //
385          // xsec & PU weights
386          //
387 <        weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
387 >        weights.w = 1;//getWeight(xstab,entrymap,chain)/ctrl.totalMC;
388 >        cout << "xsec weight: " << weights.w << endl;
389          int npu = -1;
390          for(int i=0;i<puArr->GetEntries();i++) {
391            if(puArr->At(i)->GetBunchCrossing() == 0)
# Line 302 | Line 395 | int main(int argc, char** argv)
395          evtinfo.npu;
396          weights.npuw = getPUWeight(evtinfo.npu);
397          //      cout << "weight: " << weights.w << endl;
398 +        //
399 +        // trigger
400 +        //
401 +        if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
402 +          currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
403 +          hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
404 +          hltchain->GetEntry(0);
405 +          hltTable->Clear();
406 +          fillTriggerBits(hltTable, hltTableStrings );
407 +        }
408 +        if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
409 +        fillTriggerBits( hltTable, trigMask, triggerBits );
410 +        passes_HLT_MC = true; // passed to apply as extern global, just for sync
411 +        if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
412 +          passes_HLT_MC = false;
413 +        }
414        } else {
415          //
416          // JSON
417          //
418 +        /*
419          RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
420          if( !(rlrm.HasRunLumi(rl)) )  {
421            if( ctrl.debug ) cout << "\tfails JSON" << endl;
422            continue;
423          }
424 +        */
425          //
426          // trigger
427          //
# Line 322 | Line 433 | int main(int argc, char** argv)
433            fillTriggerBits(hltTable, hltTableStrings );
434          }
435          if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
436 +        /*
437          fillTriggerBits( hltTable, trigMask, triggerBits );
438          if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
439            if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
440            continue;
441          }
442 +        */
443        }
444  
445        //
# Line 335 | Line 448 | int main(int argc, char** argv)
448        failingLeptons.clear();
449        passingLeptons.clear();
450        EventData ret4l =
451 <        apply_HZZ4L_reference_selection(ctrl, info,
452 <                              vtx,
453 <                              pfArr,
454 <                              puDArr,
455 <                              electronArr,
456 <                              &electronReferencePreSelection,
457 <                              &electronReferenceIDMVASelection,
458 <                              &electronReferenceIsoSelection,
459 <                              muonArr,
460 <                              &muonReferencePreSelection,
461 <                              &muonIDPFSelection,
462 <                              &muonReferenceIsoSelection);
463 <        /*
464 <        apply_HZZ4L_reference_selection(ctrl, info,
465 <                              vtx,
466 <                              pfArr,
467 <                              puArr,
468 <                              electronArr,
469 <                              &electronPreSelection,
470 <                              &electronIDMVASelection,
471 <                              &electronIsoMVASelection,
472 <                              muonArr,
473 <                              &muonPreSelection,
474 <                              &muonIDPFSelection,
475 <                              &muonIsoMVASelection);
476 <        */
451 >
452 >
453 >        // reference
454 >        apply_HZZ4L_reference2_selection(ctrl, info,
455 >                                         vtxArr,
456 >                                         pfArr,
457 > #ifdef HACKED_RHOS
458 >                                         rho,
459 > #else
460 >                                         puDArr,
461 > #endif
462 >                                         electronArr,
463 >                                         &electronReferencePreSelection,
464 >                                         &electronReferenceIDMVASelectionV1,
465 >                                         &electronReferenceIsoSelection,
466 >                                         muonArr,
467 >                                         &muonReferencePreSelection,
468 >                                         &muonIDPFSelection,
469 >                                         &muonReferenceIsoSelection);
470 >
471 >      /*
472 >        // evolved
473 >        apply_HZZ4L_reference2_selection(ctrl, info,
474 >                                         vtx,
475 >                                         pfArr,
476 >                                         //rho,
477 >                                         puDArr,
478 >                                         electronArr,
479 >                                         &electronReferencePreSelection,
480 >                                         //&electronPreSelection,
481 >                                         &electronReferenceIDMVASelection,
482 >                                         &electronIsoMVASelection,
483 >                                         muonArr,
484 >                                         &muonReferencePreSelection,
485 >                                         //&muonPreSelection,
486 >                                         &muonIDPFSelection,
487 >                                         &muonIsoMVASelection);
488 >      */
489  
490        if( ctrl.debug ) cout << endl;
491  
492 <      //      cout << "bits: " << ret4l.status.selectionBits << endl;
492 >
493 >      int Z1type=0, Z2type=0;
494 >      if( ret4l.status.selectionBits.to_ulong() >= PASS_ZCANDIDATE ) {
495 >        if( abs(ret4l.Z1leptons[0].type) == 11  ) Z1type = 11;
496 >        if( abs(ret4l.Z1leptons[0].type) == 13  ) Z1type = 13;
497 >      }  
498 >      if( ret4l.status.selectionBits.to_ulong() >= PASS_ZZCANDIDATE ) {
499 >        if( abs(ret4l.Z2leptons[0].type) == 11  ) Z2type = 11;
500 >        if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
501 >      }  
502 >      
503 > #ifdef SYNC
504 >      cout  << "bits: " << ret4l.status.selectionBits  << "\t"
505 >            << "Z1: " << Z1type << "\t"
506 >            << "Z2: " << Z2type << endl;
507 >      cout << endl;
508 > #endif      
509        
510        if( ret4l.status.pass() ) {
370        
511          fillAngles(ret4l,angles);
512          fillKinematics(ret4l,kinematics);
513          fillEventInfo(info,evtinfo);
# Line 399 | Line 539 | int main(int argc, char** argv)
539          pass++;
540          //      if( pass > 3 ) break;
541  
542 <      } else if( ret4l.status.selectionBits.test(PASS_ZCANDIDATE) ) { // save for fakes ...
543 <        cout << "failingLeptons : " << failingLeptons.size() << endl;
544 <        for( int f=0; f<failingLeptons.size(); f++ ) {
545 <          SimpleLepton  l = failingLeptons[f];
546 <          cout << "f: " << f << "\t"
547 <               << "type: " << l.type << "\t"
548 <               << "pT: " << l.vec->Pt()
549 <               << endl;
550 <        }
551 <        /*
552 <        cout << "passingLeptons : " << passingLeptons.GetSize() << endl;
413 <        for( int f=0; f<passingLeptons.GetSize(); f++ ) {
414 <          cout << "f: " << f << "\t"
415 <               << "type: " << passingLeptons.At(f).type << "\t"
416 <               << "pT: " << passingLeptons.At(f).vec->Pt()
417 <               << endl;
542 >      } else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
543 >        //      if(ctrl.debug) {
544 >          cout << "failingLeptons : " << failingLeptons.size() << endl;
545 >          for( int f=0; f<failingLeptons.size(); f++ ) {
546 >            SimpleLepton  l = failingLeptons[f];
547 >            cout << "f: " << f << "\t"
548 >                 << "type: " << l.type << "\t"
549 >                 << "pT: " << l.vec.Pt()  << "\t"
550 >                 << "eta: " << l.vec.Eta()
551 >                 << endl;
552 >            //    }
553          }
419        */
420        cout << endl;
554          foTree.Fill();
555        }
556      }  
# Line 426 | Line 559 | int main(int argc, char** argv)
559    foTree.getFile()->cd();
560    foTree.getTree()->Write();
561    nt.WriteClose();
429 }
562  
563 < //----------------------------------------------------------------------------
564 < //
565 < // Get primary vertices
566 < // Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
567 < // NOTE: if no PV is found from fitting tracks, the beamspot is used
568 < //
569 < //----------------------------------------------------------------------------
570 < bool setPV(ControlFlags ctrl,
571 <           const mithep::Array<mithep::Vertex> * vtxArr,
572 <           mithep::Vertex & vtx)
573 < //----------------------------------------------------------------------------
574 < {
575 <  const Vertex *bestPV = 0;
576 <  float best_sumpt=-1;
577 <
578 <  // good PV requirements
579 <  const UInt_t   fMinNTracksFit = 0;
580 <  const Double_t fMinNdof       = 4;
449 <  const Double_t fMaxAbsZ       = 24;
450 <  const Double_t fMaxRho        = 2;
451 <  
452 <  for(int i=0; i<vtxArr->GetEntries(); ++i) {
453 <    const Vertex *pv = (Vertex*)(vtxArr->At(i));
454 <    if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
455 <    
456 <    // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
457 <    //  if(!pv->IsValid()) continue;
458 <    if(pv->NTracksFit()       < fMinNTracksFit)       continue;
459 <    if(pv->Ndof()                 < fMinNdof)         continue;
460 <    if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
461 <    if(pv->Position().Rho()   > fMaxRho)              continue;    
462 <    
463 <    // take the first ...
464 <    bestPV = pv;
465 <    break;
466 <
467 <    // this never reached ...    
468 <    float tmp_sumpt=0;
469 <    for( int t=0; t<pv->NTracks(); t++ )
470 <      tmp_sumpt += pv->Trk(t)->Pt();
471 <    
472 <    if( tmp_sumpt > best_sumpt  ) {
473 <      bestPV = pv;
474 <      best_sumpt = tmp_sumpt;
475 <      if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
476 <    }
563 >  cout << endl;
564 >  cout << endl;
565 >  for( int i=0; i<cutvec.size(); i++ ) {
566 >    cout << "cut: " << i << "\t"
567 >         << "pass: " << cutvec[i];
568 >
569 >    if( i>PASS_TRIGGER&&i<=PASS_GOODZ1 )
570 >      cout << "\t2e: " << zcutvec[0][i]
571 >           << "\t2m: " << zcutvec[1][i];
572 >
573 >    if( i>PASS_ZCANDIDATE )
574 >      cout << "\t4e: " << zzcutvec[0][i]
575 >           << "\t4m: " << zzcutvec[1][i]
576 >           << "\t2e2m: " << zzcutvec[2][i];
577 >
578 >    cout << "\t" << cutstrs[i] << endl;
579 >
580 >    cout << endl;
581    }
582 <  
583 <  if(!bestPV) bestPV = vtxArr->At(0);
480 <  vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
481 <  vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
482 <  return true;
483 < };
582 >
583 > }
584  
585  
586  
# Line 502 | Line 602 | void setEffiencyWeights(EventData & evtd
602    
603    for( int k=0; k<lepvec.size(); k++ ) {
604      if( abs(lepvec[k].type) == 13 ) {
605 <      mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec->Eta()), lepvec[k].vec->Pt()) );
606 <      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec->Eta()),
607 <                                                          lepvec[k].vec->Pt() ) );
605 >      mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
606 >      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
607 >                                                          lepvec[k].vec.Pt() ) );
608      } else {
609        
610        // now deal with medium vs loose
611        //              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
612        
613 <      std::pair<float,float> tmppair(fabs(lepvec[k].vec->Eta()), lepvec[k].vec->Pt());
613 >      std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
614        evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
615        
616 <      //              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec->Eta()),
617 <      //                                                                 lepvec[k].vec->Pt() ) );
616 >      //              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
617 >      //                                                                 lepvec[k].vec.Pt() ) );
618        
619 <      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec->Eta()),
620 <                                                          lepvec[k].vec->Pt(),
619 >      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
620 >                                                          lepvec[k].vec.Pt(),
621                                                            lepvec[k].isTight ) );
622        
623      }
# Line 563 | Line 663 | double getPUWeight(unsigned npu)
663   {
664    return hpu->GetBinContent(hpu->FindBin(npu));
665   }
666 +
667 + //----------------------------------------------------------------------------
668 + void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
669 + //----------------------------------------------------------------------------
670 + {
671 +  unsigned evt;
672 +  float rho;
673 +
674 +  cout << "initialzing EvtRhoMap ";
675 +  //FILE * f = fopen("evtrho.dat", "r");
676 +  //  FILE * f = fopen("allrho.cali.uniq.dat", "r");
677 +  FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
678 +  int lcount=0;
679 +  while( fscanf( f, "%u:%f", &evt, &rho ) ) {
680 +    if( feof(f) ) break;
681 +    evtrhoMap[evt] = rho;
682 +    lcount++;
683 +    if( !(lcount%1000) ) cout << "."; flush(cout);
684 +  }
685 +  cout << " done" << endl;
686 + }
687 +
688 + //----------------------------------------------------------------------------
689 + unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
690 +                     vector<bool> &pfNoPileUpflag,
691 +                     const mithep::VertexCol      * vtxArr )
692 + //----------------------------------------------------------------------------
693 + {
694 +  for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
695 +    const PFCandidate *pf = fPFCandidates->At(i);
696 +    assert(pf);
697 +
698 +    if(pf->PFType() == PFCandidate::eHadron) {
699 +      if(pf->HasTrackerTrk() &&
700 +         vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
701 +         vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
702 +
703 +        pfNoPileUpflag.push_back(1);
704 +
705 +      } else {
706 +
707 +        Bool_t vertexFound = kFALSE;
708 +        const Vertex *closestVtx = 0;
709 +        Double_t dzmin = 10000;
710 +        
711 +        // loop over vertices
712 +        for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
713 +          const Vertex *vtx = vtxArr->At(j);
714 +          assert(vtx);
715 +          
716 +          if(pf->HasTrackerTrk() &&
717 +             vtx->HasTrack(pf->TrackerTrk()) &&
718 +             vtx->TrackWeight(pf->TrackerTrk()) > 0) {
719 +            vertexFound = kTRUE;
720 +            closestVtx = vtx;
721 +            break;
722 +          }
723 +          Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
724 +          if(dz < dzmin) {
725 +            closestVtx = vtx;
726 +            dzmin = dz;
727 +          }
728 +        }
729 +
730 +        //      if(fCheckClosestZVertex) {
731 +        if(1) {
732 +          // Fallback: if track is not associated with any vertex,
733 +          // associate it with the vertex closest in z
734 +          if(vertexFound || closestVtx != vtxArr->At(0)) {
735 +            //      pfPileUp->Add(pf);
736 +            pfNoPileUpflag.push_back(0);
737 +          } else {
738 +            pfNoPileUpflag.push_back(1);
739 +          }
740 +        } else {
741 +          if(vertexFound && closestVtx != vtxArr->At(0)) {
742 +            //      pfPileUp->Add(pf);
743 +            pfNoPileUpflag.push_back(0);
744 +          } else {
745 +            //      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
746 +            pfNoPileUpflag.push_back(1);
747 +          }
748 +        }
749 +      } //hadron & trk stuff
750 +    } else { // hadron
751 +      //      PFCandidate * ptr = pfNoPileUp->AddNew();
752 +      pfNoPileUpflag.push_back(1);
753 +    }
754 +  } // Loop over PF candidates
755 +
756 +  return pfNoPileUpflag.size();
757 + }
758 +

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines