rootEWKanalyzer/src/baseClass.C

#define baseClass_cxx
#include "baseClass.h"

#include "TH1I.h"
#include "TFile.h"

baseClass::baseClass(string * inputList, string * cutFile, string * treeName, string * outputFileName, string * cutEfficFile)
{
  std::cout << "baseClass::baseClass(): begins " << std::endl;
  inputList_ = inputList;
  cutFile_ = cutFile;
  treeName_= treeName;
  outputFileName_ = outputFileName;
  cutEfficFile_ = cutEfficFile;
  init();
  std::cout << "baseClass::baseClass(): ends " << std::endl;
}

baseClass::~baseClass()
{
  std::cout << "baseClass::~baseClass(): begin " << std::endl;
  // the following lines are from the jet prompt thing...
//   if( !writeCutHistos() )
//     {
//       STDOUT("ERROR: writeCutHistos did not complete successfully.");
//     }
//   if( !writeCutEfficFile() )
//     {
//       STDOUT("ERROR: writeStatFile did not complete successfully.");
//     }
  output_root_->Close();
  std::cout << "baseClass::~baseClass(): end " << std::endl;
}


double baseClass::multiply(double a, double b)
{
  double c;
  c = a * b;
  return c;
}

void baseClass::getHLTtable() {
  bool printTriggerTable=true;
  if(tree_ == NULL){
    std::cout << "baseClass::init(): ERROR: tree_ = NULL " << std::endl;
    return;
  }

  tree_->GetEntry();
  TFile *f = tree_->GetCurrentFile();

  cout<<f->GetName()<<endl;;
  TH1I *hlt_stats = (TH1I*)f->Get("analyze/HLTTriggerStats");
  cout<<"here 3"<<endl;

  // clear the map
  if(!HLTmap.empty()) {
    HLTmap.erase(HLTmap.begin(),HLTmap.end());
  }

  for( int i=0; i < hlt_stats->GetNbinsX(); i++ ){
    if(printTriggerTable==true) cout<<i<<": "<<hlt_stats->GetXaxis()->GetBinLabel(i+1)<<endl;
    HLTmap[hlt_stats->GetXaxis()->GetBinLabel(i+1)] = i;
  }
}

int baseClass::getHLTtriggerBit(string triggerName) {
  if(HLTmap.empty()) {
    cerr<<"HLTtrigger: no HLTmap available!"<<endl;
    return -1;
  }
  map<string, int>::iterator iter = HLTmap.find(triggerName);
  //cout<<"HLT bit: "<<iter->second<<endl;  // iter->second is the HLT trigger bit
  if(iter==HLTmap.end()) {
    cerr<<"HLTtrigger: trigger name "<<triggerName<< " not found!"<<endl;
    return -1;
  }
  else {
    return iter->second;
  }
}

bool baseClass::triggerBitSelection(int hlt[200], int l1phys[128], int l1tech[64])
{
  bool passed=false;
  if((l1tech[40]==1 || l1tech[41]==1) &&  // pass BSC
     (l1tech[36]==0 && l1tech[37]==0 && l1tech[38]==0 && l1tech[39]==0)// && // Beam Halo Veto
//      l1tech[0]==1 && // pass BPTX
//      hlt[116]==1 // pass Physics Declared Bit
     )
    {
      passed=true;
    }
  return passed;
  
}

bool baseClass::goodPrimaryVertexSelection(int ndof, double z)
{
  bool passed=false;
  if(ndof>=5 && fabs(z)<15) passed=true;
  return passed;
}

// bool baseClass::noBeamScrapingSelection(int nTracks, int goodTrack)
// {
//   passed=false;
//   if(nTracks<10) passed=true;
//   else if(ntracks>=10) {
//     for(int i=0;i<nTracks;i++) {
//       if()
//     }
//   }
//   return passed;
// }

bool myEvent::WorkingPointID(myEvent electron, vector<double> cuts)
{
  bool passed=false;
  double etaMax_B=1.44;
  double etaMax_E=2.6;
  // get the variables
  double See=electron.sigmaIetaIeta[0];
  double Dphi=fabs(electron.deltaPhi[0]);
  double Deta=fabs(electron.deltaEta[0]);
  double HoE=electron.HoverE[0];
  double TrkIso=electron.trkIso[0];
  double EcalIso=electron.ecalIso[0];
  double HcalIso=electron.hcalIso[0];
  double missingHits=electron.numberOfMissingInnerHits[0];
  // get the cuts
  double See_B=cuts[0];
  double Dphi_B=cuts[1];
  double Deta_B=cuts[2];
  double HoE_B=cuts[3];
  double TrkIso_B=cuts[4];
  double EcalIso_B=cuts[5];
  double HcalIso_B=cuts[6];
  double missingHits_B=cuts[7];
  double See_E=cuts[8];
  double Dphi_E=cuts[9];
  double Deta_E=cuts[10];
  double HoE_E=cuts[11];
  double TrkIso_E=cuts[12];
  double EcalIso_E=cuts[13];
  double HcalIso_E=cuts[14];
  double missingHits_E=cuts[15];
  if(fabs(electron.fourvector[0].Eta())<etaMax_B) {
    if(See<See_B &&
       Dphi<Dphi_B &&
       Deta<Deta_B &&
       HoE<HoE_B &&
       TrkIso<TrkIso_B &&
       EcalIso<EcalIso_B &&
       HcalIso<HcalIso_B &&
       missingHits<missingHits_B) passed=true;
  }
  else if(fabs(electron.fourvector[0].Eta())>etaMax_B && fabs(electron.fourvector[0].Eta())<etaMax_E) {
    if(See<See_E &&
       Dphi<Dphi_E &&
       Deta<Deta_E &&
       HoE<HoE_E &&
       TrkIso<TrkIso_E &&
       EcalIso<EcalIso_E &&
       HcalIso<HcalIso_E &&
       missingHits<missingHits_E) passed=true;
  }
  return passed;
}
  
