Angles/src/fillBdtTree.cc

#include <iostream>
#include <assert.h>
#include <stdio.h>     
#include <stdlib.h>    
#include <getopt.h>
#include <string.h>
#include <utility>
#include <map>
// #include <boost/regex.hpp>

#include "TFile.h"
#include "TChain.h"
#include "TList.h"
#include "TRegexp.h"
#include "TTree.h"
#include "TBranch.h"

#include "ZZMatrixElement/MELA/interface/Mela.h"
#include "ZZMatrixElement/MELA/interface/PseudoMELA.h"

#include "PwhgWrapper.h"

#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/Config.h"
#include "TMVA/MethodBDT.h"

#include "KinematicsStruct.h"
#include "AngleTuple.h"
#include "TH1F.h"
#include "filestuff.h"
#include "Angles.h"
#include "WeightStruct.h"
#include "GenInfoStruct.h"
#include "SampleWeight.h"
#include "varkeepter.h"

#ifndef CMSSW_BASE
#define CMSSW_BASE "../../"
#endif

using namespace std;

pwhegwrapper *wrap;

class CtrlFlags 
{ 
public :
  CtrlFlags():inputfile(""),
              inputtree("zznt"),
              weightfile(""),
              output("test.root"),
              mH_lo(-1),
              mH_hi(-1),
              wL(0),
              wH(0),
              varstring(""),
              makeMela(false),
              multiclass(false),
              flattree(false),
              multisigs(false),
              era(-1),
              withPtY(false),
              mH(0),
              hwidth(0)
  { };
  ~CtrlFlags() {};
  void dump();

  string inputfile;
  string inputtree;
  string weightfile;
  string output;
  float mH_lo,mH_hi;
  float wL,wH;
  TString varstring;
  bool makeMela;
  bool multiclass;
  bool flattree;
  bool multisigs;
  int era;
  bool withPtY;
  float mH;
  float hwidth;
};
void CtrlFlags::dump()
{
  cout << "mH           : " << mH << endl; 
  cout << "inputfile    : " << inputfile << endl; 
  cout << "inputtree    : " << inputtree << endl; 
  cout << "weightfile   : " << weightfile << endl; 
  cout << "output       : " << output << endl; 
  cout << "mH_lo        : " << mH_lo << endl; 
  cout << "mH_hi        : " << mH_hi << endl; 
  cout << "wL           : " << wL << endl; 
  cout << "wH           : " << wH << endl; 
  cout << "varstring    : " << varstring << endl; 
  cout << "makeMela     : " << makeMela << endl; 
  cout << "multiclass   : " << multiclass << endl; 
  cout << "flattree     : " << flattree << endl; 
  cout << "multisigs    : " << multisigs << endl; 
  cout << "withPtY      : " << withPtY << endl; 
  cout << "hwidth       : " << hwidth << endl; 
  assert(era==2011 || era==2012);
  cout << "era          : " << era << endl; 
}
PseudoMELA *psMela;
void parse_args( int, char**, CtrlFlags & );
void setFracs(float mH, float mH_lo, float mH_hi, float &wgt_lo, float &wgt_hi);
void setZzntVals(CtrlFlags &flags, varkeepter *flatoutput, filestuff *fs, TString cls, TMVA::Reader *reader=0, Mela *mela=0, TMVA::Reader *spinReader=0);
double getWeight(TH1F *hist, float val);
void setPtWgtHists(TString cls, TString fname, Long_t mH, TFile *wgtFile);
//----------------------------------------------------------------------------------------
map<TString,TH1F*> wgtHists;

