Selection/src/SelectionFuncs.cc

#include "SelectionFuncs.h"

extern vector<SimpleLepton> failingLeptons;
extern vector<SimpleLepton> passingLeptons;

extern vector<unsigned> cutvec;
extern vector<vector<unsigned> > zcutvec;
extern vector<vector<unsigned> > zzcutvec;
extern map<unsigned,float>       evtrhoMap;
extern bool passes_HLT_MC;
extern RunLumiRangeMap rlrm;

//----------------------------------------------------------------------------------------
void increment(map<TString, map<TString,int>* > &counterMap, TString cutName, int z1type, int z2type)
//----------------------------------------------------------------------------------------
// increment counters that keep track of event counts at each cut
{
  assert(counterMap.find(cutName) != counterMap.end());
  (*(counterMap[cutName]))["all"] ++;
  if(z1type!=0 && z2type!=0) {
    assert(counterMap.find(cutName) != counterMap.end());
    (*(counterMap[cutName]))[getChannel(z1type,z2type)] ++;
  }
}
//----------------------------------------------------------------------------
bool setPV(ControlFlags &ctrl, 
           const mithep::Array<mithep::Vertex> * vtxArr,
           const mithep::Vertex* &vtx) 
//----------------------------------------------------------------------------
{

  const mithep::Vertex *bestPV = 0;    
  float best_sumpt=-1;

  // good PV requirements
  const UInt_t   fMinNTracksFit = 0;
  const Double_t fMinNdof       = 4;
  const Double_t fMaxAbsZ       = 24;
  const Double_t fMaxRho        = 2;
  
  for(int i=0; i<vtxArr->GetEntries(); ++i) {
    const mithep::Vertex *pv = (mithep::Vertex*)(vtxArr->At(i));
    if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
    
    // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
    if(!pv->IsValid())                                continue;
    if(pv->NTracksFit()       < fMinNTracksFit)       continue;
    if(pv->Ndof()                 < fMinNdof)         continue;
    if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
    if(pv->Position().Rho()   > fMaxRho)              continue;    
    
    // take the first ...
    bestPV = pv;
    break;

    // this never reached ...     
    float tmp_sumpt=0;
    for( int t=0; t<pv->NTracks(); t++ )
      tmp_sumpt += pv->Trk(t)->Pt();
    
    if( tmp_sumpt > best_sumpt  ) { 
      bestPV = pv;
      best_sumpt = tmp_sumpt;
      if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
    }
  }
        
  // sync
  if(!bestPV) 
    return false;
  else {        
    vtx = bestPV;
    return true;
  }     
}
//----------------------------------------------------------------------------------------
void writeEntries(ControlFlags ctrl, unsigned total_unskimmed)
{
  ofstream ofs;
  TString ofname(ctrl.outputfile);
  ofs.open(ofname.ReplaceAll(".root",".nevents"));
  assert(ofs.is_open());
  ofs << total_unskimmed << endl;
  ofs.close();
}
//----------------------------------------------------------------------------------------
void getEATargets(ControlFlags &ctrl, mithep::MuonTools::EMuonEffectiveAreaTarget &eraMu, mithep::ElectronTools::EElectronEffectiveAreaTarget &eraEle)
{
  if( !ctrl.mc && ctrl.era == 2011 ) { 
    eraMu  = mithep::MuonTools::kMuEAData2011;
    eraEle = mithep::ElectronTools::kEleEAData2011;
  } else if( !ctrl.mc && ctrl.era == 2012 ) { 
    eraMu  = mithep::MuonTools::kMuEAData2012;
    eraEle = mithep::ElectronTools::kEleEAData2012;
  } else if( ctrl.mc && ctrl.era == 2011 ) { 
    eraMu  = mithep::MuonTools::kMuEAFall11MC;
    eraEle = mithep::ElectronTools::kEleEAFall11MC;
  } else if( ctrl.mc && ctrl.era == 2012 ) { 
    eraMu  = mithep::MuonTools::kMuEAData2012;
    eraEle = mithep::ElectronTools::kEleEAData2012;
  } else { 
    cerr << "unknown era for effective areas ... quitting." << endl;
    exit(1);
  }
}
//----------------------------------------------------------------------------
unsigned makePFnoPUArray(const mithep::Array<PFCandidate> * fPFCandidates,
                         vector<bool> &pfNoPileUpflag,
                         const mithep::Array<mithep::Vertex>      * vtxArr )