bool myEvent::WorkingPointNminus1(myEvent electron, vector<double> cuts, char *cut)
{
  bool passed=false;
  int combination[8]={0,0,0,0,0,0,0,0};
  double etaMax_B=1.44;
  double etaMax_E=2.6;
  // get the variables
  double See=electron.sigmaIetaIeta[0];
  double Dphi=fabs(electron.deltaPhi[0]);
  double Deta=fabs(electron.deltaEta[0]);
  double HoE=electron.HoverE[0];
  double TrkIso=electron.trkIso[0];
  double EcalIso=electron.ecalIso[0];
  double HcalIso=electron.hcalIso[0];
  double missingHits=electron.numberOfMissingInnerHits[0];


  // get the cuts
  double See_B=cuts[0];
  double Dphi_B=cuts[1];
  double Deta_B=cuts[2];
  double HoE_B=cuts[3];
  double TrkIso_B=cuts[4];
  double EcalIso_B=cuts[5];
  double HcalIso_B=cuts[6];
  double CombIso_B=cuts[7];
  double missingHits_B=cuts[16];
  double See_E=cuts[8];
  double Dphi_E=cuts[9];
  double Deta_E=cuts[10];
  double HoE_E=cuts[11];
  double TrkIso_E=cuts[12];
  double EcalIso_E=cuts[13];
  double HcalIso_E=cuts[14];
  double CombIso_E=cuts[15];
  double missingHits_E=cuts[16];
  double dist=cuts[17];


  // get the cut which should NOT be made
  if(!strcmp(cut,"sigmaIetaIeta")) combination[0]=1;
  else if(!strcmp(cut,"deltaPhi")) combination[1]=1;
  else if(!strcmp(cut,"deltaEta")) combination[2]=1;
  else if(!strcmp(cut,"HoverE")) combination[3]=1;
  else if(!strcmp(cut,"trkIso")) combination[4]=1;
  else if(!strcmp(cut,"ecalIso")) combination[5]=1;
  else if(!strcmp(cut,"hcalIso")) combination[6]=1;
  else if(!strcmp(cut,"conversions")) combination[7]=1;
  else cout<<"ERROR: not able to give N-1 flag for the working point!"<<endl;
  // divide in barrel and endcaps
  if(fabs(electron.fourvector[0].Eta())<etaMax_B) {
    // sigma ieta ieta
    if(combination[0]==1 &&
       Dphi<Dphi_B &&
       Deta<Deta_B &&
       HoE<HoE_B &&
       TrkIso<TrkIso_B &&
       EcalIso<EcalIso_B &&
       HcalIso<HcalIso_B &&
       missingHits<=missingHits_B) passed=true;
    // delta phi
    else if(combination[1]==1 &&
            See<See_B &&
            Deta<Deta_B &&
            HoE<HoE_B &&
            TrkIso<TrkIso_B &&
            EcalIso<EcalIso_B &&
            HcalIso<HcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // delta eta
    else if(combination[2]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            HoE<HoE_B &&
            TrkIso<TrkIso_B &&
            EcalIso<EcalIso_B &&
            HcalIso<HcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // H over E
    else if(combination[3]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            Deta<Deta_B &&
            TrkIso<TrkIso_B &&
            EcalIso<EcalIso_B &&
            HcalIso<HcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // Track isolation 
    else if(combination[4]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            Deta<Deta_B &&
            HoE<HoE_B &&
            EcalIso<EcalIso_B &&
            HcalIso<HcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // ECAL isolation 
    else if(combination[5]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            Deta<Deta_B &&
            HoE<HoE_B &&
            TrkIso<TrkIso_B &&
            HcalIso<HcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // HCAL isolation 
    else if(combination[6]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            Deta<Deta_B &&
            HoE<HoE_B &&
            TrkIso<TrkIso_B &&
            EcalIso<EcalIso_B &&
            missingHits<=missingHits_B) passed=true;
    // Conversions
    else if(combination[7]==1 &&
            See<See_B &&
            Dphi<Dphi_B &&
            Deta<Deta_B &&
            HoE<HoE_B &&
            TrkIso<TrkIso_B &&
            EcalIso<EcalIso_B &&
            HcalIso<HcalIso_B) passed=true;


  } // barrel
  else if(fabs(electron.fourvector[0].Eta())>etaMax_B && fabs(electron.fourvector[0].Eta())<etaMax_E) {
    // sigma ieta ieta
    if(combination[0]==1 &&
       Dphi<Dphi_E &&
       Deta<Deta_E &&
       HoE<HoE_E &&
       TrkIso<TrkIso_E &&
       EcalIso<EcalIso_E &&
       HcalIso<HcalIso_E &&
       missingHits<=missingHits_E) passed=true;
    // delta phi
    else if(combination[1]==1 &&
            See<See_E &&
            Deta<Deta_E &&
            HoE<HoE_E &&
            TrkIso<TrkIso_E &&
            EcalIso<EcalIso_E &&
            HcalIso<HcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // delta eta
    else if(combination[2]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            HoE<HoE_E &&
            TrkIso<TrkIso_E &&
            EcalIso<EcalIso_E &&
            HcalIso<HcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // H over E
    else if(combination[3]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            Deta<Deta_E &&
            TrkIso<TrkIso_E &&
            EcalIso<EcalIso_E &&
            HcalIso<HcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // Track isolation 
    else if(combination[4]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            Deta<Deta_E &&
            HoE<HoE_E &&
            EcalIso<EcalIso_E &&
            HcalIso<HcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // ECAL isolation 
    else if(combination[5]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            Deta<Deta_E &&
            HoE<HoE_E &&
            TrkIso<TrkIso_E &&
            HcalIso<HcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // HCAL isolation 
    else if(combination[6]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            Deta<Deta_E &&
            HoE<HoE_E &&
            TrkIso<TrkIso_E &&
            EcalIso<EcalIso_E &&
            missingHits<=missingHits_E) passed=true;
    // Conversions
    else if(combination[7]==1 &&
            See<See_E &&
            Dphi<Dphi_E &&
            Deta<Deta_E &&
            HoE<HoE_E &&
            TrkIso<TrkIso_E &&
            EcalIso<EcalIso_E &&
            HcalIso<HcalIso_E) passed=true;


  } // endcaps
  return passed;
}


testEvent baseClass::fillTestEvent() {
  // fill the event struct
  testEvent event;
  for(int i=0; i<NEles;i++) {
    TLorentzVector vec;
    vec.SetPxPyPzE(ElPx[i],ElPy[i],ElPz[i],ElE[i]);
    event.fourMomentum.push_back(vec);
  }
  return event;
}

int baseClass::ElectronPreselection(int i) {
  // function to make electron preselection
  // 0: preselected electron
  // 1: not accepted in eta
  // 2: not accepted in pt
  // 3: not accepted in H/E
  // 4: not accepted in dphi
  // 5: not accepted in deta

  //double etaMax_B=1.44;
  double etaMax_E=2.5;
  double ptMin=20;
  double HoverEMin=0.15;
  double deltaPhiMax=0.15;
  double deltaEtaMax=0.02;

  bool cutOn_Eta = true;
  bool cutOn_Pt  = true;
  bool cutOn_HoE = true;
  bool cutOn_dPhi= true;
  bool cutOn_dEta= true;

  double eta  = ElSCEta[i];
//   double eta  = ElEta[i];

//   double dPhiCorr=dphiCorrections(ElEta[i], ElPhi[i]);
//   double dEtaCorr=detaCorrections(ElEta[i], ElPhi[i]);
  double dPhiCorr=dphiCorrections(eta, ElPhi[i]);
  double dEtaCorr=detaCorrections(eta, ElPhi[i]);

  TVector3 track;
  track.SetXYZ(ElPx[i], ElPy[i], ElPz[i]);
  double trackEta=track.Eta();
//   double pt=ElCaloEnergy[i]/cosh(trackEta);
  double pt   = ElPt[i];

  double hoe  = ElHcalOverEcal[i];
  double dphi = ElDeltaPhiSuperClusterAtVtx[i] - dPhiCorr;
  double deta = ElDeltaEtaSuperClusterAtVtx[i] - dEtaCorr;
//   double dphi = ElDeltaPhiSuperClusterAtVtx[i];
//   double deta = ElDeltaEtaSuperClusterAtVtx[i];

  if(cutOn_Eta)  if(fabs(eta) > etaMax_E)  return 1;
  if(cutOn_Pt)   if(pt  < ptMin)     return 2;
  if(cutOn_HoE)  if(hoe > HoverEMin) return 3;
  if(cutOn_dPhi) if(fabs(dphi) > deltaPhiMax) return 4;
  if(cutOn_dEta) if(fabs(deta) > deltaEtaMax) return 5;

  return 0;
}

int baseClass::WPElectronID(int i, NminusOneCutLabel c, vector<double> cuts) {
  // function to ID and isolate electrons
  // return values (int):
  // 0 totally IDed & isolated electron
  // 1 See cut not passed
  // 2 dEta cut not passed
  // 3 dPhi
  // 4 HoE
  // 5 combIso
  // 6 conversion
  
  double etaMax_B=1.44;
  double etaMin_E=1.56;
  double etaMax_E=2.5;
  double SeeMax_B=cuts[0];
  double dPhiMax_B=cuts[1];
  double dEtaMax_B=cuts[2];
  double HoEMax_B=cuts[3];
  double TrkIsoMax_B=cuts[4];
  double EcalIsoMax_B=cuts[5];
  double HcalIsoMax_B=cuts[6];
  double combIsoMax_B=cuts[7];
  double SeeMax_E=cuts[8];
  double dPhiMax_E=cuts[9];
  double dEtaMax_E=cuts[10];
  double HoEMax_E=cuts[11];
  double TrkIsoMax_E=cuts[12];
  double EcalIsoMax_E=cuts[13];
  double HcalIsoMax_E=cuts[14];
  double combIsoMax_E=cuts[15];
  int hits=cuts[16];
  double dist=cuts[17];
  double cot=dist;

  bool cutOn_See=true;
  bool cutOn_dPhi=true;
  bool cutOn_dEta=true;
  bool cutOn_HoE=true;
  bool cutOn_combIso=true;
  bool cutOn_conversion=true;

  // check geometrical acceptance
  // for the isolation use the SC eta instead
  // of the electron eta
  double eta = ElSCEta[i];
//   double eta = ElEta[i];

  double dPhiCorr=dphiCorrections(eta, ElPhi[i]);
  double dEtaCorr=detaCorrections(eta, ElPhi[i]);

  double See  = ElSigmaIetaIeta[i];
  // apply the correction for the endcaps
  double dPhi = fabs(ElDeltaPhiSuperClusterAtVtx[i]-dPhiCorr);
  double dEta = fabs(ElDeltaEtaSuperClusterAtVtx[i]-dEtaCorr);
  double HoE  = ElHcalOverEcal[i];
  double trkIso  = ElDR03TkSumPt[i];
  double ecalIso = ElDR03EcalRecHitSumEt[i];
  double hcalIso = ElDR03HcalTowerSumEt[i];
  double combIso_B = (trkIso + max(0., ecalIso - 1.) + hcalIso) / ElPt[i];
  double combIso_E = (trkIso + ecalIso + hcalIso) / ElPt[i];
  bool isConversion;
  if((fabs(ElConvPartnerTrkDist[i])<dist && fabs(ElConvPartnerTrkDCot[i])<cot) || ElNumberOfMissingInnerHits[i]>hits) isConversion=true;
  else isConversion=false;

  if(c==see)  cutOn_See  = false;
  if(c==dphi) cutOn_dPhi = false;
  if(c==deta) cutOn_dEta = false;
  if(c==hoe)  cutOn_HoE  = false;
  if(c==combiso)    cutOn_combIso    = false;
  if(c==conversion) cutOn_conversion = false;


//   if(Event==27837545 || Event==73939173 || Event==15512525 || Event==191815699 || Event==29377895 || Event==29377895 || Event==74868952 || Event==76643965 || Event==62016483 || Event==4352136) {
//     cout<<eta<<endl;
//     cout<<combIso_B<<"/"<<combIso_E<<endl;
//     cout<<combIsoMax_B<<"/"<<combIsoMax_E<<endl;
//     cout<<"See: "<<See<<endl;
//     cout<<"dPhi: "<<dPhi<<endl;
//     cout<<"dEta: "<<dEta<<endl;
//     cout<<"HoE: "<<HoE<<endl;
//     cout<<"pt: "<<ElPt[i]<<endl;
//   }

  if(fabs(eta)<=etaMax_B) {
    if(cutOn_See)  if(See  > SeeMax_B)  return 1;
    if(cutOn_dPhi) if(dPhi > dPhiMax_B) return 2;
    if(cutOn_dEta) if(dEta > dEtaMax_B) return 3;
    if(cutOn_HoE)  if(HoE  > HoEMax_B)  return 4;
    if(cutOn_combIso) if(combIso_B > combIsoMax_B) return 5;
    if(cutOn_conversion) if(isConversion==true) return 6;
    return 0;
  }
  if(fabs(eta)>etaMin_E && fabs(eta)<etaMax_E) {
    if(cutOn_See)  if(See  > SeeMax_E)  return 1;
    if(cutOn_dPhi) if(dPhi > dPhiMax_E) return 2;
    if(cutOn_dEta) if(dEta > dEtaMax_E) return 3;
    if(cutOn_HoE)  if(HoE  > HoEMax_E)  return 4;
    if(cutOn_combIso) if(combIso_E > combIsoMax_E) return 5;
    if(cutOn_conversion) if(isConversion==true) return 6;
    return 0;
  }
  else {
    cout<<"WPelectronID: electron not in ECAL acceptance region!!"<<endl;
    return -1;
  }
}

double baseClass::detaCorrections(double etaEle,double phiEle) {
  // corrects for misalignement in the endcpas
  TF2 f12_fit2("f12_fit2","[0]+(TMath::TanH(y)/325.)*([1]-([2]*TMath::SinH(y)*TMath::Cos(x-[3])))",-TMath::Pi(),TMath::Pi(),-10.,10.);
  if (etaEle>1.479)
    {
      f12_fit2.SetParameter(0,0.0013);
      f12_fit2.SetParameter(1,-0.06);
      f12_fit2.SetParameter(2,0.52);
      f12_fit2.SetParameter(3,2.17);
    }
  else if (etaEle<-1.479)
    {
      f12_fit2.SetParameter(0,-0.0013);
      f12_fit2.SetParameter(1,-0.32);
      f12_fit2.SetParameter(2,0.45);
      f12_fit2.SetParameter(3,-1.58);
    }
  return f12_fit2.Eval(phiEle,etaEle);
}

double baseClass::dphiCorrections(double etaEle,double phiEle) {
  // corrects for misalignement in the endcpas
  TF2 f12_fit2("f12_fit2","[0]+[1]*(TMath::SinH(y)/325.)*(TMath::Sin([2]-x))",-TMath::Pi(),TMath::Pi(),-10.,10.);
  if (etaEle>1.479)
    {
      f12_fit2.SetParameter(0,0.0);
      f12_fit2.SetParameter(1,0.52);
      f12_fit2.SetParameter(2,2.17);
    }
  else if (etaEle<-1.479)
    {
      f12_fit2.SetParameter(0,0.0);
      f12_fit2.SetParameter(1,0.45);
      f12_fit2.SetParameter(2,-1.58);
    }
   return f12_fit2.Eval(phiEle,etaEle);
}


int baseClass::WSelection(int i, WNminusOneCutLabel c) {
  // idea: write cuts into a vector and give the vector to the function
//   double etaMax_B=1.44;
  double pre_metMin=20;
  double pre_ptMin=20;
  double metMin=30;
  double ptMin=25;
  double mtMin=60;

  bool cutOn_Preselection=true;
  bool cutOn_Met=true;
  bool cutOn_Pt=true;
  bool cutOn_Mt=true;

  double met=PFMET;
  double pt=ElPt[i];
  double delPhiMet=deltaPhi(ElPhi[i], PFMETphi);
  double mt=transverseMass(pt, met, delPhiMet);
  // check geometrical acceptance
  // for the isolation use the SC eta instead
  // of the electron eta
  //   double eta = ElSCEta[i];
//   double eta = ElEta[i];
  
  if(c==pre) cutOn_Preselection  = false;
  if(c==misset) cutOn_Met = false;
  if(c==pt)  cutOn_Pt = false;
  if(c==mt)  cutOn_Mt = false;


//   if(fabs(eta)<=etaMax_B) {
  if(cutOn_Preselection) if(met<pre_metMin || pt<pre_ptMin) return 1;
  if(cutOn_Met) if(met<metMin) return 2;
  if(cutOn_Pt) if(pt<ptMin) return 3;
  if(cutOn_Mt) if(mt<mtMin) return 4;
  return 0;
//   }
//   if(fabs(eta)>etaMax_B) {
//     if(cutOn_Preselection) if(met<pre_metMin || pt<pre_ptMin) return 1;
//     if(cutOn_Met) if(met<metMin) return 2;
//     if(cutOn_Pt) if(pt<ptMin) return 3;
//     if(cutOn_Mt) if(mt<mtMin) return 4;
//     return 0;
//   }
}


// void initializeWP(vector<double>cuts) {
//   SeeMax_B=cuts[0];
//   DphiMax_B=cuts[1];
//   DetaMax_B=cuts[2];
//   HoEMax_B=cuts[3];
//   combIsoMax_B=cuts[4];
// //   TrkIsoMax_B=cuts[4];
// //   EcalIsoMax_B=cuts[5];
// //   HcalIsoMax_B=cuts[6];
//   SeeMax_E=cuts[8];
//   DphiMax_E=cuts[9];
//   DetaMax_E=cuts[10];
//   HoEMax_E=cuts[11];
// //   TrkIsoMax_E=cuts[12];
// //   EcalIsoMax_E=cuts[13];
// //   HcalIsoMax_E=cuts[14];
//   combIsoMax_E=cuts[12];

// }

bool baseClass::ethNminus1(int isECALdriven, double deltaPhi, double deltaEta, double e, double p, double eta, char *cut, bool iso)
{
  bool passed=false;
  double oneOverEoneOverP=1/e-1/p;
  // cut values
  double etaMax_B=1.44;
  double etaMax_E=2.6;
  double deltaPhiCut_B=0.02;
  double deltaPhiCut_E=0.02;
  double deltaEtaCut_B=0.004;
  double deltaEtaCut_E=0.006;
  double oneOverEoneOverPCut_B=0.005;
  double oneOverEoneOverPCut_E=0.007;
  int combination [4] = {0,0,0,0};
  if(!strcmp(cut,"deltaPhi")) combination [0] = 1;
  else if(!strcmp(cut,"deltaEta")) combination [1] =1;
  else if(!strcmp(cut,"oneOverEminusOneOverP")) combination [2] = 1;
  else if(!strcmp(cut,"iso")) combination[3] = 1;
  else {
    cout<<"ERROR: not able to give N-1 flag for the ETH!"<<endl;
  }
  if(isECALdriven == 1) {
    // barrel
    if(fabs(eta)<etaMax_B) {
      // delta phi
      if(combination[0]==1 && fabs(deltaEta)<deltaEtaCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B && iso==true) {
        passed = true;
      } 
      // delta eta
      else if(combination[1]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B && iso==true) {
        passed = true;
      } 


      //  |1/E-1/p|
      else if(combination[2]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(deltaEta)< deltaEtaCut_B && iso==true) {
        passed = true;
      }
      // iso
      else if((combination[3]==1 && fabs(deltaPhi)<deltaPhiCut_B && fabs(deltaEta)< deltaEtaCut_B && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B)) {
        passed=true;
      } 
    }
    //endcap
    else if(fabs(eta)>etaMax_B && fabs(eta)<etaMax_E) {
      // delta phi
      if(combination[0]==1 && fabs(deltaEta)< deltaEtaCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E && iso==true) {
        passed = true;
      }
      // delta eta 
      else if(combination[1]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E && iso==true) {
        passed = true;
      }
      // |1/E-1/p|      
      else if(combination[2]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(deltaEta)< deltaEtaCut_E && iso ==true) {
        passed = true;
      }
      // iso
      else if((combination[3]==1 && fabs(deltaPhi)< deltaPhiCut_E && fabs(deltaEta)< deltaEtaCut_E && fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E)) {
        passed=true;
      } 
    } // endcap
  } // ecal driven 
  return passed;
}


bool baseClass::ethID(int isECALdriven, double deltaPhi, double deltaEta, double e, double p, double eta)
{
  bool passed=false;
  double oneOverEoneOverP=1/e-1/p;
  // cut values
  double etaMax_B=1.44;
  double etaMax_E=2.6;
  double deltaPhiCut_B=0.02;
  double deltaPhiCut_E=0.02;
  double deltaEtaCut_B=0.004;
  double deltaEtaCut_E=0.006;
  double oneOverEoneOverPCut_B=0.005;
  double oneOverEoneOverPCut_E=0.007;
  if(isECALdriven == 1 &&
     fabs(eta)<etaMax_B &&
     fabs(deltaPhi)< deltaPhiCut_B &&
     fabs(deltaEta)< deltaEtaCut_B &&
     fabs(oneOverEoneOverP)<oneOverEoneOverPCut_B) {
    passed=true;
  }
  else if(isECALdriven == 1 &&
     fabs(eta)>etaMax_B &&
     fabs(eta)<etaMax_E &&
     fabs(deltaPhi)< deltaPhiCut_E &&
     fabs(deltaEta)< deltaEtaCut_E &&
     fabs(oneOverEoneOverP)<oneOverEoneOverPCut_E) {
    passed=true;
  }  
  return passed;
}

bool baseClass::isolation(double sumpt, double e) 
{
  bool passed =false;
  double iso;
  iso = sumpt/e;
  if(iso<0.1) {
    passed =true;
  }
  return passed;
}

double baseClass::deltaPhi(double phi1, double phi2) {
  double delPhi;
  double Pi=2*acos(0);
  delPhi=fabs(phi1-phi2);
  if(delPhi>Pi) {
    delPhi=2*Pi-delPhi;
  }
  return delPhi;
}

double baseClass::transverseMass(double pt1, double pt2, double delPhi) {
  // pt2=met!!
  double mt;
  mt=sqrt(2*fabs(pt1)*fabs(pt2)*(1-cos(delPhi)));
  return mt;
}

void baseClass::init()
{
  STDOUT("begin");
    //std::cout << "baseClass::init(): begin " << std::endl;

  tree_ = NULL;
  readInputList();
//   readCutFile();   // not using cut file, need dummy cut file as input...
  if(tree_ == NULL){
    std::cout << "baseClass::init(): ERROR: tree_ = NULL " << std::endl;
    return;
  }
  Init(tree_);

  //char output_root_title[200];
  //sprintf(output_root_title,"%s%s",&std::string(*outputFileName_)[0],".root");
  //output_root_ = new TFile(&output_root_title[0],"RECREATE");

  //directly from string
  output_root_ = new TFile((*outputFileName_ + ".root").c_str(),"RECREATE");
  
  //  for (map<string, cut>::iterator it = cutName_cut_.begin(); 
  //   it != cutName_cut_.end(); it++) STDOUT("cutName_cut->first = "<<it->first)
  //  for (vector<string>::iterator it = orderedCutNames_.begin(); 
  //       it != orderedCutNames_.end(); it++) STDOUT("orderedCutNames_ = "<<*it)
  STDOUT("end");
}

void baseClass::readInputList()
{

  TChain *chain = new TChain(treeName_->c_str());
  char pName[500];
  skimWasMade_ = true;
  NBeforeSkim_ = 0;
  int NBeforeSkim;

  std::cout << "baseClass::readinputList(): inputList_ =  "<< *inputList_ << std::endl;

  ifstream is(inputList_->c_str());
  if(is.good())
    {
      cout << "baseClass::readInputList: Reading list: " << *inputList_ << " ......." << endl;
      while( is.getline(pName, 500, '\n') )
        {
          if (pName[0] == '#') continue;
          //if (pName[0] == ' ') continue; // do we want to skip lines that start with a space?
          if (pName[0] == '\n') continue;// simple protection against blank lines
          STDOUT("Adding file: " << pName);
          chain->Add(pName);
          NBeforeSkim = getGlobalInfoNstart(pName);
          NBeforeSkim_ = NBeforeSkim_ + NBeforeSkim;
          STDOUT("Initial number of events: NBeforeSkim, NBeforeSkim_ = "<<NBeforeSkim<<", "<<NBeforeSkim_);
        }

      tree_ = chain;
      cout << "baseClass::readInputList: Finished reading list: " << *inputList_ << endl;
    }
  else
    {
      cout << "baseClass::readInputList: ERROR opening inputList:" << *inputList_ << endl;
      exit (1);
    }
  is.close();

}

void baseClass::readCutFile()
{
  string s;
  STDOUT("Reading cutFile_ = "<< *cutFile_)

  ifstream is(cutFile_->c_str());
  if(is.good())
    {
      //      STDOUT("Reading file: " << *cutFile_ );
      int id=0;
      int optimize_count=0;
      while( getline(is,s) )
        {
          STDOUT("read line: " << s);
          if (s[0] == '#' || s.empty()) continue;
          vector<string> v = split(s);
          if ( v.size() == 0 ) continue;
          if (v[1]=="OPT") // add code for grabbing optimizer objects
            {
              if (optimizeName_cut_.size()>=6)
                {
                  STDOUT("WARNING:  Optimizer can only accept up to 6 variables.\nVariable "<<v[0]<<" is not being included.");
                  continue;
                }
              bool found=false;
              for (int i=0;i<optimizeName_cut_.size();++i)
                {
                  if (optimizeName_cut_[i].variableName==v[0])
                    {
                      STDOUT("ERROR:  variableName = "<<v[0]<<" is already being optimized in optimizedName_cut_.  Skipping.");
                      found=true;
                      break;
                    }
                }
              if (found) continue;

              int level_int = atoi(v[5].c_str());
              bool greaterthan=true;
              if (v[2]=="<") greaterthan=false;
              double minval=atof(v[3].c_str());
              double maxval=atof(v[4].c_str());
              Optimize opt(optimize_count,v[0],minval, maxval, greaterthan, level_int);
              optimizeName_cut_[optimize_count]=opt; // order cuts by cut #, rather than name, so that optimization histogram is consistently ordered
              ++optimize_count;
              continue;
            }

          map<string, cut>::iterator cc = cutName_cut_.find(v[0]);
          if( cc != cutName_cut_.end() )
            {
              STDOUT("ERROR: variableName = "<< v[0] << " exists already in cutName_cut_. Returning.");
              return;
            } 

          int level_int = atoi( v[5].c_str() );
          if(level_int == -1)
            {
              map<string, preCut>::iterator cc = preCutName_cut_.find(v[0]);
              if( cc != preCutName_cut_.end() )
                {
                  STDOUT("ERROR: variableName = "<< v[0] << " exists already in preCutName_cut_. Returning.");
                  return;
                } 
              preCutInfo_ << "### Preliminary cut values: " << s <<endl;
              preCut thisPreCut;
              thisPreCut.variableName =     v[0];
              thisPreCut.value1  = decodeCutValue( v[1] );
              thisPreCut.value2  = decodeCutValue( v[2] );
              thisPreCut.value3  = decodeCutValue( v[3] );
              thisPreCut.value4  = decodeCutValue( v[4] );
              preCutName_cut_[thisPreCut.variableName]=thisPreCut;
              continue;
            }
          cut thisCut;
          thisCut.variableName =     v[0];
          string m1=v[1];
          string M1=v[2];
          string m2=v[3];
          string M2=v[4];
          if( m1=="-" || M1=="-" ) 
            {
              STDOUT("ERROR: minValue1 and maxValue2 have to be provided. Returning."); 
              return; // FIXME implement exception
            } 
          if( (m2=="-" && M2!="-") || (m2!="-" && M2=="-") ) 
            {
              STDOUT("ERROR: if any of minValue2 and maxValue2 is -, then both have to be -. Returning");
              return; // FIXME implement exception
            }
          if( m2=="-") m2="+inf";
          if( M2=="-") M2="-inf";
          thisCut.minValue1  = decodeCutValue( m1 );
          thisCut.maxValue1  = decodeCutValue( M1 );
          thisCut.minValue2  = decodeCutValue( m2 );
          thisCut.maxValue2  = decodeCutValue( M2 );
          thisCut.level_int  = level_int;
          thisCut.level_str  =       v[5];
          thisCut.histoNBins = atoi( v[6].c_str() );
          thisCut.histoMin   = atof( v[7].c_str() );
          thisCut.histoMax   = atof( v[8].c_str() );
          // Not filled from file
          thisCut.id=++id;
          string s1;
          if(skimWasMade_)
            {
              s1 = "cutHisto_skim___________________" + thisCut.variableName;
            }
          else
            {
              s1 = "cutHisto_noCuts_________________" + thisCut.variableName;
            }
          string s2 = "cutHisto_allPreviousCuts________" + thisCut.variableName;
          string s3 = "cutHisto_allOthrSmAndLwrLvlCuts_" + thisCut.variableName;
          string s4 = "cutHisto_allOtherCuts___________" + thisCut.variableName;
          string s5 = "cutHisto_allCuts________________" + thisCut.variableName;
          thisCut.histo1 = TH1F (s1.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
          thisCut.histo2 = TH1F (s2.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
          thisCut.histo3 = TH1F (s3.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
          thisCut.histo4 = TH1F (s4.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
          thisCut.histo5 = TH1F (s5.c_str(),"", thisCut.histoNBins, thisCut.histoMin, thisCut.histoMax);
          thisCut.histo1.Sumw2();
          thisCut.histo2.Sumw2();
          thisCut.histo3.Sumw2();
          thisCut.histo4.Sumw2();
          thisCut.histo5.Sumw2();
          // Filled event by event
          thisCut.filled = false;
          thisCut.value = 0;
          thisCut.passed = false;
          thisCut.nEvtInput=0;
          thisCut.nEvtPassed=0;

          orderedCutNames_.push_back(thisCut.variableName);
          cutName_cut_[thisCut.variableName]=thisCut;

        }
      STDOUT( "baseClass::readCutFile: Finished reading cutFile: " << *cutFile_ );
    }
  else
    {
      STDOUT("ERROR opening cutFile:" << *cutFile_ );
      exit (1);
    }
  // make optimizer histogram
  if (optimizeName_cut_.size()>0)
    {
      h_optimizer_=new TH1F("optimizer","Optimization of cut variables",(int)pow(10,optimizeName_cut_.size()),0,
                            pow(10,optimizeName_cut_.size()));
    }

  // Create a histogram that will show events passing cuts
  int cutsize=orderedCutNames_.size()+1;
  if (skimWasMade_) ++cutsize;
  gDirectory->cd();
  eventcuts_=new TH1F("EventsPassingCuts","Events Passing Cuts",cutsize,0,cutsize);

  is.close();

}

void baseClass::resetCuts()
{
  for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++) 
    {
      cut * c = & (cc->second);
      c->filled = false;
      c->value = 0;
      c->passed = false;
    }
  combCutName_passed_.clear();
  return;
}

void baseClass::fillVariableWithValue(const string& s, const double& d)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      //STDOUT("variableName = "<< s << " not found in cutName_cut_.");
      return;
    } 
  else
    {
      cut * c = & (cc->second);
      c->filled = true;
      c->value = d;
    }
  fillOptimizerWithValue(s, d);
  return;
}

bool baseClass::variableIsFilled(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
    }
  cut * c = & (cc->second);
  return (c->filled);
}

double baseClass::getVariableValue(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_.");
    }
  cut * c = & (cc->second);
  if( !variableIsFilled(s) )
    {
      STDOUT("ERROR: requesting value of not filled variable "<<s);
    }
  return (c->value);
}

void baseClass::fillOptimizerWithValue(const string& s, const double& d)
{
  for (int i=0;i<optimizeName_cut_.size();++i)
    {
      if (optimizeName_cut_[i].variableName==s)
        {
          optimizeName_cut_[i].value=d;
          return;
        }
    }
  return;
}

void baseClass::evaluateCuts()
{
  //  resetCuts();
  combCutName_passed_.clear();
  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      if( ! ( c->filled && (c->minValue1 < c->value && c->value <= c->maxValue1 || c->minValue2 < c->value && c->value <= c->maxValue2 ) ) )
        {
          c->passed = false;
          combCutName_passed_[c->level_str.c_str()] = false;
          combCutName_passed_["all"] = false;
        }
      else
        {
          c->passed = true;
          map<string,bool>::iterator cp = combCutName_passed_.find( c->level_str.c_str() );
          combCutName_passed_[c->level_str.c_str()] = (cp==combCutName_passed_.end()?true:cp->second);
          map<string,bool>::iterator ap = combCutName_passed_.find( "all" );
          combCutName_passed_["all"] = (ap==combCutName_passed_.end()?true:ap->second);
        }
    }

  // reset optimization cut values
  //for (int i=0;i<optimizeName_cut_.size();++i)
  //  optimizeName_cut_[i].value=0;
  runOptimizer();

  if( !fillCutHistos() )
    {
      STDOUT("ERROR: fillCutHistos did not complete successfully.");
    }

  if( !updateCutEffic() )
    {
      STDOUT("ERROR: updateCutEffic did not complete successfully.");
    }

  return ;
}

void baseClass::runOptimizer()
{

  // don't run optimizer if no optimized cuts specified
  if (optimizeName_cut_.size()==0)
    return;

  // first, check that all cuts (except those to be optimized) have been passed

  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      bool ignorecut=false;
      for (unsigned int i=0; i < optimizeName_cut_.size();++i)
        {
          const string str = (const string)(*it);
          if (optimizeName_cut_[i].variableName.compare(str)==0)
            {
              ignorecut=true;
              break;
            }
        }
      if (ignorecut) continue;
      if (passedCut(*it) == false)
        return;
    }

  /*
  if (combCutName_passed_["all"] == false)
    return;
  */

  // loop over up to 6 cuts
  int counter=0;
  int thesize=optimizeName_cut_.size();
  int mysize=thesize;
  std::vector<bool> counterbins;
  for (int i=0;i<pow(10,thesize);++i) counterbins.push_back(true); // assume true
  
  // lowest-numbered cut appears first in cut ordering
  // That is, for cut:  ABCDEF
  //  A is the index of cut0, B is cut 1, etc.
  for (int cc=0;cc<thesize;++cc) // loop over all cuts, starting at cut 0
    {
      --mysize;
      for (int i=0;i<10;++i) // loop over 10 cuts for each
        {
          if (!optimizeName_cut_[cc].Compare(i)) // cut failed; set all values associated with cut to false
            {
              // loop over  all cut values starting with current cut
              for (unsigned int j=(int)(i*pow(10,mysize));j<int(pow(10,thesize));++j)
                {
                  // if relevant digit of the cut value matches the current (failed) cut, set this cut to false 
                  if ((j/int(pow(10,mysize)))%10==i)
                    counterbins[j]=false;
                  if (j>counterbins.size())
                    continue; // shouldn't ever happen
                }
            } // if (cut comparison failed)
        } // for (int i=0;i<10;++i)
    }
  // now fill histograms
  for (int i=0;i<counterbins.size();++i)
    {
      if (counterbins[i]==true)
        h_optimizer_->Fill(i,1);
      
    }

  return;
} //runOptimizer

bool baseClass::passedCut(const string& s)
{
  bool ret = false;
  //  map<string, bool>::iterator vp = cutName_passed_.find(s);
  //  if( vp != cutName_passed_.end() ) return ret = vp->second;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc != cutName_cut_.end() ) 
    {
      cut * c = & (cutName_cut_.find(s)->second);
      return (c->filled && c->passed);
    }
  map<string, bool>::iterator cp = combCutName_passed_.find(s);
  if( cp != combCutName_passed_.end() ) 
    {
      return ret = cp->second;
    }
  STDOUT("ERROR: did not find variableName = "<<s<<" neither in cutName_cut_ nor combCutName_passed_. Returning false.");
  return (ret=false);
}

bool baseClass::passedAllPreviousCuts(const string& s)
{
  //STDOUT("Examining variableName = "<<s);

  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {  
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
      return false;
    }

  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      if( c->variableName == s ) 
        {
          return true;
        }
      else
        {
          if( ! (c->filled && c->passed) ) return false;
        }
    }
  STDOUT("ERROR. It should never pass from here.");
}

bool baseClass::passedAllOtherCuts(const string& s)
{
  //STDOUT("Examining variableName = "<<s);
  bool ret = true;

  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {  
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
      return false;
    }

  for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++) 
    {
      cut * c = & (cc->second);
      if( c->variableName == s ) 
        {
          continue;
        }
      else
        {
          if( ! (c->filled && c->passed) ) return false;
        }
    }
  return ret;
}

bool baseClass::passedAllOtherSameAndLowerLevelCuts(const string& s)
{
  //STDOUT("Examining variableName = "<<s);
  bool ret = true;
  int cutLevel;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {  
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
      return false;
    }
  else
    {
      cutLevel = cc->second.level_int; 
    }

  for (map<string, cut>::iterator cc = cutName_cut_.begin(); cc != cutName_cut_.end(); cc++) 
    {
      cut * c = & (cc->second);
      if( c->level_int > cutLevel || c->variableName == s ) 
        {
          continue;
        }
      else 
        {
          if( ! (c->filled && c->passed) ) return false;
        }
    }
  return ret;
}

double baseClass::getPreCutValue1(const string& s)
{
  double ret;
  map<string, preCut>::iterator cc = preCutName_cut_.find(s);
  if( cc == preCutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
    }
  preCut * c = & (cc->second);
  if(c->value1 == -9999999) STDOUT("ERROR: value1 of preliminary cut "<<s<<" was not provided.");
  return (c->value1);
}

double baseClass::getPreCutValue2(const string& s)
{
  double ret;
  map<string, preCut>::iterator cc = preCutName_cut_.find(s);
  if( cc == preCutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
    }
  preCut * c = & (cc->second);
  if(c->value2 == -9999999) STDOUT("ERROR: value2 of preliminary cut "<<s<<" was not provided.");
  return (c->value2);
}

double baseClass::getPreCutValue3(const string& s)
{
  double ret;
  map<string, preCut>::iterator cc = preCutName_cut_.find(s);
  if( cc == preCutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
    }
  preCut * c = & (cc->second);
  if(c->value3 == -9999999) STDOUT("ERROR: value3 of preliminary cut "<<s<<" was not provided.");
  return (c->value3);
}

double baseClass::getPreCutValue4(const string& s)
{
  double ret;
  map<string, preCut>::iterator cc = preCutName_cut_.find(s);
  if( cc == preCutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in preCutName_cut_. Returning");
    }
  preCut * c = & (cc->second);
  if(c->value4 == -9999999) STDOUT("ERROR: value4 of preliminary cut "<<s<<" was not provided.");
  return (c->value4);
}


double baseClass::getCutMinValue1(const string& s)
{
  double ret;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
    }
  cut * c = & (cc->second);
  return (c->minValue1);
}

double baseClass::getCutMaxValue1(const string& s)
{
  double ret;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
    }
  cut * c = & (cc->second);
  return (c->maxValue1);
}

double baseClass::getCutMinValue2(const string& s)
{
  double ret;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
    }
  cut * c = & (cc->second);
  return (c->minValue2);
}

double baseClass::getCutMaxValue2(const string& s)
{
  double ret;
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() ) 
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning");
    }
  cut * c = & (cc->second);
  return (c->maxValue2);
}