//------------------------------------------------------------
int main(int argc, char ** argv)  
//------------------------------------------------------------
{
  CtrlFlags flags;
  parse_args( argc, argv, flags);
  flags.dump();
  assert(flags.varstring!="");
  assert(fopen(flags.weightfile.c_str(), "r") || flags.makeMela);

  TFile *pt4lWgts = TFile::Open("data/pt4l/pt4l_weights.root");
  wrap = new pwhegwrapper;
  
  TString cmsswpath(CMSSW_BASE + TString("/src"));
  string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
  SimpleTable xstab(xspath.c_str());

  //
  // setup input
  //

  cout << "flags: ----------------------------------------------------------------------------------------" << endl;
  cout << "inputfile:      " << flags.inputfile << endl;
  cout << "output:      " << flags.output << endl;

  TString spinWeightFile("");//ntupledir;
  vector<TString> classes,fileNames;
  vector<double> fileFracs,fileMasses;
  vector<TFile*> inFiles;
  vector<TTree*> inTrees;
  ifstream ifs(flags.inputfile.c_str());
  assert(ifs.is_open());
  string line;
  while(getline(ifs,line)) {
    if(line[0]=='#') continue;
    stringstream ss(line);
    if(line[0]=='^') {
      TString dummy;
      // if(TString(line).Contains("^ntupledir")) ss >> dummy >> ntupledir;
      if(TString(line).Contains("^spinWeightFile")) { ss >> dummy >> spinWeightFile; cout << "spinWeightFile: " << spinWeightFile << endl; }
      continue; 
    }

    TString cls,fname,type;
    double frac,fileMass;
    ss >> cls >> type >> frac >> fname >> fileMass;
    assert(cls=="gfH" || cls=="vbfH" || cls=="vttH" || cls=="qqH" || cls=="qqZZ" || cls=="ggZZ" || cls=="lljj" || cls=="data");
    classes.push_back(cls);
    fileNames.push_back(/*ntupledir+"/"+*/fname);
    fileFracs.push_back(frac);
    fileMasses.push_back(fileMass);
    inFiles.push_back(TFile::Open(fname));
    assert(inFiles.back()->IsOpen());
    inTrees.push_back((TTree*)inFiles.back()->Get(flags.inputtree.c_str()));
    assert(inTrees.back());
    cout << "    pushing back: " << cls << " " << setw(12) << frac << " " << fname << " " << fileMass << endl;
  }
  ifs.close();

  // mass cut based on weightfile name ... 
  float mH = flags.mH;
  // boost::regex expression2(".*_mH(.*)_.*");
  // boost::cmatch match2;
  // if(boost::regex_match(flags.weightfile.c_str(), match2, expression2))
  //   mH = std::strtof(match2[1].first,NULL);

  // assert(flags.wL>=0 && flags.wL < 1000 && flags.wH >= 0 && flags.wH < 1000);
  char cutstr[512];
  sprintf(cutstr, "m4l>%.3f&&m4l<%.3f", flags.wL, flags.wH);
  cout << "mH: "   << setw(8) << mH 
       << "  wL: " << setw(8) << flags.wL 
       << "  wH: " << setw(8) << flags.wH 
       << "  mH_lo: " << setw(8) << flags.mH_lo 
       << "  mH_hi: " << setw(8) << flags.mH_hi 
       << endl;

  float frac_lo=1,frac_hi=1;
  map<double,double> mH_fracs;
  mH_fracs[-1] = 1; // this means zz
  if(flags.mH_lo != 0) {
    setFracs(mH, flags.mH_lo,  flags.mH_hi, frac_lo, frac_hi);
    mH_fracs[flags.mH_lo] = frac_lo;
    mH_fracs[flags.mH_hi] = frac_hi;
  }
  mH_fracs[mH]    = 1;

  unsigned nMin=9999999;
  vector<double> passFracs;
  for(unsigned ifile=0; ifile<fileNames.size(); ifile++) {
    cout << classes[ifile] << ":  "
         << setw(8) << inTrees[ifile]->GetEntries(cutstr) << " / " << setw(8) << inTrees[ifile]->GetEntries() << " events pass cut in file: " << fileNames[ifile] << endl;
    int nEntriesWithCut = inTrees[ifile]->GetEntries(cutstr);
    if(nEntriesWithCut < nMin) nMin = nEntriesWithCut;
    passFracs.push_back( double(nEntriesWithCut)/inTrees[ifile]->GetEntries() );
  }

  cout << "twice the min number of entries over all the files, 2*nMin: " << 2*nMin << endl;

  TString tmpfname(flags.output);
  tmpfname.ReplaceAll("/","-");
  tmpfname = "/tmp/"+tmpfname+"-tmpfile.root";
  TFile tmpfile(tmpfname,"recreate");
  vector<TTree*> TreeCopies; // copies of the input trees, but with only the number of events we want
  TList *treeList = new TList;
  vector<long> nActuallyAdded,nActuallyAddedCumulative;
  cout << "making tree copies: " << endl;
  for(unsigned ifile=0; ifile<fileNames.size(); ifile++) {
    if(mH_fracs.find(fileMasses[ifile]) == mH_fracs.end())
      cout << "ERROR: filemasses[ifile]: " << fileMasses[ifile] << endl;
    assert(mH_fracs.find(fileMasses[ifile]) != mH_fracs.end());
    double fracThisFile,nTarget;
    long nToRequest;
    if(fileFracs[ifile] == 1) {            // use all of this sample (eg for background)
      fracThisFile = 1;
      nTarget      = passFracs[ifile]*inTrees[ifile]->GetEntries();
      nToRequest   = 1000000000; // use default value for TTree::CopyTree(): all of 'em
    } else {
      if(fileFracs[ifile] == -1)
        fracThisFile = mH_fracs[fileMasses[ifile]]; // get fraction for this sample from function
      else
        fracThisFile = fileFracs[ifile];            // get fraction for this sample from txt file
      nTarget      = 2*fracThisFile*nMin;
      nToRequest   = nTarget / passFracs[ifile];
    }
    cout << "  adding " << setw(12) << nTarget << " from file " << fileNames[ifile]
         << ", which means asking for: " << nToRequest << endl;
    TreeCopies.push_back(inTrees[ifile]->CopyTree(cutstr,"",nToRequest));

    long nCumulative=0;
    for(unsigned j=0; j<nActuallyAdded.size(); j++) nCumulative += nActuallyAdded[j];
    nActuallyAddedCumulative.push_back(nCumulative);
    nActuallyAdded.push_back(long(TreeCopies.back()->GetEntries()));
    cout << "actually added : " << nActuallyAdded.back() << ", cumulative: " << nActuallyAddedCumulative.back() << endl;

    treeList->Add(TreeCopies.back());
  }

  TTree * t_tmp = TTree::MergeTrees(treeList);
  tmpfile.Write();
  tmpfile.Close();

  //
  // setup output
  //
  filestuff fs(tmpfname,tmpfname,"",classes[0]=="data");
  float f_bdt; 
  TBranch *br_bdt=0;
  varkeepter *flatoutput=0;
  TString varfname("zzntvars");
  if(flags.multisigs) varfname += "-multisigs";
  flatoutput = new varkeepter("Angles/data/"+varfname+".txt","w",flags.output.c_str(),42,"zznt");

  //
  // setup TMVA
  //
  float reduced_mZ1, reduced_mZ2;
  float reduced_zzdotz1, reduced_zzdotz2;
  float reduced_ZZpt, reduced_Z1pt, reduced_Z2pt;

  map<string,bool> varmap;
  char * dovar = strtok(const_cast<char *>(flags.varstring.Data()), ":");
  if( dovar != NULL ) varmap[string(dovar)]=true;
  while( (dovar = strtok(NULL, ":")) != NULL ) {
    varmap[string(dovar)]=true;
  }

  int run,lumi,evt;
  TMVA::Reader *reader=0,*spinReader=0;
  Mela *mela=0;
  psMela=0;
  if(!flags.makeMela) {
    reader = new TMVA::Reader( "V" );
    if(varmap[string("costheta1")])     reader->AddVariable( "costheta1",               &fs.angles->costheta1);
    if(varmap[string("costheta2")])     reader->AddVariable( "costheta2",               &fs.angles->costheta2);
    if(varmap[string("costhetastar")])  reader->AddVariable( "costhetastar",            &fs.angles->costhetastar);
    if(varmap[string("Phi")])           reader->AddVariable( "Phi",                     &fs.angles->Phi);
    if(varmap[string("Phi1")])          reader->AddVariable( "Phi1",                    &fs.angles->Phi1);
    if(varmap[string("mZ1")])           reader->AddVariable( "mZ1",                     &reduced_mZ1);
    if(varmap[string("mZ2")])           reader->AddVariable( "mZ2",                     &reduced_mZ2);
    if(varmap[string("pt4l")])          reader->AddVariable( "ZZpt/m4l",                &reduced_ZZpt);
    if(varmap[string("zzdotz1")])       reader->AddVariable( "ZZdotZ1/(m4l*mZ1)",       &reduced_zzdotz1 );
    if(varmap[string("zzdotz2")])       reader->AddVariable( "ZZdotZ2/(m4l*mZ2)",       &reduced_zzdotz2 );
    if(varmap[string("dphi1")])         reader->AddVariable( "ZZptCosDphiZ1pt",         &fs.kine->ZZptZ1ptCosDphi );
    if(varmap[string("dphi2")])         reader->AddVariable( "ZZptCosDphiZ2pt",         &fs.kine->ZZptZ2ptCosDphi );
    if(varmap[string("Z1pt")])          reader->AddVariable( "Z1pt/m4l",                &reduced_Z1pt );
    if(varmap[string("Z2pt")])          reader->AddVariable( "Z2pt/m4l",                &reduced_Z2pt );
    if(varmap[string("y4l")])           reader->AddVariable( "ZZy",                     &fs.kine->ZZy);
    if(varmap[string("m4l")])             reader->AddVariable( "m4l",                     &fs.kine->m4l);
    else                                  reader->AddSpectator("m4l",                     &fs.kine->m4l);
    reader->BookMVA("BDTG", flags.weightfile.c_str());
    if(spinWeightFile != "") {
      spinReader = new TMVA::Reader( "V" );
      if(varmap[string("costheta1")])           spinReader->AddVariable( "costheta1",           &fs.angles->costheta1);
      if(varmap[string("costheta2")])           spinReader->AddVariable( "costheta2",           &fs.angles->costheta2);
      if(varmap[string("costhetastar")])        spinReader->AddVariable( "costhetastar",        &fs.angles->costhetastar);
      if(varmap[string("Phi")])                 spinReader->AddVariable( "Phi",                 &fs.angles->Phi);
      if(varmap[string("Phi1")])                spinReader->AddVariable( "Phi1",                &fs.angles->Phi1);
      if(varmap[string("mZ1")])                 spinReader->AddVariable( "mZ1",                 &reduced_mZ1);
      if(varmap[string("mZ2")])                 spinReader->AddVariable( "mZ2",                 &reduced_mZ2);
      // if(varmap[string("pt4l")])             spinReader->AddVariable( "ZZpt/m4l",            &reduced_ZZpt);
      // if(varmap[string("zzdotz1")])          spinReader->AddVariable( "ZZdotZ1/(m4l*mZ1)",   &reduced_zzdotz1 );
      // if(varmap[string("zzdotz2")])          spinReader->AddVariable( "ZZdotZ2/(m4l*mZ2)",   &reduced_zzdotz2 );
      // if(varmap[string("dphi1")])            spinReader->AddVariable( "ZZptCosDphiZ1pt",     &fs.kine->ZZptZ1ptCosDphi );
      // if(varmap[string("dphi2")])            spinReader->AddVariable( "ZZptCosDphiZ2pt",     &fs.kine->ZZptZ2ptCosDphi );
      // if(varmap[string("Z1pt")])             spinReader->AddVariable( "Z1pt/m4l",            &reduced_Z1pt );
      // if(varmap[string("Z2pt")])             spinReader->AddVariable( "Z2pt/m4l",            &reduced_Z2pt );
      // if(varmap[string("y4l")])                      spinReader->AddVariable( "ZZy",                 &fs.kine->ZZy);
      if(varmap[string("m4l")])                 spinReader->AddVariable( "m4l",                 &fs.kine->m4l);
      else                                      spinReader->AddSpectator("m4l",                 &fs.kine->m4l);
      spinReader->AddSpectator("run",                 &run);
      spinReader->AddSpectator("lumi",                &lumi);
      spinReader->AddSpectator("evt",                 &evt);
      spinReader->BookMVA("BDTG", spinWeightFile);
    }
  } else {
    // mela = new Mela::Mela(bool usePowhegTemplate=false, int LHCsqrts=8);
    mela = new Mela(false, flags.era==2011 ? 7 : 8);
    if(flags.multisigs)
      psMela = new PseudoMELA();
  }

  //
  // evaluate
  //
  vector<unsigned> runv,evtv; // for removing duplicate events
  int ifile=-1;
  int nent=fs.get_zz_chain()->GetEntries();
  cout << "beginning tree from mass " << mH << "(" << flags.wL << "," << flags.wH << ") with " << nent << "  entries" << endl;
  for( int i=0; i<nent; i++ ) { 
    if( !(i%5000) ) cerr << "processed:  " << i << "/" << nent << endl;
    fs.getentry(i,"","zznt");

    if(classes[0]=="data") { // the rest of the classes better be data as well...
      bool isDuplic=false;
      for(unsigned ievt=0; ievt<evtv.size(); ievt++) {
        if((fs.info->run == runv[ievt]) && (fs.info->evt == evtv[ievt])) {
          isDuplic = true;
          break;
        }
      }
      if(isDuplic) {
        cout << "   skipping duplicate event: " << fs.info->run << " " << fs.info->evt << endl;
        continue;
      }
    }
    runv.push_back(fs.info->run);
    evtv.push_back(fs.info->evt);

    if(ifile+1 < (int)nActuallyAddedCumulative.size() && i >= nActuallyAddedCumulative[ifile+1]) {
      ifile++;
      cout << " ientry: " << i << ", switching file to (" << classes[ifile] << "): " << fileMasses[ifile]
           << " so weight multiplies by: " << ((classes[ifile]=="gfH" || classes[ifile]=="vbfH" || classes[ifile]=="qqH") ? mH_fracs[fileMasses[ifile]] : 1) << endl;
      if(classes[ifile]=="gfH" || classes[ifile]=="qqZZ")
        setPtWgtHists(classes[ifile],fileNames[ifile],mH,pt4lWgts);
    }

    if( fs.kine->m4l>flags.wH || fs.kine->m4l<flags.wL ) { cout << "ERROR! shouldn't be here (already applied the cut): " << fs.kine->m4l << endl; continue; }

    reduced_mZ1         = fs.kine->mZ1;
    reduced_mZ2         = fs.kine->mZ2;
    reduced_ZZpt        = fs.kine->ZZpt/fs.kine->m4l;
    reduced_Z1pt        = fs.kine->Z1pt/fs.kine->m4l;
    reduced_Z2pt        = fs.kine->Z2pt/fs.kine->m4l;
    reduced_zzdotz1     = fs.kine->ZZdotZ1/(fs.kine->m4l*fs.kine->mZ1);
    reduced_zzdotz2     = fs.kine->ZZdotZ2/(fs.kine->m4l*fs.kine->mZ2);

    if(classes[ifile]=="gfH" || classes[ifile]=="vbfH" || classes[ifile]=="qqH")
      fs.weights->w *= mH_fracs[fileMasses[ifile]];

    if(!flags.makeMela && !flags.flattree) { // shouldn't ever use this at the point
      assert(0);
      fs.outtuple->Fill();
    } else {
      setZzntVals(flags, flatoutput, &fs, classes[ifile], reader, mela, spinReader);
      flatoutput->fill();
    }
  }

  if(!flags.makeMela && !flags.flattree) {
    fs.outtuple->getTree()->Print();
    fs.outtuple->WriteClose();
  } else {
    flatoutput->gettree()->Print();
    flatoutput->writeclose();
  }
}
//----------------------------------------------------------
void parse_args( int argc, char** argv, CtrlFlags &flags ) 
//----------------------------------------------------------
{
  int c;
  int digit_optind = 0;

  while (1) {
    int this_option_optind = optind ? optind : 1;
    int option_index = 0;
    static struct option long_options[] = {
      {"inputfile",  1, 0, 'a'},
      {"inputtree",  1, 0, 'b'},
      {"weightfile", 1, 0, 'c'},
      {"output", 1, 0, 'd'},
      {"mH_lo", 1, 0, 'e'},
      {"mH_hi", 1, 0, 'f'},
      {"wL", 1, 0, 'g'},
      {"wH", 1, 0, 'h'},
      {"varstring", 1, 0, 'i'},
      {"makeMela", 0, 0, 'j'},
      {"multiclass", 0, 0, 'k'},
      {"flattree", 0, 0, 'l'},
      {"multisigs", 0, 0, 'm'},
      {"era", 1, 0, 'n'},
      {"withPtY", 0, 0, 'o'},
      {"mH", 1, 0, 'p'},
      {"hwidth", 1, 0, 'q'},
      {0, 0, 0, 0}
    };
    
    c = getopt_long (argc, argv, "a:e",
                     long_options, &option_index);
    if (c == -1)
      break;
    
    switch (c) {
    case 0:
      printf ("option %s", long_options[option_index].name);
      if (optarg)
        printf (" with arg %s", optarg);
      printf ("\n");
      break;
      
    case 'a':
      flags.inputfile = std::string(optarg);
      break;
    case 'b':
      flags.inputtree = std::string(optarg);
      break;
    case 'c':
      flags.weightfile = std::string(optarg);
      break;
    case 'd':
      flags.output = std::string(optarg);
      break;
    case 'e':
      flags.mH_lo = strtof(optarg,NULL);
      break;
    case 'f':
      flags.mH_hi = strtof(optarg,NULL);
      break;
    case 'g':
      flags.wL = strtof(optarg,NULL);
      break;
    case 'h':
      flags.wH = strtof(optarg,NULL);
      break;
    case 'i':
      flags.varstring = TString(optarg);
      break;
    case 'j':
      flags.makeMela = true;
      break;
    case 'k':
      flags.multiclass = true;
      break;
    case 'l':
      flags.flattree = true;
      break;
    case 'm':
      flags.multisigs = true;
      break;
    case 'n':
      flags.era = atoi(optarg);
      assert(flags.era==2011 || flags.era==2012);
      break;
    case 'o':
      flags.withPtY = true;
      break;
    case 'p':
      flags.mH = strtof(optarg,NULL);
      break;
    case 'q':
      flags.hwidth = strtof(optarg,NULL);
      break;
    case '2':
      if (digit_optind != 0 && digit_optind != this_option_optind)
        printf ("digits occur in two different argv-elements.\n");
      digit_optind = this_option_optind;
      printf ("option %c\n", c);
      break;
      
    case '?':
      break;
      
    default:
      printf ("?? getopt returned character code 0%o ??\n", c);
    }
  }
  
  if( flags.inputfile.empty() ) {
    cerr << "please specify an inputfile with --inputfile " << endl;
    exit(1);
  }

  if (optind < argc) {
    printf ("non-option ARGV-elements: ");
    while (optind < argc)
      printf ("%s ", argv[optind++]);
    printf ("\n");
  }
}
//----------------------------------------------------------------------------------------
void setFracs(float mH, float mH_lo, float mH_hi, float &frac_lo, float &frac_hi)
{
  assert(mH_lo<=mH && mH_hi>=mH && mH_lo>0 && mH_hi>0);
  double delta = mH_hi - mH_lo;
  frac_lo = (mH_hi - mH)/delta;
  frac_hi = (mH - mH_lo)/delta;
}
//----------------------------------------------------------------------------------------
void setZzntVals(CtrlFlags &flags, varkeepter *flatoutput, filestuff *fs, TString cls, TMVA::Reader *reader, Mela *mela, TMVA::Reader *spinReader)
{
  flatoutput->setval("run",             (UInt_t)fs->info->run);
  flatoutput->setval("evt",             (UInt_t)fs->info->evt);
  flatoutput->setval("channel",         (UInt_t)fs->kine->channel);
  flatoutput->setval("ZZpt",            (Double_t)fs->kine->ZZpt);
  flatoutput->setval("m4l",             (Double_t)fs->kine->m4l);
  flatoutput->setval("w",               (Double_t)fs->weights->w);
  flatoutput->setval("won",             (Double_t)fs->weights->won);
  flatoutput->setval("woff",            (Double_t)fs->weights->woff);
  flatoutput->setval("npuw",            (Double_t)fs->weights->npuw);

  // set pt4l up/down weights
  Double_t pt4lDown(1),pt4lUp(1),pt4lPdfDown(1),pt4lPdfUp(1),pt4lHresDown(1),pt4lHresUp(1);
  if(cls=="qqZZ" || (cls=="gfH" && !fs->fname_.Contains(TRegexp("[02][mp]")))) {
    assert(fs->get_zz_chain()->FindBranch("geninfo"));
    assert(wgtHists["renfac_Lo"]);
    assert(wgtHists["renfac_Hi"]);
    assert(wgtHists["pdf_Lo"]);
    assert(wgtHists["pdf_Hi"]);
    pt4lDown = getWeight(wgtHists["renfac_Lo"], fs->gen->pt_zz);
    pt4lUp   = getWeight(wgtHists["renfac_Hi"], fs->gen->pt_zz);
    pt4lPdfDown = getWeight(wgtHists["pdf_Lo"],    fs->gen->pt_zz);
    pt4lPdfUp   = getWeight(wgtHists["pdf_Hi"],    fs->gen->pt_zz);
    // cout << "pt4l: " << fs->gen->pt_zz << " renfac: " << pt4lDown << " " << pt4lUp << " pdf: " << pt4lPdfDown << " " << pt4lUp;
    if(cls != "qqZZ" ) {
      assert(wgtHists["hres_Lo"]);
      assert(wgtHists["hres_Hi"]);
      pt4lHresDown = getWeight(wgtHists["hres_Lo"],   fs->gen->pt_zz);
      pt4lHresUp   = getWeight(wgtHists["hres_Hi"],   fs->gen->pt_zz);
      // cout << " hres: " << pt4lHresDown << " " << pt4lHresUp;
    }
    // cout << endl;
  }
  flatoutput->setval("pt4lDown",        (Double_t)pt4lDown);
  flatoutput->setval("pt4lUp",          (Double_t)pt4lUp);
  flatoutput->setval("pt4lPdfDown",     (Double_t)pt4lPdfDown);
  flatoutput->setval("pt4lPdfUp",       (Double_t)pt4lPdfUp);
  flatoutput->setval("pt4lHresDown",    (Double_t)pt4lHresDown);
  flatoutput->setval("pt4lHresUp",      (Double_t)pt4lHresUp);

  // set mH weights
  Double_t wmH(1);
    if((cls=="gfH" && !fs->fname_.Contains(TRegexp("[02][mp]"))) || cls=="vbfH") {
    assert(fs->get_zz_chain()->FindBranch("geninfo"));
    GenInfoStruct gen(*fs->gen);
    TLorentzVector genZ1,genZ2;
    genZ1.SetPtEtaPhiM( gen.vpt_a, gen.veta_a, gen.vphi_a, gen.vmass_a);
    genZ2.SetPtEtaPhiM( gen.vpt_b, gen.veta_b, gen.vphi_b, gen.vmass_b);
    double genM4l = (genZ1 + genZ2).M();
    wmH = wrap->getweight(flags.mH, flags.hwidth, 172.5, genM4l, 1);
    // cout << "m4l weights: " << genM4l << " " << flatoutput->getval("wmH") << endl;
  }
  flatoutput->setval("wmH",             (Double_t)wmH);

  if(flags.makeMela) {
    assert(mela);
    TLorentzVector Z1_lept1, Z1_lept2, Z2_lept1, Z2_lept2;
    KinematicsStruct kine(*(fs->kine));

    Z1_lept1.SetPtEtaPhiM( kine.l1pt, kine.l1eta, kine.l1phi, abs(kine.l1type)==13 ? 105.658369e-3 : 0.51099892e-3);
    Z1_lept2.SetPtEtaPhiM( kine.l2pt, kine.l2eta, kine.l2phi, abs(kine.l2type)==13 ? 105.658369e-3 : 0.51099892e-3);
    Z2_lept1.SetPtEtaPhiM( kine.l3pt, kine.l3eta, kine.l3phi, abs(kine.l3type)==13 ? 105.658369e-3 : 0.51099892e-3);
    Z2_lept2.SetPtEtaPhiM( kine.l4pt, kine.l4eta, kine.l4phi, abs(kine.l4type)==13 ? 105.658369e-3 : 0.51099892e-3);

    // NOTE: it's p.d.g. ID, so opposite sign to the charge
    int Z1_lept1Id(-1*kine.l1q*kine.l1type),Z1_lept2Id(-1*kine.l2q*kine.l2type),Z2_lept1Id(-1*kine.l3q*kine.l3type),Z2_lept2Id(-1*kine.l4q*kine.l4type);
    float costhetastar,costheta1,costheta2,phi,phistar1,kd,psig,pbkg,spinKd;
    mela->computeKD(Z1_lept1, Z1_lept1Id, Z1_lept2, Z1_lept2Id, Z2_lept1, Z2_lept1Id, Z2_lept2, Z2_lept2Id,
                    costhetastar, costheta1, costheta2, phi, phistar1, kd, psig, pbkg, flags.withPtY, flags.withPtY);
  
    flatoutput->setval("melaLD",                (Float_t)kd);

    if(flags.multisigs) {
      assert(psMela);
      psMela->computeKD(Z1_lept1, Z1_lept1Id, Z1_lept2, Z1_lept2Id, Z2_lept1, Z2_lept1Id, Z2_lept2, Z2_lept2Id,
                        spinKd, psig, pbkg);
      flatoutput->setval("SpinBDT",     (Float_t)spinKd);
    }
  } else {
    assert(reader);
    if(flags.multiclass) {
      flatoutput->setval("BDT",         (Float_t)((reader->EvaluateMulticlass("BDTG"))[0])); // signal component
      flatoutput->setval("bkgBDT",      (Float_t)((reader->EvaluateMulticlass("BDTG"))[1])); // bkg (ZZ) component
      flatoutput->setval("bkg2BDT",     (Float_t)((reader->EvaluateMulticlass("BDTG"))[2])); // bkg2 (Z+jets) component
    } else {
      flatoutput->setval("BDT",         (Float_t)(reader->EvaluateMVA("BDTG")));
    }
    if(flags.multisigs) {
      assert(spinReader);
      flatoutput->setval("SpinBDT",     (Float_t)spinReader->EvaluateMVA("BDTG"));
    }
  }
}
//----------------------------------------------------------------------------------------
double getWeight(TH1F *hist, float val)
{
  int bin = hist->FindBin(val);
  // return 1 for out-of-bounds
  if(bin<1 || bin>hist->GetNbinsX())
    return 1;
  return hist->GetBinContent(bin);
}
//----------------------------------------------------------------------------------------
TString findClosestPtWgts(TFile *wgtFile, int mH, TString scaleType) // scaleType is renfac, pdf, or hres
{
  Long_t closestMass(-1);
  TString closestMassStr;
  TList *keys = wgtFile->GetListOfKeys();
  for(unsigned ikey=0; ikey<keys->GetEntries(); ikey++) {
    TString name(keys->At(ikey)->GetName());
    TString tmpName(name(TRegexp("h[0-9][0-9]*zz4l_"+scaleType)));
    TString massStr(tmpName(TRegexp("[0-9][0-9]*")));
    int mass = atoi(massStr);
    if(mass<100 || mass > 1200) continue;
    if(abs(mass - mH) < abs(closestMass - mH)) {
      closestMass = mass;
      closestMassStr = massStr;
    }
  }

  cout << "  using closest mass of: " << closestMassStr << endl;
  return closestMassStr;
}
//----------------------------------------------------------------------------------------
void findHist(TString cls, Long_t mH, TString fname, TFile *wgtFile, TString scaleType, TString hiLo)
{
  assert(hiLo=="Lo" || hiLo=="Hi");
  TString procStr;
  if(cls=="gfH") {
    procStr = TString("h")+mH+"zz4l";
  } else if(cls=="qqZZ") {
    if(fname.Contains("zz4e")) procStr = "zz4e";
    if(fname.Contains("zz4m")) procStr = "zz4m";
    if(fname.Contains("zz2e2m")) procStr = "zz2e2m";
  } else assert(0);
  assert(procStr!="");

  TString histName(hiLo+"Wgts_"+procStr+"_"+scaleType);
  TH1F *htest = (TH1F*)wgtFile->Get(histName);
  TString closestMass;
  if(!htest) {
    cout << "  hist name " << histName << " not found" << endl;
    if(cls=="gfH") {
      closestMass = findClosestPtWgts(wgtFile, mH, scaleType);
      histName.ReplaceAll(TString("")+mH, closestMass);
    } else {
      histName.ReplaceAll("zz4e","zz2e2m"); // only made 2e2m for a few of them
      histName.ReplaceAll("zz4m","zz2e2m");
    }
    cout << "  trying with " << histName << endl;
  }

  wgtHists[scaleType+"_"+hiLo] = (TH1F*)wgtFile->Get(histName); assert(wgtHists[scaleType+"_"+hiLo]);
}
//----------------------------------------------------------------------------------------
void setPtWgtHists(TString cls, TString fname, Long_t mH, TFile *wgtFile)
{
  assert(wgtFile->IsOpen());

  findHist(cls, mH, fname, wgtFile, "renfac",   "Lo");
  findHist(cls, mH, fname, wgtFile, "renfac",   "Hi");
  findHist(cls, mH, fname, wgtFile, "pdf",      "Lo");
  findHist(cls, mH, fname, wgtFile, "pdf",      "Hi");
  if(cls=="gfH") {
    findHist(cls, mH, fname, wgtFile, "hres",    "Lo");
    findHist(cls, mH, fname, wgtFile, "hres",    "Hi");
  }
}
Revision:	1.2
Committed:	Mon Dec 17 12:34:49 2012 UTC (12 years, 5 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+147 -75 lines
Log Message:	* empty log message *