//----------------------------------------------------------------------------
//
// Make a vector of ints, one for each pf candidate, where we only include it in
// the isolation sum if it's 1
//
{
  assert(pfNoPileUpflag.size() == 0); // make sure we remembered to clear it after the last event
  for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
    const PFCandidate *pf = fPFCandidates->At(i);
    assert(pf);

    if(pf->PFType() == PFCandidate::eHadron) {
      if(pf->HasTrackerTrk() && 
         vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
         vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {

        pfNoPileUpflag.push_back(1);

      } else { 

        Bool_t vertexFound = kFALSE;
        const Vertex *closestVtx = 0;
        Double_t dzmin = 10000;
        
        // loop over vertices
        for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
          const Vertex *vtx = vtxArr->At(j);
          assert(vtx);
          
          if(pf->HasTrackerTrk() && 
             vtx->HasTrack(pf->TrackerTrk()) &&
             vtx->TrackWeight(pf->TrackerTrk()) > 0) {
            vertexFound = kTRUE;
            closestVtx = vtx;
            break;
          }
          Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
          if(dz < dzmin) {
            closestVtx = vtx;
            dzmin = dz;
          }
        }

        //      if(fCheckClosestZVertex) {
        if(1) {
          // Fallback: if track is not associated with any vertex,
          // associate it with the vertex closest in z
          if(vertexFound || closestVtx != vtxArr->At(0)) {
            //      pfPileUp->Add(pf);
            pfNoPileUpflag.push_back(0); // i.e., don't use it in the isolation sum
          } else {
            pfNoPileUpflag.push_back(1); // use it
          }
        } else {
          if(vertexFound && closestVtx != vtxArr->At(0)) {
            //      pfPileUp->Add(pf);
            pfNoPileUpflag.push_back(0);
          } else {
            //      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
            pfNoPileUpflag.push_back(1);
          }
        }
      } //hadron & trk stuff
    } else { // hadron
      //      PFCandidate * ptr = pfNoPileUp->AddNew();
      pfNoPileUpflag.push_back(1);
    }
  } // Loop over PF candidates

  return pfNoPileUpflag.size();
}
//--------------------------------------------------------------------------------
void setEfficiencyWeights(ControlFlags &ctrl,
                        unsigned era, 
                        EventData & evtdat, 
                        std::bitset<1024> triggerBits,
                        mithep::TriggerTable *hltTable,
                        mithep::Array<mithep::TriggerObject> *hltObjArr, 
                        mithep::TriggerObjectsTable *fTrigObjs,
                        WeightStruct & weights ) 