const TH1F& baseClass::getHisto_noCuts_or_skim(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }

  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histo1);
}

const TH1F& baseClass::getHisto_allPreviousCuts(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }

  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histo2);
}

const TH1F& baseClass::getHisto_allOthrSmAndLwrLvlCuts(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }

  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histo3);
}

const TH1F& baseClass::getHisto_allOtherCuts(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }

  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histo4);
}

const TH1F& baseClass::getHisto_allCuts(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }

  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histo5);
}

int baseClass::getHistoNBins(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }
  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histoNBins);
}

double baseClass::getHistoMin(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }
  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histoMin);
}

double baseClass::getHistoMax(const string& s)
{
  map<string, cut>::iterator cc = cutName_cut_.find(s);
  if( cc == cutName_cut_.end() )
    {
      STDOUT("ERROR: did not find variableName = "<<s<<" in cutName_cut_. Returning false.");
    }
  cut * c = & (cutName_cut_.find(s)->second);
  return (c->histoMax);
}


bool baseClass::fillCutHistos()
{
  bool ret = true;
  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      if( c->filled )
        {
          c->histo1.Fill( c->value );
          if( passedAllPreviousCuts(c->variableName) )                c->histo2.Fill( c->value );
          if( passedAllOtherSameAndLowerLevelCuts(c->variableName) )  c->histo3.Fill( c->value );
          if( passedAllOtherCuts(c->variableName) )                   c->histo4.Fill( c->value );
          if( passedCut("all") )                                      c->histo5.Fill( c->value );
        }
    }
  return ret;
}

bool baseClass::writeCutHistos()
{ 
  bool ret = true;
  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      c->histo1.Write();
      c->histo2.Write();
      c->histo3.Write();
      c->histo4.Write();
      c->histo5.Write();
    }

  // Any failure mode to implement?
  return ret;
}

bool baseClass::updateCutEffic()
{
  bool ret = true;
  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      if( passedAllPreviousCuts(c->variableName) )  
        {
          c->nEvtInput++;
          if( passedCut(c->variableName) )              c->nEvtPassed++;
        }
    }
  return ret;
}


bool baseClass::writeCutEfficFile()
{

  bool ret = true;
  
  // Set bin labels for event counter histogram
  int bincounter=1;
  eventcuts_->GetXaxis()->SetBinLabel(bincounter,"NoCuts");
  ++bincounter;
  if (skimWasMade_)
    {
      eventcuts_->GetXaxis()->SetBinLabel(bincounter,"Skim");
      ++bincounter;
    }
  for (int i=0;i<orderedCutNames_.size();++i)
    {
      eventcuts_->GetXaxis()->SetBinLabel(bincounter,orderedCutNames_[i].c_str());
      ++bincounter;
    }

  bincounter=1;

  int nEntRoottuple = fChain->GetEntriesFast();
  int nEntTot = (skimWasMade_ ? NBeforeSkim_ : nEntRoottuple );
  string cutEfficFile = *cutEfficFile_ + ".dat";
  ofstream os(cutEfficFile.c_str());

  os << "################################## Preliminary Cut Values ###################################################################\n"
     << "########################### variableName                        value1          value2          value3          value4          level\n"
     << preCutInfo_.str();

  int cutIdPed=0;
  os.precision(4); 
  os << "################################## Cuts #####################################################################################\n" 
     <<"#id             variableName           min1           max1           min2           max2          level              N          Npass         EffRel      errEffRel         EffAbs      errEffAbs"<<endl
     << fixed
     << setw(3) << cutIdPed 
     << setw(25) << "nocut" 
     << setprecision(4) 
     << setw(15) << "-"
     << setw(15) << "-"
     << setw(15) << "-"
     << setw(15) << "-"
     << setw(15) << "-"
     << setw(15) << nEntTot
     << setw(15) << nEntTot
     << setprecision(11) 
     << setw(15) << "1.00000000000"
     << setw(15) << "0.00000000000"
     << setw(15) << "1.00000000000"
     << setw(15) << "0.00000000000"
     << endl;

  double effRel;
  double effRelErr;
  double effAbs;
  double effAbsErr;

  eventcuts_->SetBinContent(bincounter,nEntTot);
  if (optimizeName_cut_.size())
    h_optimizer_->SetBinContent(0, nEntTot);

  if(skimWasMade_)
    {
      ++bincounter;
      eventcuts_->SetBinContent(bincounter, nEntRoottuple);
      effRel = (double) nEntRoottuple / (double) NBeforeSkim_;
      effRelErr = sqrt( (double) effRel * (1.0 - (double) effRel) / (double) NBeforeSkim_ );
      effAbs = effRel;
      effAbsErr = effRelErr;
      os << fixed
         << setw(3) << ++cutIdPed
         << setw(25) << "skim" 
         << setprecision(4) 
         << setw(15) << "-"
         << setw(15) << "-"
         << setw(15) << "-"
         << setw(15) << "-"
         << setw(15) << "-"
         << setw(15) << NBeforeSkim_
         << setw(15) << nEntRoottuple
         << setprecision(11) 
         << setw(15) << effRel
         << setw(15) << effRelErr
         << setw(15) << effAbs
         << setw(15) << effAbsErr
         << endl;
    }
  for (vector<string>::iterator it = orderedCutNames_.begin(); 
       it != orderedCutNames_.end(); it++) 
    {
      cut * c = & (cutName_cut_.find(*it)->second);
      ++bincounter;
      eventcuts_->SetBinContent(bincounter, c->nEvtPassed);
      effRel = (double) c->nEvtPassed / (double) c->nEvtInput;
      effRelErr = sqrt( (double) effRel * (1.0 - (double) effRel) / (double) c->nEvtInput );
      effAbs = (double) c->nEvtPassed / (double) nEntTot;
      effAbsErr = sqrt( (double) effAbs * (1.0 - (double) effAbs) / (double) nEntTot );

      std::stringstream ssm1, ssM1, ssm2,ssM2;
      ssm1 << fixed << setprecision(4) << c->minValue1;
      ssM1 << fixed << setprecision(4) << c->maxValue1;
      if(c->minValue2 == -9999999) 
        {
          ssm2 << "-inf";
        }
      else
        {
          ssm2 << fixed << setprecision(4) << c->minValue2;
        }
      if(c->maxValue2 ==  9999999) 
        {
          ssM2 << "+inf";
        }
      else
        {
          ssM2 << fixed << setprecision(4) << c->maxValue2;
        }
      os << setw(3) << cutIdPed+c->id 
         << setw(25) << c->variableName 
         << setprecision(4)
         << fixed 
         << setw(15) << ( ( c->minValue1 == -9999999.0 ) ? "-inf" : ssm1.str() )
         << setw(15) << ( ( c->maxValue1 ==  9999999.0 ) ? "+inf" : ssM1.str() )
         << setw(15) << ( ( c->minValue2 > c->maxValue2 ) ? "-" : ssm2.str() )
         << setw(15) << ( ( c->minValue2 > c->maxValue2 ) ? "-" : ssM2.str() )
         << setw(15) << c->level_int
         << setw(15) << c->nEvtInput
         << setw(15) << c->nEvtPassed
         << setprecision(11) 
         << setw(15) << effRel
         << setw(15) << effRelErr
         << setw(15) << effAbs
         << setw(15) << effAbsErr
         << endl;
    }

  // Write optimization histograms
  if (optimizeName_cut_.size())
    {
      gDirectory->mkdir("Optimizer");
      gDirectory->cd("Optimizer");
      h_optimizer_->Write();
      for (int i=0;i<optimizeName_cut_.size();++i)
        {
          stringstream x;
          x<<"Cut"<<i<<"_"<<optimizeName_cut_[i].variableName;
          if (optimizeName_cut_[i].testgreater==true)
            x<<"_gt_";
          else
            x<<"_lt_";
          x<<optimizeName_cut_[i].minvalue<<"_to_"<<optimizeName_cut_[i].maxvalue;
          TObjString test(x.str().c_str());
          test.Write();
        }
      gDirectory->cd("..");
    }

  eventcuts_->Write(); // write here, since WriteCutHistos is called before WriteCutEfficFile
  // Any failure mode to implement?
  return ret;
} // writeCutEffFile