//--------------------------------------------------------------------------------
{
  //
  // offline weights
  //
  vector<SimpleLepton> lepvec = evtdat.Z1leptons;
  lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());

  // get eff for each leg, for each lepton
  vector< pair <double,double> > wlegs; // pair here is eff & err
  for( int k=0; k<lepvec.size(); k++ ) {
    double absEta = fabs(lepvec[k].vec.Eta());
    double pt = lepvec[k].vec.Pt();
    if( abs(lepvec[k].type) == 13 )
      wlegs.push_back( muonPerLegOfflineEfficiencyWeight( ctrl, era, absEta, pt ) );
    else
      wlegs.push_back( elePerLegOfflineEfficiencyWeight(  ctrl, era, absEta, pt, lepvec[k].isTight ) );
  }
  // combine weights and propagate uncertainties
  pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs ); 
  weights.woff    = offpair.first;
  weights.werroff = offpair.second;
  
  //
  // online weights
  //
  vector<SimpleLepton> muvec, elvec; 
  if(abs(evtdat.Z1leptons[0].type) == 11 ) { 
    elvec.push_back(evtdat.Z1leptons[0]);
    elvec.push_back(evtdat.Z1leptons[1]);
  }  else {
    muvec.push_back(evtdat.Z1leptons[0]);
    muvec.push_back(evtdat.Z1leptons[1]);
  }
  if(abs(evtdat.Z2leptons[0].type) == 11 ) { 
    elvec.push_back(evtdat.Z2leptons[0]);
    elvec.push_back(evtdat.Z2leptons[1]);
  } else {
    muvec.push_back(evtdat.Z2leptons[0]);
    muvec.push_back(evtdat.Z2leptons[1]);
  }
  pair<double,double> onpair  = getOnlineEfficiencyWeight( ctrl, triggerBits, hltTable, hltObjArr, fTrigObjs,  muvec, elvec );
  weights.won    = onpair.first;
  weights.werron = onpair.second;
} 
//----------------------------------------------------------------------------
void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
//----------------------------------------------------------------------------
{
  unsigned evt;
  float rho;

  cout << "initialzing EvtRhoMap ";
  //FILE * f = fopen("evtrho.dat", "r");
  //  FILE * f = fopen("allrho.cali.uniq.dat", "r");
  // FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
  assert(fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat","r")); FILE * f = fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
  int lcount=0;
  while( fscanf( f, "%u:%f", &evt, &rho ) ) {
    if( feof(f) ) break;
    evtrhoMap[evt] = rho;
    lcount++;
    if( !(lcount%1000) ) cout << "."; flush(cout);
  }
  cout << " done" << endl;
}

//----------------------------------------------------------------------------------------
void setEra(string fname, ControlFlags &ctrl) 
{
//----------------------------------------------------------------------------------------
  if( ctrl.mc && (fname.find("f11-") != string::npos)  ) ctrl.era=2011 ;
  else if( ctrl.mc && (fname.find("s12-") != string::npos)  ) ctrl.era=2012 ;
  else if( !ctrl.mc && (fname.find("r11") != string::npos)  ) ctrl.era=2011 ;
  else if( !ctrl.mc && (fname.find("r12") != string::npos)  ) ctrl.era=2012 ;
  else assert(0);
}

//----------------------------------------------------------------------------------------
unsigned getNPU(mithep::Array<mithep::PileupInfo> *puArr, int bx)
//----------------------------------------------------------------------------------------
{

  // Identify bunch crossing 0 in PU info
  int ibx0=-1, ibxm=-1, ibxp=-1;
  assert(puArr);    
  for(unsigned i=0; i<puArr->GetEntries(); ++i) {
    if(puArr->At(i)->GetBunchCrossing() == 0) ibx0=i;
    if(puArr->At(i)->GetBunchCrossing() ==-1) ibxm=i;
    if(puArr->At(i)->GetBunchCrossing() == 1) ibxp=i;
  }

  if(bx==0)       return puArr->At(ibx0)->GetPU_NumInteractions();
  else if(bx==-1) return puArr->At(ibxm)->GetPU_NumInteractions();
  else if(bx== 1) return puArr->At(ibxp)->GetPU_NumInteractions();
  else assert(0);

  return 0;
}


//----------------------------------------------------------------------------
void initPUWeights() 
//----------------------------------------------------------------------------
{
  // TFile * puf;

  // puf = new TFile("data/pileup/PUWeights_F11To2011.root");
  // hpu_2011 = (TH1D*)(puf->Get("puWeights"));
  // hpu_2011->SetDirectory(0);
  // puf->Close();

  // puf = new TFile("data/pileup/PUWeights_S12To2012_190456-199011.root");
  // hpu_2012 = (TH1D*)(puf->Get("puWeights"));
  // hpu_2012->SetDirectory(0);
  // puf->Close();
}

//----------------------------------------------------------------------------
double getPUWeight(int era, TString fname, unsigned npu)
//----------------------------------------------------------------------------
{
  if(npu > 59) npu = 59; // hardcoded crap in the pu functions
  
  if(era==2011) 
    return weightTrue2011(npu);
  else assert(era==2012);

  // assume that if it's in /029/ that we want to reweigh to 53x data
  TString vers,dataVers,mcVers;
  if(fname.Contains("/028/"))      dataVers = mcVers = vers = "52x";
  else if(fname.Contains("/029/")) dataVers = mcVers = vers = "53x";
  else {
    // if(ctrl.debug)
    //   cout << "WARNING in getPUWeight(): couldn't figure out if this is 52x or 53x, so just using 53x" << endl;
    dataVers = mcVers = vers = "53x";
  }

  if(dataVers=="52x" && mcVers=="52x")
    return weightTruePileupV07toIchep52X(npu);
  else if(dataVers=="53x" && mcVers=="53x")
    return weightTruePileupSummer12toMoriond(npu);
//   else if(dataVers=="53x" && mcVers=="53x")
//     return weightTruePileupV10toHcp53X(npu);
  // else if(dataVers=="53x" && mcVers=="52x")  // umm... can't really do this a.t.m..., but this would be the proper function
  //   return weightTruePileupV07toHcp53X(npu);
  else { cout << "Couldn't work out which npu to use!" << endl; assert(0); };

  return 1;
}
//----------------------------------------------------------------------------------------
void initJets(ControlFlags ctrl,
              string cmsswpath,
              vector<TString> &fJetCorrParsv,
              vector<JetCorrectorParameters> &correctionParameters,
              FactorizedJetCorrector *&fJetCorrector,
              JetIDMVA *&fJetIDMVA)
//----------------------------------------------------------------------------------------
{
  //
  // jet energy corrections
  //
  vector<TString> corrTypes;
  corrTypes.push_back("L1FastJet");
  corrTypes.push_back("L2Relative");
  corrTypes.push_back("L3Absolute");
  if(!ctrl.mc)
    corrTypes.push_back("L2L3Residual");

  TString tag;
  if(ctrl.era==2011) {
    cout << "\n\nNOTE: dont have the right JEC tags yet!\n\n" << endl;
    
    //these are the right tags for 42X monte carlo
    if(ctrl.mc) tag = "START42_V14B"; 
    else        tag = "GR_R_42_V23"; //same as GR_R_44_V15C and GR_R_42_V25;
  } else if(ctrl.era==2012) {
    cout << "\n\nNOTE: dont have the right JEC tags yet!\n\n" << endl;
    if(ctrl.mc) tag = "START53_V7F";   //same as START53_V10
    else        tag = "GR_P_V41_AN2";  //same as GR_P_V41_AN1
    assert(TString(cmsswpath).Contains("CMSSW_5_3"));
  } else assert(0);

  TString dataPath = TString(cmsswpath) + "/MitPhysics/data/";
  for(unsigned icorr=0; icorr<corrTypes.size(); icorr++) {
    if(ctrl.era==2011)       corrTypes[icorr] = "AK5PF_" + corrTypes[icorr];
    else if(ctrl.era==2012)  corrTypes[icorr] = corrTypes[icorr] + "_AK5PF";
    else assert(0);
    TString fname(dataPath + tag + "_" + corrTypes[icorr] + ".txt");
    if( !fopen(fname, "r") ) { cout << " no such file! " << fname << endl; assert(0); }
    addJetCorr(fname, fJetCorrParsv);
  }

  for(unsigned icorr=0; icorr<fJetCorrParsv.size(); icorr++)
    correctionParameters.push_back(JetCorrectorParameters(fJetCorrParsv[icorr].Data()));
  fJetCorrector = new FactorizedJetCorrector(correctionParameters);

  //
  // jet ID
  //
  fJetIDMVA = new JetIDMVA();
  JetIDMVA::MVAType mvaType;
  TString loPtXml(dataPath+"mva_JetID_lowpt.weights.xml");
  TString cfi(cmsswpath+"/MitPhysics/Utils/python/JetIdParams_hzz4l_moriond_cfi.py");
  //TString cfi(cmsswpath+"/MitPhysics/Utils/python/JetIdParams_cfi.py");
  TString Xml;
  if(ctrl.era==2011) {
    mvaType = JetIDMVA::k42;
    Xml = dataPath+"TMVAClassification_PuJetIdOptMVA.weights.xml";
  } else if(ctrl.era==2012) {
    mvaType = JetIDMVA::k52;
    Xml = dataPath+"TMVAClassification_5x_BDT_fullPlusRMS.weights.xml";
  } else assert(0);

  fJetIDMVA->Initialize(JetIDMVA::kLoose, loPtXml, Xml, mvaType, cfi);
}
//----------------------------------------------------------------------------------------
int getGenSampleType(TString ofname)
{
  if(ofname.Contains(TRegexp("-zz4[em]-")))               return ESampleType::kVV;  // powheg ZZ
  else if(ofname.Contains(TRegexp("-zz2e2m-")))           return ESampleType::kVV;  // ggZZ
  else if(ofname.Contains(TRegexp("-zz4l-gf")))           return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV;  // ggZZ
  else if(ofname.Contains(TRegexp("-zz2l2l-gf")))         return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV;  // ggZZ
  else if(ofname.Contains(TRegexp("-ggzz4l-")))           return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV;  // ggZZ
  else if(ofname.Contains(TRegexp("-ggzz2l2l-")))         return -1; // damn ggZZ has differences so MC parsing code fails ESampleType::kVV;  // ggZZ
  else if(ofname.Contains(TRegexp("-x12[0-9]zz4l-")))     return ESampleType::kHZZJHU; // spin samples
  else if(ofname.Contains(TRegexp("-h[0-9][0-9]*zz4l-"))) return ESampleType::kHZZ; // powheg HZZ
  else if(ofname.Contains(TRegexp("-h[0-9][0-9]*-")))     return ESampleType::kHZZ; // powheg HZZ where I fucked up the name
  else if(ofname.Contains(TRegexp("-zllm50-")))           return ESampleType::kZ;   // drell-yan
  else if(ofname.Contains(TRegexp("-zllm1050-")))         return ESampleType::kZ;   // drell-yan
  else if(ofname.Contains(TRegexp("-zjets-")))            return ESampleType::kZ;   // drell-yan
  else assert(0);
}
//----------------------------------------------------------------------------------------
TString getChannel(int z1type, int z2type)
{
  if(      abs(z1type) == 11 && abs(z2type) == 11 )  return "4e";//0;
  else if( abs(z1type) == 13 && abs(z2type) == 13 )  return "4m";//1;
  else                                               return "2e2m";//2;
}  
//----------------------------------------------------------------------------------------
// START: functions from Mangano
//----------------------------------------------------------------------------------------
float weightTruePileupSummer12toMoriond(float input){


  double w[60] = {
    0.252009, 
    0.327043, 
    0.335502, 
    0.352291, 
    0.324698, 
    0.582754, 
    0.455286, 
    0.441035, 
    0.607629, 
    0.930019, 
    1.3379, 
    1.69521, 
    1.74041, 
    1.54857, 
    1.32193, 
    1.15754, 
    1.07437, 
    1.05152, 
    1.07105, 
    1.11463, 
    1.15493, 
    1.17791, 
    1.18516, 
    1.17815, 
    1.1528, 
    1.10551, 
    1.03652, 
    0.948613, 
    0.844304, 
    0.728083, 
    0.607242, 
    0.489813, 
    0.381767, 
    0.287126, 
    0.207777, 
    0.144102, 
    0.0957851, 
    0.0611744, 
    0.0376984, 
    0.0226007, 
    0.0133203, 
    0.00782018, 
    0.00464555, 
    0.00284065, 
    0.00182028, 
    0.00123555, 
    0.000891118, 
    0.000679799, 
    0.000543107, 
    0.000449514, 
    0.000382089, 
    0.000331034, 
    0.000290923, 
    0.00025824, 
    0.000230472, 
    0.000206238, 
    0.000184523, 
    0.000164717, 
    0.000146364, 
    0.000271879};

  //return w[(int)floor(input+0.5)];
  return w[(int)floor(input)];    
}


float weightTruePileupV07toIchep52X(float input){
  float w[60] = {
    2.59113e-07,
    1.77107e-05,
    0.0263017,
    1.66396,
    0.16421,
    0.754166,
    2.84622,
    5.57103,
    8.66558,
    11.5716,
    11.8712,
    12.8102,
    10.3421,
    8.91019,
    7.614,
    6.10397,
    4.52745,
    3.26031,
    2.39558,
    1.83297,
    1.47821,
    1.26728,
    1.14716,
    1.06707,
    0.98066,
    0.860877,
    0.708281,
    0.539789,
    0.37652,
    0.237298,
    0.1338,
    0.0671236,
    0.0299236,
    0.0118723,
    0.00420968,
    0.00134235,
    0.000389563,
    0.000104892,
    2.69214e-05,
    6.79674e-06,
    1.73307e-06,
    4.52553e-07,
    1.21124e-07,
    3.29924e-08,
    9.10616e-09,
    2.53998e-09,
    7.16146e-10,
    2.03786e-10,
    5.84308e-11,
    1.68192e-11,
    4.8434e-12,
    1.38959e-12,
    3.96112e-13,
    1.11358e-13,
    3.17245e-14,
    5.34916e-15,
    0,
    0,
    0,
    0};
  
  //return w[(int)floor(input+0.5)];
  return w[(int)floor(input)];
    
    /*   
         TH1F h("boh","boh",60,0.,60.); 
         for(int k=0;k<60;k++){
         h.SetBinContent(k+1,w[k]);
         } 
         return h.GetBinContent(h.FindBin(input));
    */
}


float weightTruePileupV07toHcp53X(float input){
  float w[60] = {
    0.0447136,     
    0.11785,       
    0.23825,
    1.08447,
    0.102575,
    0.454605,
    1.79761,
    4.00271,
    6.83281,
    9.83701,
    10.7966,
    12.2356,
    10.0247,
    8.49395,
    7.1125,
    5.69527,
    4.31256,
    3.19305,
    2.42035,
    1.91666,
    1.58485,
    1.36297,
    1.21166,
    1.09466,
    0.978941,
    0.84653,
    0.699235,
    0.548996,
    0.408673,
    0.288194,
    0.193367,
    0.124653,
    0.0781124,
    0.0479268,
    0.0287763,
    0.0167744,
    0.00941834,
    0.00507877,
    0.00264364,
    0.00134612,
    0.000682678,
    0.000351412,
    0.0001864,
    0.00010259,
    5.87818e-05,
    3.5033e-05,
    2.17116e-05,
    1.39777e-05,
    9.36123e-06,
    6.53328e-06,
    4.76598e-06,
    3.64139e-06,
    2.92018e-06,
    2.4602e-06,
    2.17291e-06,
    2.01107e-06,
    1.94392e-06,
    1.9598e-06,
    2.0583e-06,
    2.24895e-06};

  return w[(int)floor(input)];
}


float weightTruePileupV10toIchep53X(float input){
  float w[60] = {
    2.35693e-06,
    7.51928e-05,
    0.0263529,
    0.609947,
    0.737917,
    1.29365,
    0.994503,
    0.85454,
    1.01559,
    1.33243,
    1.72454,
    2.01264,
    2.00573,
    1.80333,
    1.56328,
    1.37452,
    1.24753,
    1.16481,
    1.11738,
    1.09701,
    1.08843,
    1.08796,
    1.09768,
    1.10763,
    1.09328,
    1.0339,
    0.92408,
    0.771537,
    0.59283,
    0.41266,
    0.256892,
    0.14188,
    0.0692543,
    0.029902,
    0.0114564,
    0.00391383,
    0.00120625,
    0.000341485,
    9.09127e-05,
    2.34008e-05,
    5.95438e-06,
    1.5122e-06,
    3.82094e-07,
    9.51794e-08,
    2.32205e-08,
    5.51698e-09,
    1.27267e-09,
    2.84346e-10,
    6.12799e-11,
    1.26731e-11,
    2.50309e-12,
    4.69797e-13,
    8.35153e-14,
    1.39452e-14,
    2.24718e-15,
    2.03841e-16,
    0,
    0,
    0,
    0};

  return w[(int)floor(input)];
}


float weightTruePileupV10toHcp53X(float input){
  float w[60] = {
    0.409409,
    0.527276,
    0.39328,
    0.507892,
    0.48029,
    0.787701,
    0.632356,
    0.618033,
    0.806089,
    1.14018,
    1.5788,
    1.93507,
    1.957,
    1.73004,
    1.46737,
    1.28278,
    1.18189,
    1.13388,
    1.12578,
    1.14415,
    1.16048,
    1.1618,
    1.15318,
    1.13405,
    1.09239,
    1.01915,
    0.914837,
    0.786744,
    0.644879,
    0.502039,
    0.371688,
    0.263586,
    0.18067,
    0.120472,
    0.0780184,
    0.0486113,
    0.0289039,
    0.0163367,
    0.00879674,
    0.00456046,
    0.0023098,
    0.00115977,
    0.000583207,
    0.000294815,
    0.000149865,
    7.62892e-05,
    3.87537e-05,
    1.96105e-05,
    9.87744e-06,
    4.95418e-06,
    2.47913e-06,
    1.23919e-06,
    6.19751e-07,
    3.10125e-07,
    1.54934e-07,
    7.71425e-08,
    3.8182e-08,
    1.87455e-08,
    9.10765e-09,
    9.19802e-09};
  return w[(int)floor(input)];
}


float weightTrue2011(float input){
  if(input>50) 
    return 1;

    
  float w[50];


w[0]= 0.212929;
w[1]= 0.0208114;
w[2]= 0.0584048;
w[3]= 0.538898;
w[4]= 1.357;
w[5]= 1.49913;
w[6]= 1.42247;
w[7]= 1.35904;
w[8]= 1.29946;
w[9]= 1.27925;
w[10]= 1.37845;
w[11]= 1.71246;
w[12]= 1.5291;
w[13]= 1.35234;
w[14]= 1.22215;
w[15]= 1.0155;
w[16]= 1.01137;
w[17]= 0.395465;
w[18]= 0.230984;
w[19]= 0.109883;
w[20]= 0.0433739;
w[21]= 0.0111497;
w[22]= 0.00408801;
w[23]= 0.00115678;
w[24]= 0.000365505;
w[25]= 0.000112391;
w[26]= 3.83894e-05;
w[27]= 1.60651e-05;
w[28]= 4.81412e-06;
w[29]= 1.39717e-06;
w[30]= 1.92368e-06;
w[31]= 4.10748e-06;
w[32]= 2.33157e-05;
w[33]= 4.0181e-05;
w[34]= 4.87786e-05;
w[35]= 0.00194128;
w[36]= 8.97414e-05;
w[37]= 1;
w[38]= 1;
w[39]= 0.000162709;
w[40]= 1;
w[41]= 1;
w[42]= 1;
w[43]= 1;
w[44]= 1;
w[45]= 1;
w[46]= 1;
w[47]= 1;
w[48]= 1;
w[49]= 1;


 TH1F h("boh","boh",50,0.,50.);
 
 for(int k=0;k<50;k++){
   h.SetBinContent(k+1,w[k]);
 }
 
 return h.GetBinContent(h.FindBin(input));

}
//----------------------------------------------------------------------------------------
float weightTrue2011to2012(float input){
  if(input>50) 
    return 1;
    
  float w[50];

w[0]= 0.000443112;
w[1]= 0.000248044;
w[2]= 0.000273111;
w[3]= 0.00109511;
w[4]= 0.00195699;
w[5]= 0.00480746;
w[6]= 0.027013;
w[7]= 0.074795;
w[8]= 0.166231;
w[9]= 0.309545;
w[10]= 0.577657;
w[11]= 1.12488;
w[12]= 1.36899;
w[13]= 1.56925;
w[14]= 1.89846;
w[15]= 2.20828;
w[16]= 3.14112;
w[17]= 1.87712;
w[18]= 1.97062;
w[19]= 2.07067;
w[20]= 2.17791;
w[21]= 1.7176;
w[22]= 2.10953;
w[23]= 2.0805;
w[24]= 2.29498;
w[25]= 2.42189;
w[26]= 2.80303;
w[27]= 3.94091;
w[28]= 3.67917;
w[29]= 2.26081;
w[30]= 2.99726;
w[31]= 3.76553;
w[32]= 11.285;
w[33]= 10.2781;
w[34]= 6.73407;
w[35]= 148.182;
w[36]= 3.88144;
w[37]= 1;
w[38]= 1;
w[39]= 1.48128;
w[40]= 1;
w[41]= 1;
w[42]= 1;
w[43]= 1;
w[44]= 1;
w[45]= 1;
w[46]= 1;
w[47]= 1;
w[48]= 1;
w[49]= 1;


 TH1F h("boh","boh",50,0.,50.);
 
 for(int k=0;k<50;k++){
   h.SetBinContent(k+1,w[k]);
 }
 
 return h.GetBinContent(h.FindBin(input));

}
//----------------------------------------------------------------------------------------
// END: functions from Mangano
//----------------------------------------------------------------------------------------
Revision:	1.20
Committed:	Wed Feb 6 22:31:34 2013 UTC (12 years, 3 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.19:	+1 -1 lines
Log Message:	fixed problems so MitHzz4l compiles