bool baseClass::sortCuts(const cut& X, const cut& Y)
{
  return X.id < Y.id;
}


vector<string> baseClass::split(const string& s)
{
  vector<string> ret;
  string::size_type i =0;
  while (i != s.size()){
    while (i != s.size() && isspace(s[i]))
      ++i;
    string::size_type j = i;
    while (j != s.size() && !isspace(s[j]))
      ++j;
    if (i != j){
      ret.push_back(s.substr(i, j -i));
      i = j;
    }
  }
  return ret;
}

double baseClass::decodeCutValue(const string& s)
{
  //STDOUT("s = "<<s);
  double ret;
  if( s == "inf" || s == "+inf" )
    {
       ret = 9999999;
    }
  else if ( s == "-inf" || s == "-" )
    {
       ret = -9999999;
    }
  else
    {
       ret = atof( s.c_str() );
    }
  return ret;
}

int baseClass::getGlobalInfoNstart(char *pName)
{
  int NBeforeSkim = 0;
  STDOUT(pName<<"  "<< NBeforeSkim)
  TFile *f = new TFile(pName);
  string s1 = "LJFilter/EventCount/EventCounter";
  string s2 = "LJFilterPAT/EventCount/EventCounter";
  TH1I* hCount1 = (TH1I*)f->Get(s1.c_str());
  TH1I* hCount2 = (TH1I*)f->Get(s2.c_str());
  if( !hCount1 && !hCount2 ) 
    {
      STDOUT("Skim filter histogram(s) not found. Will assume skim was not made for ALL files.");
      skimWasMade_ = false;
      return NBeforeSkim;
    }
    
  if (hCount1) NBeforeSkim = (int)hCount1->GetBinContent(1);
  else NBeforeSkim = (int)hCount2->GetBinContent(1);
  
//   STDOUT(pName<<"  "<< NBeforeSkim) 
  f->Close();

  return NBeforeSkim;
}


Revision:	1.2
Committed:	Tue Jul 20 15:32:43 2010 UTC (14 years, 9 months ago) by jueugste
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+358 -25 lines
Log Message:	many changes...