rootEWKanalyzer/src/analysisClass_EWK.C


#define analysisClass_cxx
#include "analysisClass.h"

#include <TH2.h>
#include <TF2.h>
#include <TH1F.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TVector2.h>
#include <TVector3.h>
#include <TTree.h>

// avg. inst. lumi.
#include "readLumiInfo/readLumiInfo.h"

analysisClass::analysisClass(string * inputList, string * cutFile, string * treeName, string * outputFileName, string * cutEfficFile)
  :baseClass(inputList, cutFile, treeName, outputFileName, cutEfficFile)
{
  std::cout << "analysisClass::analysisClass(): begins " << std::endl;

  std::cout << "analysisClass::analysisClass(): ends " << std::endl;
}

analysisClass::~analysisClass()
{
  std::cout << "analysisClass::~analysisClass(): begins " << std::endl;

  std::cout << "analysisClass::~analysisClass(): ends " << std::endl;
}


void analysisClass::Loop()
{
  std::cout << "analysisClass::Loop() begins" <<std::endl;   
  
  if (fChain == 0) return;
  

  double PtCorrection = 1.; //0.99;


  int counterSize=15;
  int* counter = new int[counterSize];
  memset(counter,0,counterSize*sizeof(int));

  int El_counterSize=15;
  int* El_count = new int[El_counterSize];
  memset(El_count,0,El_counterSize*sizeof(int));
  
  FILE *out=fopen("out.txt","w");
  FILE *out_Z=fopen("out_z.txt","w");
  fprintf(out_Z,"Run \tEvent \tLS \tDiLeptonMass \tPt_1 \tPt2 \tEta_1 \tEta_2 \tDiLeptonPt \tsumJetPt \n");

  // open the *.root file containing the lumi info
//   readLumiInfo a("readLumiInfo/lumiByLSRun2010.root");


  // --------------------------------------------------------------------------------------------------------------------
  // Analysis input parameters
  // --------------------------------------------------------------------------------------------------------------------
  int workingPoint = 80;   // -1 takes hard coded values (which are the old robust tight values...)
  double sampleSize=-1;
 
  bool isdata;
  bool useElectronTiggers = true;


  bool ethIDandIso = false;

  bool writeIdPlots=true;
  bool writeN1Plots=true;
  bool writeWenuPlots=true;
  bool applyDeltaCorrection=false;
  // --------------------------------------------------------------------------------------------------------------------

  double Pi=2*acos(0);  // pi=3.14159...

  int MAXE=6;

  vector<double> WorkingPointCuts;
  vector<double> EthCuts;


  // --------------------------------------------------------------------------------------------------------------------
  // Index vectors
  // --------------------------------------------------------------------------------------------------------------------
  vector<int> i_preselected;
  vector<int> i_isolated;
  vector<int> i_identified;
  vector<int> i_goodElectrons;
  vector<int> i_goodBarrelElectrons;
  vector<int> i_goodEndcapElectrons;

  vector<int> i_jets;
  // --------------------------------------------------------------------------------------------------------------------


  // --------------------------------------------------------------------------------------------------------------------
  // tree variables
  // --------------------------------------------------------------------------------------------------------------------
//   TTree *subTree;
  TTree *subTree = new TTree("tree", "tree");

  int tRun;
  int tEvent;
  int tLumiSection;
  // electrons
  int ElMax = 20;
  int tElN;
  double tElPt[ElMax];
  double tElPx[ElMax];
  double tElPy[ElMax];
  double tElPz[ElMax];
  double tElE[ElMax];
  double tElEta[ElMax];
  double tElPhi[ElMax];
  int tElCharge[ElMax];
  double tElDz[ElMax];
  // MET
  double tCaloMet;
  double tTcMet;
  double tPfMet;
  // mass
  double tElPfMt[ElMax];
  double tElCaloMt[ElMax];
  double tElTcMt[ElMax];
  // HF stuff
  double tSumEHF_minus;
  double tSumEHF_plus;
  double tSumEHF_Aminus;
  double tSumEHF_Aplus;
  double tSumEHF_Bminus;
  double tSumEHF_Bplus;
  double tSumETHF_minus;
  double tSumETHF_plus;
  int tNtowers_plus;
  int tNtowers_minus;
  // calo towers
  int TowMax = 1000;
  int tCaloTowersN;
  double tCaloTowersE[TowMax];
  double tCaloTowersEt[TowMax];
  double tCaloTowersEta[TowMax];
  double tCaloTowersEmE[TowMax];
  double tCaloTowersHadE[TowMax];
  int tCaloTowersHasEB[TowMax];
  int tCaloTowersHasEE[TowMax];
  int tCaloTowersHasHB[TowMax];
  int tCaloTowersHasHE[TowMax];
  int tCaloTowersHasHF[TowMax];
  int tCaloTowersId[TowMax];

  // jets
  int JetMax = 20;
  int tJetN;
  double tJetPt[JetMax];
  double tJetEta[JetMax];
  double tJetPhi[JetMax];

  // diffractive variables
  double tEtaMax;
  double tEtaMin;
  double tEminusPz;
  double tEplusPz;
  double tSumE;
  double tSumPz;

  // PFlow candidates
  int PfMax = 1000;
  int tPfN;
  double tPfE[PfMax];
  double tPfEt[PfMax];
  double tPfEta[PfMax];
  double tPfPt[PfMax];


  // vertices
  int VertexMax = 25;
  int tVertexN;
//   int tVertexNdof[VertexMax];
  double tVertexNdof[VertexMax];
  double tVertexChi2[VertexMax];
  double tVertexX[VertexMax];
  double tVertexY[VertexMax];
  double tVertexZ[VertexMax];
  int tVertexTrackIndex[25][400];
//   float tVertexTrackWeight[25][400];

  // pile-up
  int tPUnumInteractions;
  double tPUzPositions[15];

  // beam spot
  double tBeamSpotX;
  double tBeamSpotY;
  double tBeamSpotZ;

  // Tracks
  int TrackMax = 400;
  int tTracksN;
  double tTrackPt[TrackMax];
  double tTrackPx[TrackMax];
  double tTrackPy[TrackMax];
  double tTrackPz[TrackMax];
  double tTrackDzBs[TrackMax];
  double tTrackEta[TrackMax];
  double tTrackPhi[TrackMax];
//   double tTrackdsz[TrackMax];
//   double tTrackdxy[TrackMax];


//   tRun           = -999;
//   tEvent         = -999;
//   tLumiSection   = -999;
  
//   double *tElPt = new double[ElMax];
//   memset(tElPt, 0, ElMax*sizeof(double));

//     int MAXWcounter=16;
//     int *Wcounter = new int[MAXWcounter];
//     memset(Wcounter,0,MAXWcounter*sizeof(int));


  // Run infos
  subTree->Branch("Run"            ,&tRun,         "Run/I");
  subTree->Branch("Event"          ,&tEvent,       "Event/I");
  subTree->Branch("LumiSection"    ,&tLumiSection, "LumiSection/I");
  // electrons
  subTree->Branch("ElN",            &tElN,         "ElN/I");
  subTree->Branch("ElPt",           &tElPt,        "ElPt[ElN]/D");
  subTree->Branch("ElPx",           &tElPx,        "ElPx[ElN]/D");
  subTree->Branch("ElPy",           &tElPy,        "ElPy[ElN]/D");
  subTree->Branch("ElPz",           &tElPz,        "ElPz[ElN]/D");
  subTree->Branch("ElE",            &tElE,         "ElE[ElN]/D");
  subTree->Branch("ElEta",          &tElEta,       "ElEta[ElN]/D");
  subTree->Branch("ElPhi",          &tElPhi,       "ElPhi[ElN]/D");
  subTree->Branch("ElCharge",       &tElCharge,    "ElCharge[ElN]/I");
  subTree->Branch("ElDz",           &tElDz,        "ElDz[ElN]/D");

  subTree->Branch("PfMet",          &tPfMet,       "PfMet/D");
  subTree->Branch("CaloMet",        &tCaloMet,     "CaloMet/D");
  subTree->Branch("TcMet",          &tTcMet,       "TcMet/D");

  subTree->Branch("ElPfMt",         &tElPfMt,      "ElPfMt[ElN]/D");
  subTree->Branch("ElCaloMt",       &tElCaloMt,    "ElCaloMt[ElN]/D");
  subTree->Branch("ElTcMt",         &tElTcMt,      "ElTcMt[ElN]/D");

  subTree->Branch("SumEHF_minus",   &tSumEHF_minus,"SumEHF_minus/D");
  subTree->Branch("SumEHF_plus",    &tSumEHF_plus, "SumEHF_plus/D");
  subTree->Branch("SumEHF_Aminus",  &tSumEHF_Aminus,"SumEHF_Aminus/D");
  subTree->Branch("SumEHF_Aplus",   &tSumEHF_Aplus, "SumEHF_Aplus/D");
  subTree->Branch("SumEHF_Bminus",  &tSumEHF_Bminus,"SumEHF_Bminus/D");
  subTree->Branch("SumEHF_Bplus",   &tSumEHF_Bplus, "SumEHF_Bplus/D");
  subTree->Branch("SumETHF_minus",  &tSumETHF_minus,"SumETHF_minus/D");
  subTree->Branch("SumETHF_plus",   &tSumETHF_plus,"SumETHF_plus/D");
  subTree->Branch("Ntowers_plus",   &tNtowers_plus,"Ntowers_plus/I");
  subTree->Branch("Ntowers_minus",  &tNtowers_minus,"Ntowers_minus/I");

  subTree->Branch("CaloTowersN",     &tCaloTowersN,     "CaloTowersN/I");
  subTree->Branch("CaloTowersE",     &tCaloTowersE,     "CaloTowersE[CaloTowersN]/D");
  subTree->Branch("CaloTowersEt",    &tCaloTowersEt,    "CaloTowersEt[CaloTowersN]/D");
  subTree->Branch("CaloTowersEta",   &tCaloTowersEta,   "CaloTowersEta[CaloTowersN]/D");
  subTree->Branch("CaloTowersEmE",   &tCaloTowersEmE,   "CaloTowersEmE[CaloTowersN]/D");
  subTree->Branch("CaloTowersHadE",  &tCaloTowersHadE,  "CaloTowersHadE[CaloTowersN]/D");
  subTree->Branch("CaloTowersHasEB", &tCaloTowersHasEB, "CaloTowersHasEB[CaloTowersN]/I");
  subTree->Branch("CaloTowersHasEE", &tCaloTowersHasEE, "CaloTowersHasEE[CaloTowersN]/I");
  subTree->Branch("CaloTowersHasHB", &tCaloTowersHasHB, "CaloTowersHasHB[CaloTowersN]/I");
  subTree->Branch("CaloTowersHasHE", &tCaloTowersHasHE, "CaloTowersHasHE[CaloTowersN]/I");
  subTree->Branch("CaloTowersHasHF", &tCaloTowersHasHF, "CaloTowersHasHF[CaloTowersN]/I");
  subTree->Branch("CaloTowersId",    &tCaloTowersId,    "CaloTowersId[CaloTowersN]/I");


  subTree->Branch("JetN",           &tJetN,        "JetN/I");
  subTree->Branch("JetPt",          &tJetPt,       "JetPt[JetN]/D");
  subTree->Branch("JetEta",         &tJetEta,      "JetEta[JetN]/D");
  subTree->Branch("JetPhi",         &tJetPhi,      "JetPhi[JetN]/D");

  subTree->Branch("EtaMax",         &tEtaMax,      "EtaMax/D");
  subTree->Branch("EtaMin",         &tEtaMin,      "EtaMin/D");

  subTree->Branch("EminusPz",       &tEminusPz,    "EminusPz/D");
  subTree->Branch("EplusPz",        &tEplusPz,     "EplusPz/D");
  subTree->Branch("SumE",           &tSumE,        "SumE/D");
  subTree->Branch("SumPz",          &tSumPz,       "SumPz/D");

  // pflow
  subTree->Branch("PfN",            &tPfN,         "PfN/I");
  subTree->Branch("PfE",            &tPfE,         "PfE[PfN]/D");
  subTree->Branch("PfEt",           &tPfEt,        "PfEt[PfN]/D");
  subTree->Branch("PfEta",          &tPfEta,       "PfEta[PfN]/D");
  subTree->Branch("PfPt",           &tPfPt,        "PfPt[PfN]/D");
  

  subTree->Branch("VertexN",        &tVertexN,     "VertexN/I");
//   subTree->Branch("VertexNdof",     &tVertexNdof,  "VertexNdof[VertexN]/I");
  subTree->Branch("VertexNdof",     &tVertexNdof,  "VertexNdof[VertexN]/D");
  subTree->Branch("VertexChi2",     &tVertexChi2,  "VertexChi2[VertexN]/I");
  subTree->Branch("VertexX",        &tVertexX,     "VertexX[VertexN]/D");
  subTree->Branch("VertexY",        &tVertexY,     "VertexY[VertexN]/D");
  subTree->Branch("VertexZ",        &tVertexZ,     "VertexZ[VertexN]/D");
  subTree->Branch("VertexTrackIndex",     &tVertexTrackIndex    ,  "VertexTrackIndex[25][400]/I");
//   subTree->Branch("VertexTrackWeight",    &fTVertexTrackWeight   ,  "VertexTrackWeight[25][400]/F");

  // pile-up
  subTree->Branch("PUnumInteractions",  &tPUnumInteractions,  "PUnumInteractions/I");
  subTree->Branch("PUzPositions",  &tPUzPositions,  "PUzPositions[PUnumInteractions]/D");


  subTree->Branch("BeamSpotX",      &tBeamSpotX,   "BeamSpotX/D");
  subTree->Branch("BeamSpotY",      &tBeamSpotY,   "BeamSpotY/D");
  subTree->Branch("BeamSpotZ",      &tBeamSpotZ,   "BeamSpotZ/D");


  subTree->Branch("TracksN",        &tTracksN,     "TracksN/I");
  subTree->Branch("TrackPt",        &tTrackPt,     "TrackPt[TracksN]/D");
  subTree->Branch("TrackPx",        &tTrackPx,     "TrackPx[TracksN]/D");
  subTree->Branch("TrackPy",        &tTrackPy,     "TrackPy[TracksN]/D");
  subTree->Branch("TrackPz",        &tTrackPz,     "TrackPz[TracksN]/D");
  subTree->Branch("TrackDzBs",        &tTrackDzBs,     "TrackDzBs[TracksN]/D");
  subTree->Branch("TrackEta",       &tTrackEta,    "TrackEta[TracksN]/D");
  subTree->Branch("TrackPhi",       &tTrackPhi,    "TrackPhi[TracksN]/D");
//   subTree->Branch("Trackdsz",       &tTrackdsz,    "Trackdsz[TracksN]/D");
//   subTree->Branch("Trackdxy",       &tTrackdxy,    "Trackdxy[TracksN]/D");
//   subTree->Branch("",        &,     "[]/D");


//   subTree->Branch("");
//   subTree->Branch("");
//   subTree->Branch("");
//   subTree->Branch("");
//   subTree->Branch("");
//   subTree->Branch("");
//   subTree->Branch("");

//     double tElN;
//   double tElPt;
//   double tElPx;
//   double tElPy;
//   double tElPz;
//   double tElEta;
//   double tElPhi;


  // --------------------------------------------------------------------------------------------------------------------


  // --------------------------------------------------------------------------------------------------------------------
  // Electron Cuts
  // --------------------------------------------------------------------------------------------------------------------
  
  // ETH ID
  double ETH_oneOverEminusOneOverP_B=0.005;
  double ETH_deltaEtaMax_B=0.004;
  double ETH_deltaPhiMax_B=0.02;
  double ETH_iso_B=0.1;
  int ETH_missingHits_B=1;
  double ETH_oneOverEminusOneOverP_E=0.007;
  double ETH_deltaEtaMax_E=0.006;
  double ETH_deltaPhiMax_E=0.02;
  double ETH_iso_E=0.1;
  int ETH_missingHits_E=1;

  
  int debug=0;
  int debug2=0;
  int debug3=0;
  int debug4=0;


  vector<int> diff_run;
  vector<int> diff_event;
  vector<int> diff_ls;
  vector<double> diff_sumehf_plus;
  vector<double> diff_sumehf_minus;
  vector<double> diff_sumethf_plus;
  vector<double> diff_sumethf_minus;
  vector<double> diff_ntowershf_plus;
  vector<double> diff_ntowershf_minus;

  vector<int> sel_WOjet_run;
  vector<int> sel_WOjet_event;
  vector<int> sel_WOjet_ls;
  vector<int> sel_ZOjet_run;
  vector<int> sel_ZOjet_event;
  vector<int> sel_ZOjet_ls;
  
  int previousRun = 0;
  bool firstEvent=true;
  // fill the eth cuts vector
  EthCuts.push_back(ETH_deltaPhiMax_B);
  EthCuts.push_back(ETH_deltaEtaMax_B);
  EthCuts.push_back(ETH_oneOverEminusOneOverP_B);
  EthCuts.push_back(ETH_iso_B);
  EthCuts.push_back(ETH_missingHits_B);
  EthCuts.push_back(ETH_deltaPhiMax_E);
  EthCuts.push_back(ETH_deltaEtaMax_E);
  EthCuts.push_back(ETH_oneOverEminusOneOverP_E);
  EthCuts.push_back(ETH_iso_E);
  EthCuts.push_back(ETH_missingHits_E);


  // --------------------------------------------------------------------------------------------------------------------
  // read the cut values from include/workingPoints.txt
  // --------------------------------------------------------------------------------------------------------------------

  // check if it is a possible working point
  if(workingPoint>0 && (
     workingPoint!=60 &&
     workingPoint!=70 &&
     workingPoint!=80 &&
     workingPoint!=85 &&
     workingPoint!=90 &&
     workingPoint!=95 )
     ) {
    cerr<<"ERROR: "<<workingPoint<<" is NOT a valid working point!!!"<<endl;
    exit(0);
  }
  // read the file
  if(workingPoint>0) {
    const char *file="include/workingPoints.txt";
    cout<<"Opening file: "<<file<<"..."<<endl;
    ifstream toread (file);
    if(!toread.is_open()) {
      cerr<<"File "<<file<<" not found!!"<<endl;
      //goto END;
    }
    else {
      cout<<"File opened..."<<endl;
    }
    
    // build the working point variable
    char wpB [20];
    char wpE [20];
    sprintf(wpB,"B%d",workingPoint);
    sprintf(wpE,"E%d",workingPoint);
    
    
    // read every line in the file
    if(toread.is_open()) {
      string line;
      float value;
      char name[100];
      //cout<<"Reading working point input file"<<endl;
      while(!toread.eof()) {
        getline (toread, line);
        
        if(line.substr(0,1)=="#") continue;
        //cout<<line<<endl;
        
        char wp[20];
        float v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10;
        sscanf (line.c_str(), "%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f", wp, &v0, &v1, &v2, &v3, &v4, &v5, &v6, &v7, &v8, &v9, &v10);
        // barrel values
        if(!strcmp(wp, wpB)) {
          WorkingPointCuts.push_back(v0); // 0 see
          WorkingPointCuts.push_back(v1); // 1 dphi
          WorkingPointCuts.push_back(v2); // 2 deta
          WorkingPointCuts.push_back(v3); // 3 hoe
          WorkingPointCuts.push_back(v4); // 4 tkIso
          WorkingPointCuts.push_back(v5); // 5 ecalIso
          WorkingPointCuts.push_back(v6); // 6 hcalIso
          WorkingPointCuts.push_back(v7); // 7 combIso
        }
        // endcap values
        if(!strcmp(wp, wpE)) {
          WorkingPointCuts.push_back(v0); // 8 see
          WorkingPointCuts.push_back(v1); // 9 dphi
          WorkingPointCuts.push_back(v2); // 10 deta
          WorkingPointCuts.push_back(v3); // 11 hoe
          WorkingPointCuts.push_back(v4); // 12 tkIso
          WorkingPointCuts.push_back(v5); // 13 ecalIso
          WorkingPointCuts.push_back(v6); // 14 hcalIso
          WorkingPointCuts.push_back(v7); // 15 combIso
          WorkingPointCuts.push_back(v8); // 16 hits
          WorkingPointCuts.push_back(v9); // 17 dist
        }
        
      } // read line
    }  // read file
  } // get working points
  
  cout<<WorkingPointCuts[4]<<endl;
  cout<<WorkingPointCuts[5]<<endl;
  cout<<WorkingPointCuts[6]<<endl;
  cout<<WorkingPointCuts[7]<<endl;
  cout<<WorkingPointCuts[8]<<endl;
  cout<<WorkingPointCuts[9]<<endl;


  cout<<"Running the analysis with working point "<<workingPoint<<endl;
  cout<<"Value for barrel / endcap"<<endl;
  char *cutName[8] = {"sigmaIetaIeta", "delta phi at vtx", "delta eta at vtx", "H/E", "trkIso", "ecalIso", "hcalIso", "# inner hits"};
  for(int i=0; i< 8;i++) {
    cout<<cutName[i]<<" < "<<WorkingPointCuts[i]<<" / "<<WorkingPointCuts[8+i]<<endl;
  }

  // --------------------------------------------------------------------------------------------------------------------
  

  // --------------------------------------------------------------------------------------------------------------------
  // number of events
  // this histogram is filled for normalization purposes i.e. to keep track of the MC sample sizes
  // --------------------------------------------------------------------------------------------------------------------
  TH1I *h_numberOfEvents;
  h_numberOfEvents=new TH1I("numberOfEvents", "numberOfEvents", 1, -0.01, 1.01);
  TH1I *h_HLTInfo;
  h_HLTInfo = new TH1I("HLTInfo", "HLTInfo", 6, -.5, 5.5);
  h_HLTInfo->GetXaxis()->SetBinLabel(1, "No trigger found");
  h_HLTInfo->GetXaxis()->SetBinLabel(2, "Ele15_SW");
  h_HLTInfo->GetXaxis()->SetBinLabel(2, "Ele15_LW");
  h_HLTInfo->GetXaxis()->SetBinLabel(2, "Ele15_SW_CaloEleId");
  h_HLTInfo->GetXaxis()->SetBinLabel(2, "Ele17_SW_TightEleId");
  h_HLTInfo->GetXaxis()->SetBinLabel(2, "Ele17_SW_CaloEleId");
  
   // --------------------------------------------------------------------------------------------------------------------


   // --------------------------------------------------------------------------------------------------------------------
   // Electron plots
   // --------------------------------------------------------------------------------------------------------------------
   char name [50];
   int nbins=12;
   double start=0.;
   double end=120.;
   TH1D *hs_El_ElPt[MAXE];
   TH1D *hs_El_ElEta[MAXE];
   TH1D *hs_El_ElPhi[MAXE];
   TH1D *hs_El_ElDeltaPhi[MAXE];
   TH1D *hs_El_ElDeltaEta[MAXE];
   TH1D *hs_El_ElOneOverEminusOneOverP[MAXE];
   TH1D *hs_El_ElMet[MAXE];
   TH1D *hs_El_ElTransMass[MAXE];
   TH1D *hs_El_ElCharge[MAXE];
   TH1D *hs_El_Njets[MAXE];
   TH1D *hs_El_ElEthIso[MAXE];
   TH1D *hs_El_ElTrkIso[MAXE];
   TH1D *hs_El_ElEcalIso[MAXE];
   TH1D *hs_El_ElHcalIso[MAXE];
   TH1D *hs_El_ElCombIso[MAXE];
   TH1D *hs_El_ElHoverE[MAXE];
   TH1D *hs_El_ElSigmaIetaIeta[MAXE];
   TH1D *hs_El_ElNumberOfMissingInnerHits[MAXE];
    for(int i=0; i<MAXE; i++) {
      sprintf(name,"s_El_ElPt_%d",i);
      hs_El_ElPt[i]=new TH1D(name, name, 60, 0, 120);
      hs_El_ElPt[i]->SetXTitle("p_{T} [GeV]");
      hs_El_ElPt[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_El_ElEta_%d",i);
      hs_El_ElEta[i]=new TH1D(name, name, 25, -2.5, 2.5);
      hs_El_ElEta[i]->SetXTitle("#eta");
      hs_El_ElEta[i]->SetYTitle("events / 0.2");
      sprintf(name,"s_El_ElPhi_%d",i);
      hs_El_ElPhi[i]=new TH1D(name, name, 30, -Pi,Pi);
      hs_El_ElPhi[i]->SetXTitle("#phi");
      hs_El_ElPhi[i]->SetYTitle("events / 0.21");
      sprintf(name,"s_El_ElDeltaPhi_%d",i);
      hs_El_ElDeltaPhi[i]=new TH1D(name, name, 30, -0.05001, 0.05001);
      hs_El_ElDeltaPhi[i]->SetXTitle("#Delta #phi");
      hs_El_ElDeltaPhi[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElDeltaEta_%d",i);
      hs_El_ElDeltaEta[i]=new TH1D(name, name, 30, -0.05001, 0.05001);
      hs_El_ElDeltaEta[i]->SetXTitle("#Delta #eta");
      hs_El_ElDeltaEta[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElMet_%d",i);
      hs_El_ElMet[i]=new TH1D(name, name, 60, 0, 120);
      hs_El_ElMet[i]->SetXTitle("#slash{E}_{T} [GeV]");
      hs_El_ElMet[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_El_ElTransMass_%d",i);
      hs_El_ElTransMass[i]=new TH1D(name, name, 60, 0, 120);
      hs_El_ElTransMass[i]->SetXTitle("m_{T} [GeV]");
      hs_El_ElTransMass[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_El_ElCharge_%d",i);
      hs_El_ElCharge[i]=new TH1D(name, name, 2, -2, 2);
      hs_El_ElCharge[i]->SetXTitle("Charge");
      hs_El_ElCharge[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_Njets_%d",i);
      hs_El_Njets[i]=new TH1D(name, name, 10, 0, 10);
      hs_El_Njets[i]->SetXTitle("number of jets");
      hs_El_Njets[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElOneOverEminusOneOverP_%d",i);
      hs_El_ElOneOverEminusOneOverP[i]=new TH1D(name, name, 30, -0.01, 0.5);
      hs_El_ElOneOverEminusOneOverP[i]->SetXTitle("|1/E-1/p|");
      hs_El_ElOneOverEminusOneOverP[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElEthIso_%d",i);
//       hs_El_ElEthIso[i]=new TH1D(name, name, 50, -0.01, 5.);
      hs_El_ElEthIso[i]=new TH1D(name, name, 20, -0.01, 1.1);
      hs_El_ElEthIso[i]->SetXTitle("ETH Iso");
      hs_El_ElEthIso[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElTrkIso_%d",i);
//       hs_El_ElTrkIso[i]=new TH1D(name, name, 100, -0.01, 50.);
      hs_El_ElTrkIso[i]=new TH1D(name, name, 20, -0.01, 1.1);
      hs_El_ElTrkIso[i]->SetXTitle("Trk Iso");
      hs_El_ElTrkIso[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElEcalIso_%d",i);
      hs_El_ElEcalIso[i]=new TH1D(name, name, 20, -0.01, 1.1);
//       hs_El_ElEcalIso[i]=new TH1D(name, name, 40, -0.01, 20.);
      hs_El_ElEcalIso[i]->SetXTitle("ECAL Iso");
      hs_El_ElEcalIso[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElHcalIso_%d",i);
      hs_El_ElHcalIso[i]=new TH1D(name, name, 20, -0.01, 10.);
      hs_El_ElHcalIso[i]->SetXTitle("HCAL Iso");
      hs_El_ElHcalIso[i]->SetYTitle("events / bin");

      sprintf(name,"s_El_ElCombIso_%d",i);
      hs_El_ElCombIso[i]=new TH1D(name, name, 20, -0.01, 3.1);
      hs_El_ElCombIso[i]->SetXTitle("Combined Iso");
      hs_El_ElCombIso[i]->SetYTitle("events / bin");

      sprintf(name,"s_El_ElHoverE_%d",i);
      hs_El_ElHoverE[i]=new TH1D(name, name, 30, -0.01, 0.5);
      hs_El_ElHoverE[i]->SetXTitle("H/E");
      hs_El_ElHoverE[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElSigmaIetaIeta_%d",i);
      hs_El_ElSigmaIetaIeta[i]=new TH1D(name, name, 50, -0.01, 0.08);
      hs_El_ElSigmaIetaIeta[i]->SetXTitle("#sigma_{i#eta i#eta}");
      hs_El_ElSigmaIetaIeta[i]->SetYTitle("events / bin");
      sprintf(name,"s_El_ElNumberOfMissingInnerHits_%d",i);
      hs_El_ElNumberOfMissingInnerHits[i]=new TH1D(name, name, 10, 0., 10.);
      hs_El_ElNumberOfMissingInnerHits[i]->SetXTitle("Number of missing inner hits");
      hs_El_ElNumberOfMissingInnerHits[i]->SetYTitle("events / bin");
    }

    TH1D *hc_El_counter;
    hc_El_counter=new TH1D("c_El_counter", "c_El_counter",15, 0, 15);
    hc_El_counter->GetXaxis()->SetBinLabel(1,"# e preselection");
    hc_El_counter->GetXaxis()->SetBinLabel(2,"# e presel, barrel");
    hc_El_counter->GetXaxis()->SetBinLabel(3,"# e see");
    hc_El_counter->GetXaxis()->SetBinLabel(4,"# e dphi");
    hc_El_counter->GetXaxis()->SetBinLabel(5,"# e deta");
    hc_El_counter->GetXaxis()->SetBinLabel(6,"# e hoe");
    hc_El_counter->GetXaxis()->SetBinLabel(7,"# e combIso");
    hc_El_counter->GetXaxis()->SetBinLabel(8,"# e conversion");
    hc_El_counter->GetXaxis()->SetBinLabel(9,"# e presel, endcap");
    hc_El_counter->GetXaxis()->SetBinLabel(10,"# e see");
    hc_El_counter->GetXaxis()->SetBinLabel(11,"# e dphi");
    hc_El_counter->GetXaxis()->SetBinLabel(12,"# e deta");
    hc_El_counter->GetXaxis()->SetBinLabel(13,"# e hoe");
    hc_El_counter->GetXaxis()->SetBinLabel(14,"# e combIso");
    hc_El_counter->GetXaxis()->SetBinLabel(15,"# e conversion");
    
    
    // --------------------------------------------------------------------------------------------------------------------


    // --------------------------------------------------------------------------------------------------------------------
    // Electron ID N-1 plots
    // --------------------------------------------------------------------------------------------------------------------
    int MAXN=4;
    TH1D *hn1_El_ElDeltaPhi[MAXN];
    TH1D *hn1_El_ElDeltaEta[MAXN];
    TH1D *hn1_El_ElOneOverEminusOneOverP[MAXN];
    TH1D *hn1_El_ElEthIso[MAXN];
    TH1D *hn1_El_ElCombIso[MAXN];
    TH1D *hn1_El_ElHoverE[MAXN];
    TH1D *hn1_El_ElSigmaIetaIeta[MAXN];
    TH1D *hn1_El_ElNumberOfMissingInnerHits[MAXN];
    

    for(int i=0; i<MAXN;i++) {
      sprintf(name,"n1_El_ElDeltaPhi_%d",i);
      hn1_El_ElDeltaPhi[i]=new TH1D(name, name, 30, -0.05001, 0.05001);    // plotting absolut value
//      hn1_El_ElDeltaPhi[i]=new TH1D(name, name, 15, -0.05001, 0.05001);     // plotting real value
      hn1_El_ElDeltaPhi[i]->SetXTitle("#Delta #phi");
      hn1_El_ElDeltaPhi[i]->SetYTitle("events / bin");
      sprintf(name,"n1_El_ElDeltaEta_%d",i);
      hn1_El_ElDeltaEta[i]=new TH1D(name, name, 30, -0.02001, 0.02001);     // plotting absolut value
//       hn1_El_ElDeltaEta[i]=new TH1D(name, name, 15, -0.05001, 0.05001);     // plotting real value
      hn1_El_ElDeltaEta[i]->SetXTitle("#Delta #eta");
      hn1_El_ElDeltaEta[i]->SetYTitle("events / bin");
      sprintf(name,"n1_El_ElOneOverEminusOneOverP_%d",i);
      hn1_El_ElOneOverEminusOneOverP[i]=new TH1D(name, name, 15, -0.0001, 0.05001);
      hn1_El_ElOneOverEminusOneOverP[i]->SetXTitle("|1/E-1/p|");
      hn1_El_ElOneOverEminusOneOverP[i]->SetYTitle("events / bin");
      sprintf(name,"n1_El_ElEthIso_%d",i);
      hn1_El_ElEthIso[i]=new TH1D(name, name, 50, -0.01, 5.);
      hn1_El_ElEthIso[i]->SetXTitle("Isolation");
      hn1_El_ElEthIso[i]->SetYTitle("events / bin");


      sprintf(name,"n1_El_ElCombIso_%d",i);
      hn1_El_ElCombIso[i]=new TH1D(name, name, 20, -0.01, 0.5);
//       hn1_El_ElCombIso[i]=new TH1D(name, name, 20, -0.01, 10.);
      hn1_El_ElCombIso[i]->SetXTitle("Combined Iso");
      hn1_El_ElCombIso[i]->SetYTitle("events / bin");

      sprintf(name,"n1_El_ElHoverE_%d",i);
      hn1_El_ElHoverE[i]=new TH1D(name, name, 30, -0.01, 0.2);
      hn1_El_ElHoverE[i]->SetXTitle("H/E");
      hn1_El_ElHoverE[i]->SetYTitle("events / bin");
      sprintf(name,"n1_El_ElSigmaIetaIeta_%d",i);
      hn1_El_ElSigmaIetaIeta[i]=new TH1D(name, name, 50, -0.01, 0.08);
      hn1_El_ElSigmaIetaIeta[i]->SetXTitle("#sigma_{i#eta i#eta}");
      hn1_El_ElSigmaIetaIeta[i]->SetYTitle("events / bin");
      sprintf(name,"n1_El_ElNumberOfMissingInnerHits_%d",i);
      hn1_El_ElNumberOfMissingInnerHits[i]=new TH1D(name, name, 10, 0., 10.);
      hn1_El_ElNumberOfMissingInnerHits[i]->SetXTitle("Number of missing inner hits");
      hn1_El_ElNumberOfMissingInnerHits[i]->SetYTitle("events / bin");


    }
    
    // --------------------------------------------------------------------------------------------------------------------
    

   // --------------------------------------------------------------------------------------------------------------------
   // W Selection plots
   // --------------------------------------------------------------------------------------------------------------------
   int MAXS=21;
   TH1D *hs_W_ElPt[MAXS];
   TH1D *hs_W_ElEta[MAXS];
   TH1D *hs_W_ElPhi[MAXS];
   TH1D *hs_W_ElDeltaPhi[MAXS];
   TH1D *hs_W_ElMet[MAXS];
   TH1D *hs_W_ElTransMass[MAXS];
   TH1D *hs_W_ElCharge[MAXS];
   TH1D *hs_W_Njets[MAXS];
   TH1D *hs_W_jetPt[MAXS];
    for(int i=0; i<MAXS; i++) {
      sprintf(name,"s_W_ElPt_%d",i);
      hs_W_ElPt[i]=new TH1D(name, name, 60, 0, 120);
      hs_W_ElPt[i]->SetXTitle("p_{T} [GeV]");
      hs_W_ElPt[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_W_ElEta_%d",i);
      hs_W_ElEta[i]=new TH1D(name, name, 50, -5., 5.);
      hs_W_ElEta[i]->SetXTitle("#eta");
      hs_W_ElEta[i]->SetYTitle("events / 0.2");
      sprintf(name,"s_W_ElPhi_%d",i);
      hs_W_ElPhi[i]=new TH1D(name, name, 30, -Pi,Pi);
      hs_W_ElPhi[i]->SetXTitle("#phi");
      hs_W_ElPhi[i]->SetYTitle("events / 0.21");
      sprintf(name,"s_W_ElDeltaPhi_%d",i);
      hs_W_ElDeltaPhi[i]=new TH1D(name, name, 15, 0., Pi);
      hs_W_ElDeltaPhi[i]->SetXTitle("#Delta #phi_{e,#slash{E}_{T}}");
      hs_W_ElDeltaPhi[i]->SetYTitle("events / bin");
      sprintf(name,"s_W_ElMet_%d",i);
      hs_W_ElMet[i]=new TH1D(name, name, 60, 0, 120);
      hs_W_ElMet[i]->SetXTitle("#slash{E}_{T} [GeV]");
      hs_W_ElMet[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_W_ElTransMass_%d",i);
      hs_W_ElTransMass[i]=new TH1D(name, name, 100, 0, 200);
      hs_W_ElTransMass[i]->SetXTitle("m_{T} [GeV]");
      hs_W_ElTransMass[i]->SetYTitle("events / 2 GeV");
      sprintf(name,"s_W_ElCharge_%d",i);
      hs_W_ElCharge[i]=new TH1D(name, name, 2, -2, 2);
      hs_W_ElCharge[i]->SetXTitle("Charge");
      hs_W_ElCharge[i]->SetYTitle("events / bin");
      sprintf(name,"s_W_Njets_%d",i);
      hs_W_Njets[i]=new TH1D(name, name, 10, 0, 10);
      hs_W_Njets[i]->SetXTitle("number of jets");
      hs_W_Njets[i]->SetYTitle("events / bin");
      sprintf(name,"s_W_jetPt_%d",i);
      hs_W_jetPt[i]=new TH1D(name, name, 40, 0, 200);
      hs_W_jetPt[i]->SetXTitle("p_{T} [GeV]");
      hs_W_jetPt[i]->SetYTitle("events / bin");
    }

    int MAXWcounter=16;
    int *Wcounter = new int[MAXWcounter];
    memset(Wcounter,0,MAXWcounter*sizeof(int));


    TH1D *hc_W_counter;
    hc_W_counter=new TH1D("c_W_counter", "c_W_counter",16, 0, 16 );
    hc_W_counter->GetXaxis()->SetBinLabel(1,"# e, N jets, Preselection");
    hc_W_counter->GetXaxis()->SetBinLabel(2,"# W, N jets");
    hc_W_counter->GetXaxis()->SetBinLabel(3,"# W+, N jets");
    hc_W_counter->GetXaxis()->SetBinLabel(4,"# W-, N jets");
    hc_W_counter->GetXaxis()->SetBinLabel(5,"# e, 0 jets, Preselection");
    hc_W_counter->GetXaxis()->SetBinLabel(6,"# W, 0 jets");
    hc_W_counter->GetXaxis()->SetBinLabel(7,"# W+, 0 jets");
    hc_W_counter->GetXaxis()->SetBinLabel(8,"# W-, 0 jets");
    hc_W_counter->GetXaxis()->SetBinLabel(9,"# e, 1 jet, Preselection");
    hc_W_counter->GetXaxis()->SetBinLabel(10,"# W, 1 jet");
    hc_W_counter->GetXaxis()->SetBinLabel(11,"# W+, 1 jet");
    hc_W_counter->GetXaxis()->SetBinLabel(12,"# W-, 1 jet");
    hc_W_counter->GetXaxis()->SetBinLabel(13,"# e, >1 jets, Preselection");
    hc_W_counter->GetXaxis()->SetBinLabel(14,"# W, >1 jets");
    hc_W_counter->GetXaxis()->SetBinLabel(15,"# W+, >1 jets");
    hc_W_counter->GetXaxis()->SetBinLabel(16,"# W-, >1jets");

   // --------------------------------------------------------------------------------------------------------------------


    // --------------------------------------------------------------------------------------------------------------------
    // Z plots 
    // --------------------------------------------------------------------------------------------------------------------    
    int MAXZ=1;
    TH1D *h_Z_invMass[MAXZ];
    TH1D *h_Z_invMass2[MAXZ];
    TH1D *h_Z_ptAll[MAXZ];
    TH1D *h_Z_etaAll[MAXZ];
    TH1D *h_Z_phiAll[MAXZ];
    TH1D *h_Z_deltaPhi[MAXZ];
    TH1D *h_Z_ptDiLepton[MAXZ];
    TH1D *h_Z_ptDiLepton2[MAXZ];
    TH1D *h_Z_etaDiLepton[MAXZ];
    TH1D *h_Z_phiDiLepton[MAXZ];

    TH1D *h_Z_leadingJetPt[MAXZ];
    TH1D *h_Z_summedJetPt[MAXZ];
    TH1D *h_Z_leadingJetEta[MAXZ];
    TH1D *h_Z_summedJetEta[MAXZ];
    TH1D *h_Z_leadingJetPhi[MAXZ];
    TH1D *h_Z_summedJetPhi[MAXZ];
    TH1D *h_Z_nJets[MAXZ];
    TH1D *h_Z_leadingJetDeltaPt[MAXZ];
    TH1D *h_Z_leadingJetDeltaEta[MAXZ];
    TH1D *h_Z_leadingJetDeltaPhi[MAXZ];
    TH1D *h_Z_summedJetDeltaPt[MAXZ];
    TH1D *h_Z_summedJetDeltaEta[MAXZ];
    TH1D *h_Z_summedJetDeltaPhi[MAXZ];
    TH1D *h_Z_recoilPt[MAXZ];
    TH1D *h_Z_recoilEta[MAXZ];
    TH1D *h_Z_recoilPhi[MAXZ];
    TH1D *h_Z_metPt[MAXZ];
    TH1D *h_Z_metPhi[MAXZ];


    for(int i=0; i<MAXZ; i++) {
      sprintf(name,"h_Z_invMass_%d",i);
      h_Z_invMass[i] = new TH1D(name,name,100,0,200);
      h_Z_invMass[i]->SetXTitle("m_{ee} [GeV]");
      h_Z_invMass[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_invMass2_%d",i);
      h_Z_invMass2[i] = new TH1D(name,name,30,60,120);
      h_Z_invMass2[i]->SetXTitle("m_{ee} [GeV]");
      h_Z_invMass2[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_ptAll_%d",i);
      h_Z_ptAll[i] = new TH1D(name,name,20,0,200);
      h_Z_ptAll[i]->SetXTitle("p_{T} [GeV]");
      h_Z_ptAll[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_etaAll_%d",i);
      h_Z_etaAll[i] = new TH1D(name,name,12,-3.,3.);
      h_Z_etaAll[i]->SetXTitle("#eta");
      h_Z_etaAll[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_phiAll_%d",i);
      h_Z_phiAll[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_phiAll[i]->SetXTitle("#phi");
      h_Z_phiAll[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_deltaPhi_%d",i);
      h_Z_deltaPhi[i] = new TH1D(name,name,15,0,Pi);
      h_Z_deltaPhi[i]->SetXTitle("#Delta #phi_{ee}");
      h_Z_deltaPhi[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_ptDiLepton_%d",i);
      h_Z_ptDiLepton[i] = new TH1D(name,name,40,0,200);
      h_Z_ptDiLepton[i]->SetXTitle("p_{T} [GeV]");
      h_Z_ptDiLepton[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_ptDiLepton2_%d",i);
      h_Z_ptDiLepton2[i] = new TH1D(name,name,50,0,50);
      h_Z_ptDiLepton2[i]->SetXTitle("p_{T} [GeV]");
      h_Z_ptDiLepton2[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_etaDiLepton_%d",i);
      h_Z_etaDiLepton[i] = new TH1D(name,name,16,-8.,8.);
      h_Z_etaDiLepton[i]->SetXTitle("#eta");
      h_Z_etaDiLepton[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_phiDiLepton_%d",i);
      h_Z_phiDiLepton[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_phiDiLepton[i]->SetXTitle("#phi");
      h_Z_phiDiLepton[i]->SetYTitle("events / bin");


      sprintf(name,"h_Z_leadingJetPt_%d",i);
      h_Z_leadingJetPt[i] = new TH1D(name,name,20,0,200);
      h_Z_leadingJetPt[i]->SetXTitle("p_{T} [GeV]");
      h_Z_leadingJetPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_summedJetPt_%d",i);
      h_Z_summedJetPt[i] = new TH1D(name,name,20,0,200);
      h_Z_summedJetPt[i]->SetXTitle("p_{T} [GeV]");
      h_Z_summedJetPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_leadingJetEta_%d",i);
      h_Z_leadingJetEta[i] = new TH1D(name,name,12,-3,3);
      h_Z_leadingJetEta[i]->SetXTitle("#eta");
      h_Z_leadingJetEta[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_summedJetEta_%d",i);
      h_Z_summedJetEta[i] = new TH1D(name,name,16,-8,8);
      h_Z_summedJetEta[i]->SetXTitle("#eta");
      h_Z_summedJetEta[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_leadingJetPhi_%d",i);
      h_Z_leadingJetPhi[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_leadingJetPhi[i]->SetXTitle("#phi");
      h_Z_leadingJetPhi[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_summedJetPhi_%d",i);
      h_Z_summedJetPhi[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_summedJetPhi[i]->SetXTitle("#phi");
      h_Z_summedJetPhi[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_nJets_%d",i);
      h_Z_nJets[i] = new TH1D(name,name,11,0,11);
      h_Z_nJets[i]->SetXTitle("jet number");
      h_Z_nJets[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_leadingJetDeltaPt_%d",i);
      h_Z_leadingJetDeltaPt[i] = new TH1D(name,name,60,-150,150);
      h_Z_leadingJetDeltaPt[i]->SetXTitle("#Delta p_{T,jZ} [GeV]");
      h_Z_leadingJetDeltaPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_leadingJetDeltaEta_%d",i);
      h_Z_leadingJetDeltaEta[i] = new TH1D(name,name,16,-8,8);
      h_Z_leadingJetDeltaEta[i]->SetXTitle("#Delta #eta_{jZ}");
      h_Z_leadingJetDeltaEta[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_leadingJetDeltaPhi_%d",i);
      h_Z_leadingJetDeltaPhi[i] = new TH1D(name,name,15,0,Pi);
      h_Z_leadingJetDeltaPhi[i]->SetXTitle("#Delta #phi_{jZ}");
      h_Z_leadingJetDeltaPhi[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_summedJetDeltaPt_%d",i);
      h_Z_summedJetDeltaPt[i] = new TH1D(name,name,60,-150,150);
      h_Z_summedJetDeltaPt[i]->SetXTitle("#Delta p_{T,jZ} [GeV]");
      h_Z_summedJetDeltaPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_summedJetDeltaEta_%d",i);
      h_Z_summedJetDeltaEta[i] = new TH1D(name,name,16,-8,8);
      h_Z_summedJetDeltaEta[i]->SetXTitle("#Delta #eta_{jZ}");
      h_Z_summedJetDeltaEta[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_summedJetDeltaPhi_%d",i);
      h_Z_summedJetDeltaPhi[i] = new TH1D(name,name,15,0,Pi);
      h_Z_summedJetDeltaPhi[i]->SetXTitle("#Delta #phi_{jZ}");
      h_Z_summedJetDeltaPhi[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_recoilPt_%d",i);
      h_Z_recoilPt[i] = new TH1D(name,name,20,0,200);
      h_Z_recoilPt[i]->SetXTitle("p_{T,j+Z} [GeV]");
      h_Z_recoilPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_recoilEta_%d",i);
      h_Z_recoilEta[i] = new TH1D(name,name,16,-8,8);
      h_Z_recoilEta[i]->SetXTitle("#eta_{j+Z}");
      h_Z_recoilEta[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_recoilPhi_%d",i);
      h_Z_recoilPhi[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_recoilPhi[i]->SetXTitle("#phi_{j+Z}");
      h_Z_recoilPhi[i]->SetYTitle("events / bin");

      sprintf(name,"h_Z_metPt_%d",i);
      h_Z_metPt[i] = new TH1D(name,name,20,0,200);
      h_Z_metPt[i]->SetXTitle("#slash{E}_{T} [GeV]");
      h_Z_metPt[i]->SetYTitle("events / bin");
      sprintf(name,"h_Z_metPhi_%d",i);
      h_Z_metPhi[i] = new TH1D(name,name,30,-Pi,Pi);
      h_Z_metPhi[i]->SetXTitle("#phi_{#slash{E}_{T}}");
      h_Z_metPhi[i]->SetYTitle("events / bin");
      
//       TH2D *h_Z_summedJetDeltaPhi_VS_deltaPhi("summedJetDeltaPhi_VS_deltaPhi", "summedJetDeltaPhi_VS_deltaPhi", 30,-Pi,Pi,15,0,Pi);

    }
      // --------------------------------------------------------------------------------------------------------------------


    ////////////////////////////////////////////////////////////////////
    ////////////////////Loop over entries starts here //////////////////
    ////////////////////////////////////////////////////////////////////
    

    TH1D *h_jetPt = new TH1D("h_jetPt", "jet_Pt",100,0,100);
    
    /////////initialize variables
    
    Long64_t nentries = fChain->GetEntriesFast();
    std::cout << "analysisClass::Loop(): nentries = " << nentries << std::endl;   
    
    ////// The following ~7 lines have been taken from rootNtupleClass->Loop() /////
    ////// If the root version is updated and rootNtupleClass regenerated,     /////
    ////// these lines may need to be updated.                                 /////    
    Long64_t nb = 0;
    

    for (Long64_t jentry=0; jentry<nentries;jentry++) {  // loop over all events (entries) in the tree/chain


      Long64_t ientry = fChain->LoadTree(jentry);  // ientry: this is a counter over each single file in the input list

      if (ientry < 0) {cerr<<"WARNING: ientry = "<<ientry<<" is smaller than 0! (this is ok if the end of the file is reached...)"<<endl; break; }
      
      nb=fChain->GetEntry(jentry);
      //fChain->GetEntry(jentry);
      
      if(jentry < 10 || jentry%1000 == 0) {
        std::cout << "analysisClass::Loop(): Events processed = " << jentry << std::endl;   
      }
      
      if(sampleSize!=-1 && sampleSize==jentry) {
        goto END;
      }
      h_numberOfEvents->Fill(1);
      
     ////////////////////////////////////////////////////////////////////
     ////////////////////// User's code starts here /////////////////////
     ////////////////////////////////////////////////////////////////////
     
     // check the methods defined in baseClass and clean it up!!
     

      // --------------------------------------------------------------------------------------------------------------------
      // the MC samples have always run number ==1,
      // exploit this fact to set the isdata flag
      // --------------------------------------------------------------------------------------------------------------------
      if(firstEvent) {
        firstEvent=false;
        cout<<"check run number "<<Run<<endl;
        if(Run!=1) isdata=true;
        else if(Run==1) isdata=false;
        if(isdata) cout<<"Running over real data!!!"<<endl;
        else if(!isdata) cout<<"Running over MC!!!"<<endl;
      }
      // --------------------------------------------------------------------------------------------------------------------
      

      // check for every event if a new run has started and get the trigger info
      if(Run != previousRun) {
        previousRun = Run;
        cout<<"Now processing run number: "<<Run<<endl;
        GetHLTNames(Run);
        cout<<"got them..."<<endl;
      }
      

      i_preselected.clear();
      i_isolated.clear();
      i_identified.clear();
      i_goodElectrons.clear();
      i_goodBarrelElectrons.clear();
      i_goodEndcapElectrons.clear();
      i_jets.clear();

      // reset the branch arrays
      tRun         = -999;
      tEvent       = -999;
      tLumiSection = -999;

      tElN = 0;
      memset(tElPt,     -999.999, ElMax*sizeof(double));
      memset(tElPx,     -999.999, ElMax*sizeof(double));
      memset(tElPy,     -999.999, ElMax*sizeof(double));
      memset(tElPz,     -999.999, ElMax*sizeof(double));
      memset(tElE,    -999.999, ElMax*sizeof(double));
      memset(tElEta,    -999.999, ElMax*sizeof(double));
      memset(tElPhi,    -999.999, ElMax*sizeof(double));
      memset(tElCharge, -999,     ElMax*sizeof(int));
      memset(tElDz, -999,     ElMax*sizeof(double));

      tPfMet   = -999.999;
      tCaloMet = -999.999;
      tTcMet   = -999.999;

      memset(tElPfMt,   -999.999, ElMax*sizeof(double));
      memset(tElCaloMt, -999.999, ElMax*sizeof(double));
      memset(tElTcMt,   -999.999, ElMax*sizeof(double));

      tSumEHF_minus  = -999.999;
      tSumEHF_plus   = -999.999;
      tSumEHF_Aminus  = -999.999;
      tSumEHF_Aplus   = -999.999;
      tSumEHF_Bminus  = -999.999;
      tSumEHF_Bplus   = -999.999;
      tSumETHF_minus = -999.999;
      tSumETHF_plus  = -999.999;
      tNtowers_plus  = 0.;
      tNtowers_minus  = 0.;

      tCaloTowersN = 0;
      memset(tCaloTowersE,     -999.999, TowMax*sizeof(double));
      memset(tCaloTowersEt,     -999.999, TowMax*sizeof(double));
      memset(tCaloTowersEta,     -999.999, TowMax*sizeof(double));
      memset(tCaloTowersEmE,   -999.999, TowMax*sizeof(double));
      memset(tCaloTowersHadE,  -999.999, TowMax*sizeof(double));
      memset(tCaloTowersHasEB, -999.999, TowMax*sizeof(int));
      memset(tCaloTowersHasEE, -999.999, TowMax*sizeof(int));
      memset(tCaloTowersHasHB, -999.999, TowMax*sizeof(int));
      memset(tCaloTowersHasHE, -999.999, TowMax*sizeof(int));
      memset(tCaloTowersHasHF, -999.999, TowMax*sizeof(int));
      memset(tCaloTowersId,    -999.999, TowMax*sizeof(int));


      tJetN = 0;
      memset(tJetPt,   -999.999, ElMax*sizeof(double));
      memset(tJetEta,  -999.999, ElMax*sizeof(double));
      memset(tJetPhi,  -999.999, ElMax*sizeof(double));

      tEtaMax   = -9.9;
      tEtaMin   =  9.9;
      tEminusPz = 0.;
      tEplusPz  = 0.;
      
      tSumE  = 0.;
      tSumPz = 0.;

      // pflow
      tPfN = 0;
      memset(tPfE,    -999.999, 1000*sizeof(double));
      memset(tPfEt,   -999.999, 1000*sizeof(double));
      memset(tPfEta,  -999.999, 1000*sizeof(double));
      memset(tPfPt,   -999.999, 1000*sizeof(double));


      tVertexN = 0;
//       memset(tVertexNdof,   -999.999, 25*sizeof(int));
      memset(tVertexNdof,   -999.999, 25*sizeof(double));
      memset(tVertexChi2,   -999.999, 25*sizeof(double));
      memset(tVertexX,      -999.999, 25*sizeof(double));
      memset(tVertexY,      -999.999, 25*sizeof(double));
      memset(tVertexZ,      -999.999, 25*sizeof(double));
      // reset the vertex tracks array
      // FIXME: hardcoded array size!!!
      for (int ivtx=0;ivtx<25;++ivtx)
        {
          for (int itrk=0;itrk<400;++itrk)
            {
              tVertexTrackIndex[ivtx][itrk]=-1;
//            fTVertexTrackWeight[ivtx][itrk]=-1;
            }
    }
      
      //pile-up
      tPUnumInteractions = 0;
      memset(tPUzPositions,   -999.999, 15*sizeof(double));

      tBeamSpotX = -999.999;
      tBeamSpotY = -999.999;
      tBeamSpotZ = -999.999;

      tTracksN = 0;
      memset(tTrackPt,      -999.999, TrackMax*sizeof(double));
      memset(tTrackPx,      -999.999, TrackMax*sizeof(double));
      memset(tTrackPy,      -999.999, TrackMax*sizeof(double));
      memset(tTrackPz,      -999.999, TrackMax*sizeof(double));
      memset(tTrackDzBs,      -999.999, TrackMax*sizeof(double));
      memset(tTrackEta,      -999.999, TrackMax*sizeof(double));
      memset(tTrackPhi,      -999.999, TrackMax*sizeof(double));
//       memset(tTrackdsz,      -999.999, 1000*sizeof(double));
//       memset(tTrackdxy,      -999.999, 1000*sizeof(double));


     // --------------------------------------------------------------------------------
     // get the MET
     // --------------------------------------------------------------------------------
     double met;
     double metPhi;
     double delPhiMet;
     met=PFMET;
     metPhi=PFMETphi;

     // --------------------------------------------------------------------------------
     // get the muons
     // --------------------------------------------------------------------------------
//      for(int m=0; m<NMus; m++) {
       
//        cout<<MuPt[m]<<endl;
       
//      }

     // --------------------------------------------------------------------------------


     // --------------------------------------------------------------------------------
     // loop over all electrons in one event
     // --------------------------------------------------------------------------------
     for(int i=0; i<NEles; i++) {
       
       // --------------------------------------------------------------------------------
       // get the indices of the electrons in each event
       // --------------------------------------------------------------------------------
       bool pre=false;
       bool iso=false;
       bool id=false;

       if(ElectronPreselection(i)==0) {i_preselected.push_back(i); debug++; pre=true;}

       if(!ethIDandIso) {

//        if(WPElectronID(i, iso, WorkingPointCuts)==0) {i_isolated.push_back(); iso=true;}
//        else continue;
         if(WPElectronID(i, all, WorkingPointCuts)==0) {i_identified.push_back(i); debug2++; id=true;}
//       if(ElIDsimpleWP80relIso[i] == 7) {i_identified.push_back(i); debug2++; id=true;}
         
         cout<<WPElectronID(i, all, WorkingPointCuts)<<endl;

         if(ElIDsimpleWP80relIso[i] == 5) cout<<WPElectronID(i, all, WorkingPointCuts)<<endl;

         if(id==true) cout<<ElPt[i]<<endl;


         if(pre==true && id==true && (fabs(ElSCEta[i]) < 1.4442 || fabs(ElSCEta[i]) > 1.566) && fabs(ElSCEta[i]) < 2.5) {i_goodElectrons.push_back(i);}
         // fiducial cuts
         if(fabs(ElSCEta[i]) < 1.4442 && pre==true && id==true)                          {i_goodBarrelElectrons.push_back(i);}
         if(fabs(ElSCEta[i]) > 1.566 && fabs(ElSCEta[i]) < 2.5 && pre==true && id==true) {i_goodEndcapElectrons.push_back(i);}
       }

       if(ethIDandIso) {
         if(ETHElectronID(i, all)) {i_identified.push_back(i); debug2++; id=true;}
         if(fabs(ElSCEta[i]) < 1.4442 && pre==true && id==true)                          {i_goodBarrelElectrons.push_back(i);}
         if(fabs(ElSCEta[i]) > 1.566 && fabs(ElSCEta[i]) < 2.5 && pre==true && id==true) {i_goodEndcapElectrons.push_back(i);}

       }
       // --------------------------------------------------------------------------------
       
       
           // the N-1 plots
//         if(fabs(accEta)<etaBarrelMax) {
//           if(WPElectronID(i, see, WorkingPointCuts)==0)  hn1_El_ElSigmaIetaIeta[2]->Fill(See);
//           if(WPElectronID(i, dphi, WorkingPointCuts)==0) hn1_El_ElDeltaPhi[2]->Fill(delPhi);
//           if(WPElectronID(i, deta, WorkingPointCuts)==0) hn1_El_ElDeltaEta[2]->Fill(delEta);
//           if(WPElectronID(i, hoe, WorkingPointCuts)==0)  hn1_El_ElHoverE[2]->Fill(Hoe);
//           if(WPElectronID(i, combiso, WorkingPointCuts)==0)    hn1_El_ElCombIso[2]->Fill(combIso_B);
//           if(WPElectronID(i, conversion, WorkingPointCuts)==0) hn1_El_ElNumberOfMissingInnerHits[2]->Fill(missingHits);
//         }
//         if(fabs(accEta)>etaEndcapMin && fabs(accEta)<etaEndcapMax) {
//           if(WPElectronID(i, see, WorkingPointCuts)==0)  hn1_El_ElSigmaIetaIeta[3]->Fill(See);
//           if(WPElectronID(i, dphi, WorkingPointCuts)==0) hn1_El_ElDeltaPhi[3]->Fill(delPhi);
//           if(WPElectronID(i, deta, WorkingPointCuts)==0) hn1_El_ElDeltaEta[3]->Fill(delEta);
//           if(WPElectronID(i, hoe, WorkingPointCuts)==0)  hn1_El_ElHoverE[3]->Fill(Hoe);
//           if(WPElectronID(i, combiso, WorkingPointCuts)==0)    hn1_El_ElCombIso[3]->Fill(combIso_E);
//           if(WPElectronID(i, conversion, WorkingPointCuts)==0) hn1_El_ElNumberOfMissingInnerHits[3]->Fill(missingHits);
//         } 
           

     } // loop over all electrons in one event

     
     // ******************************************************************************
     // this is the actual code
     // ******************************************************************************

     double etaMax_B = 1.4442;
     double etaMin_E = 1.566;
     double etaMax_E = 2.5;

     //--------------------------------------------------------------------------------------------------------------------
     // check fiducial region of the preselected electrons
     //--------------------------------------------------------------------------------------------------------------------
     for(int i=0; i<i_preselected.size(); i++) {
       if(fabs(ElSCEta[i_preselected[i]]) < 1.4442) {
         debug3++;

       }
     }
     //--------------------------------------------------------------------------------------------------------------------


     //--------------------------------------------------------------------------------------------------------------------
     // HLT Trigger
     //--------------------------------------------------------------------------------------------------------------------
     bool isTriggered = false;
     if(isdata && useElectronTiggers) {
       int Photon15_L1R         = getHLTtriggerBit("HLT_Photon15_L1R");
       int Photon15_Cleaned_L1R = getHLTtriggerBit("HLT_Photon15_Cleaned_L1R");
       int Photon20_Cleaned_L1R = getHLTtriggerBit("HLT_Photon20_Cleaned_L1R");
       int Ele10_LW_L1R         = getHLTtriggerBit("HLT_Ele10_LW_L1R");
       int Ele10_SW_L1R         = getHLTtriggerBit("HLT_Ele10_SW_L1R");
       int Ele15_SW_L1R         = getHLTtriggerBit("HLT_Ele15_SW_L1R");
       int Ele15_SW_CaloEleId_L1R  = getHLTtriggerBit("HLT_Ele15_SW_CaloEleId_L1R");
       int Ele15_SW_TightEleId_L1R = getHLTtriggerBit("HLT_Ele15_SW_TightEleId_L1R");
       int Ele17_SW_CaloEleId_L1R  = getHLTtriggerBit("HLT_Ele17_SW_CaloEleId_L1R");
       int Ele17_SW_EleId_L1R      = getHLTtriggerBit("HLT_Ele17_SW_EleId_L1R");
       int Ele17_SW_TightEleId_L1R = getHLTtriggerBit("HLT_Ele17_SW_TightEleId_L1R");
       int Ele17_SW_TighterEleId_L1R_v1 = getHLTtriggerBit("HLT_Ele17_SW_TighterEleId_L1R_v1");
       int Ele17_SW_TighterEleId_L1R_v2 = getHLTtriggerBit("HLT_Ele17_SW_TighterEleId_L1R_v2");
       int Ele17_SW_TighterEleId_L1R_v3 = getHLTtriggerBit("HLT_Ele17_SW_TighterEleId_L1R_v3");

//        check if at least one trigger was found
         if(Photon15_L1R == -1 &&
            Photon15_Cleaned_L1R == -1 &&
            Photon20_Cleaned_L1R == -1 &&
            Ele10_LW_L1R == -1 &&
            Ele10_SW_L1R == -1 &&
            Ele15_SW_L1R == -1 &&
            Ele15_SW_CaloEleId_L1R  == -1 &&
            Ele15_SW_TightEleId_L1R == -1 &&
            Ele17_SW_CaloEleId_L1R  == -1 && 
            Ele17_SW_EleId_L1R      == -1 &&
            Ele17_SW_TightEleId_L1R == -1 &&
            Ele17_SW_TighterEleId_L1R_v1 == -1 &&
            Ele17_SW_TighterEleId_L1R_v2 == -1 &&
            Ele17_SW_TighterEleId_L1R_v3 == -1
            ) {
           cerr<<"None of the trigger bits found!!!"<<endl;
           h_HLTInfo->Fill(0);
         }
       
//        if(isdata  && Run <  141956 && HLTResults[Ele15_LW] == 1)           isTriggered=true;
//        if(isdata  && Run >= 141956 && HLTResults[Ele15_SW_CaloEleId] == 1) isTriggered=true;
       //cout<<"HLT result: "<<HLTResults[Ele17_SW_CaloEleId]<<endl;
       //        if(isdata &&                   HLTResults[Ele17_SW_CaloEleId]==1)   isTriggered=true;     
       if(isdata && (
                     HLTResults[Photon15_L1R]            == 1 ||
                     HLTResults[Photon15_Cleaned_L1R]    == 1 ||
                     HLTResults[Photon20_Cleaned_L1R]    == 1 ||
                     HLTResults[Ele10_LW_L1R]            == 1 ||
                     HLTResults[Ele10_SW_L1R]            == 1 ||
                     HLTResults[Ele15_SW_L1R]            == 1 ||
                     HLTResults[Ele15_SW_CaloEleId_L1R]  == 1 ||
                     HLTResults[Ele15_SW_TightEleId_L1R] == 1 ||
                     HLTResults[Ele17_SW_CaloEleId_L1R]  == 1 ||
                     HLTResults[Ele17_SW_EleId_L1R]      == 1 ||
                     HLTResults[Ele17_SW_TightEleId_L1R] == 1 ||
                     HLTResults[Ele17_SW_TighterEleId_L1R_v1] == 1 ||
                     HLTResults[Ele17_SW_TighterEleId_L1R_v2] == 1 ||
                     HLTResults[Ele17_SW_TighterEleId_L1R_v3] == 1
                     )) isTriggered=true;     
     }

     else if(isdata && !useElectronTiggers) {
       int HLT_ZeroBias        = getHLTtriggerBit("HLT_ZeroBias");
       if(isdata && (
                     HLTResults[HLT_ZeroBias]==1
                     )) isTriggered=true;     
     }
     
     else if(!isdata) {
       int Ele15           = getHLTtriggerBit("HLT_Ele15_LW_L1R");
       //cout<<"HLT result: "<<HLTResults[ Ele15]<<endl;
       if(HLTResults[Ele15]==1)             isTriggered=true;
     }
//      //if(HLTResults(Ele15_SW)==1)                              isTriggered=true;
// //      if(isdata  && Run <  141956 && HLTResults[Ele15_LW] == 1)           isTriggered=true;
// //      if(isdata  && Run >= 141956 && HLTResults[Ele15_SW_CaloEleId] == 1) isTriggered=true;
// //      if(!isdata &&                  HLTResults[Ele15_SW]==1)             isTriggered=true;

       isTriggered=true;

//      cout<<"here0"<<endl;


     //--------------------------------------------------------------------------------------------------------------------
     
     //--------------------------------------------------------------------------------------------------------------------
     // get the number of jets per event */*
     //--------------------------------------------------------------------------------------------------------------------
     int nJets=0;
     for(int j=0; j<PFNJets; j++) {
       double scale = 1; //JScale[j];
       double jpx = PFJPx[j] / scale;
       double jpy = PFJPy[j] / scale;
       double jpz = PFJPz[j] / scale;
       double jpt = PFJPt[j] / scale;
       double jeta = PFJEta[j];
       TVector3 jet;
       jet.SetXYZ(jpx, jpy, jpz);
       double deltar;
       bool notMatched=true;
       // match jet with an electron
       for(int e=0; e<i_goodElectrons.size(); e++) {
         TVector3 el;
         el.SetXYZ(ElPx[i_goodElectrons[e]], ElPy[i_goodElectrons[e]], ElPz[i_goodElectrons[e]]);
         deltar = jet.DeltaR(el);
         if(fabs(deltar) < 0.3) notMatched=false;
       }
       // apply loose jet ID
       bool jetID=false;
       if(JEMfrac[j] > 0.01 &&
          JID_HPD[j] < 0.989 &&
          JID_n90Hits[j] > 1) jetID=true;
       // check if matched and good jet
       if(notMatched && jpt > 30 && fabs(jeta) < 2.5 && jetID) {
         nJets++;
         i_jets.push_back(j);
       }
     } // end: number of jets

     
     //--------------------------------------------------------------------------------------------------------------------

     //--------------------------------------------------------------------------------------------------------------------
     // check if diffractive event
     //--------------------------------------------------------------------------------------------------------------------
     bool isDiffractive=false;
//      if(SumEHF_minus == 0  || SumEHF_plus == 0 ||
//      SumETHF_minus == 0 || SumETHF_plus == 0 ) {
//        isDiffractive=true;
//      }
     
     //--------------------------------------------------------------------------------------------------------------------


     //--------------------------------------------------------------------------------------------------------------------
     //--------------------------------------------------------------------------------------------------------------------
     // ELECTRONS and general things
     //--------------------------------------------------------------------------------------------------------------------
     //--------------------------------------------------------------------------------------------------------------------

     // count total number of events
     El_count[0]++;

     if(!isTriggered) continue;
     
     // count triggered events
     El_count[1]++;

//      cout<<"bla"<<endl;

//      if(i_preselected.size() > 0) {
//        tRun = Run;
//        tEvent = Event;
//        tLumiSection = LumiSection;

//      }

//      cout<<tRun<<endl;
//      cout<<tEvent<<endl;
//      cout<<tLumiSection<<endl;


//      cout<<"here1"<<endl;


     // look at preselected electrons
     for(int i=0; i< i_preselected.size(); i++) {
       // calculate the variables to fill
       double deltaPhiElMet=deltaPhi(ElPhi[i_preselected[i]], PFMETphi);
       double mt = transverseMass(ElPt[i_preselected[i]] * PtCorrection, met, deltaPhiElMet);
       int charge = ElCharge[i_preselected[i]];
       int dz = ElDzBS[i_preselected[i]];
       double pt = ElPt[i_preselected[i]] * PtCorrection;
       double phi = ElPhi[i_preselected[i]];
       double eta = ElEta[i_preselected[i]];
       double Hoe=ElHcalOverEcal[i_preselected[i]];
       double trkIso  = ElDR03TkSumPt[i_preselected[i]];
       double ecalIso = ElDR03EcalRecHitSumEt[i_preselected[i]];
       double hcalIso = ElDR03HcalTowerSumEt[i_preselected[i]];
       double combIso_B=(trkIso + max(0., ecalIso - 1.) + hcalIso) / ( ElPt[i_preselected[i]] * PtCorrection);
       double combIso_E=(trkIso + ecalIso + hcalIso) / (ElPt[i_preselected[i]] * PtCorrection);
       double See=ElSigmaIetaIeta[i_preselected[i]];
       double missingHits=ElNumberOfMissingInnerHits[i_preselected[i]];
       double dEtaCorrection=0;
       double dPhiCorrection=0;
       if(isdata==true && applyDeltaCorrection==true) {
         dEtaCorrection=detaCorrections(ElSCEta[i_preselected[i]], ElPhi[i_preselected[i]]);
         dPhiCorrection=dphiCorrections(ElSCEta[i_preselected[i]], ElPhi[i_preselected[i]]);
       }
       double delPhi=ElDeltaPhiSuperClusterAtVtx[i_preselected[i]]-dphiCorrections(ElEta[i_preselected[i]], ElPhi[i_preselected[i]]);
       double delEta=ElDeltaEtaSuperClusterAtVtx[i_preselected[i]]-detaCorrections(ElEta[i_preselected[i]], ElPhi[i_preselected[i]]);          
       
       // barrel
       if(fabs(ElSCEta[i_preselected[i]]) < etaMax_B) {
         hs_El_ElTransMass[0]->Fill(mt);
         hs_El_ElPt[0]->Fill(pt);
         hs_El_ElPhi[0]->Fill(phi);
         hs_El_ElEta[0]->Fill(eta);
         hs_El_ElMet[0]->Fill(met);
         hs_El_ElCharge[0]->Fill(charge);
         hs_El_Njets[0]->Fill(nJets);
         El_count[2]++;
       }
       // endcap
       if(fabs(ElSCEta[i_preselected[i]]) > etaMin_E && fabs(ElSCEta[i_preselected[i]]) < etaMax_E) {
         hs_El_ElTransMass[1]->Fill(mt);
         hs_El_ElPt[1]->Fill(pt);
         hs_El_ElPhi[1]->Fill(phi);
         hs_El_ElEta[1]->Fill(eta);
         hs_El_ElMet[1]->Fill(met);
         hs_El_ElCharge[1]->Fill(charge);
         hs_El_Njets[1]->Fill(nJets);
         El_count[3]++;
       }
       // barrel + endcap
       if(fabs(ElSCEta[i_preselected[i]]) < etaMax_B || (fabs(ElSCEta[i_preselected[i]]) > etaMin_E && fabs(ElSCEta[i_preselected[i]]) < etaMax_E)) {
         hs_El_ElTransMass[2]->Fill(mt);
         hs_El_ElPt[2]->Fill(pt);
         hs_El_ElPhi[2]->Fill(phi);
         hs_El_ElEta[2]->Fill(eta);
         hs_El_ElMet[2]->Fill(met);
         hs_El_ElCharge[2]->Fill(charge);
         hs_El_Njets[2]->Fill(nJets);
         El_count[4]++;
       }
       
       // make the N-1 ID & isolation plots (VBTF WP)
       if(fabs(ElSCEta[i_preselected[i]])<etaMax_B) {
         if(WPElectronID(i_preselected[i], see,        WorkingPointCuts)==0) hn1_El_ElSigmaIetaIeta[0]->Fill(See);
         if(WPElectronID(i_preselected[i], dphi,       WorkingPointCuts)==0) hn1_El_ElDeltaPhi[0]->Fill(delPhi);
         if(WPElectronID(i_preselected[i], deta,       WorkingPointCuts)==0) hn1_El_ElDeltaEta[0]->Fill(delEta);
         if(WPElectronID(i_preselected[i], hoe,        WorkingPointCuts)==0) hn1_El_ElHoverE[0]->Fill(Hoe);
         if(WPElectronID(i_preselected[i], combiso,    WorkingPointCuts)==0) hn1_El_ElCombIso[0]->Fill(combIso_B);
         if(WPElectronID(i_preselected[i], conversion, WorkingPointCuts)==0) hn1_El_ElNumberOfMissingInnerHits[0]->Fill(missingHits);
       }
       if(fabs(ElSCEta[i_preselected[i]])>etaMin_E && fabs(ElSCEta[i_preselected[i]])<etaMax_E) {
         if(WPElectronID(i_preselected[i], see,        WorkingPointCuts)==0) hn1_El_ElSigmaIetaIeta[1]->Fill(See);
         if(WPElectronID(i_preselected[i], dphi,       WorkingPointCuts)==0) hn1_El_ElDeltaPhi[1]->Fill(delPhi);
         if(WPElectronID(i_preselected[i], deta,       WorkingPointCuts)==0) hn1_El_ElDeltaEta[1]->Fill(delEta);
         if(WPElectronID(i_preselected[i], hoe,        WorkingPointCuts)==0) hn1_El_ElHoverE[1]->Fill(Hoe);
         if(WPElectronID(i_preselected[i], combiso,    WorkingPointCuts)==0) hn1_El_ElCombIso[1]->Fill(combIso_E);
         if(WPElectronID(i_preselected[i], conversion, WorkingPointCuts)==0) hn1_El_ElNumberOfMissingInnerHits[1]->Fill(missingHits);
       }

       // look at IDed and isolated Electrons
       // barrel
       if(fabs(ElSCEta[i_preselected[i]]) < etaMax_B && 
          WPElectronID(i_preselected[i], all, WorkingPointCuts)==0) {
         hs_El_ElTransMass[3]->Fill(mt);
         hs_El_ElPt[3]->Fill(pt);
         hs_El_ElPhi[3]->Fill(phi);
         hs_El_ElEta[3]->Fill(eta);
         hs_El_ElMet[3]->Fill(met);
         hs_El_ElCharge[3]->Fill(charge);
         hs_El_Njets[3]->Fill(nJets);
         El_count[5]++;
       }
       // endcap
       if(fabs(ElSCEta[i_preselected[i]]) > etaMin_E && fabs(ElSCEta[i_preselected[i]]) < etaMax_E &&
          WPElectronID(i_preselected[i], all, WorkingPointCuts)==0) {
         hs_El_ElTransMass[1]->Fill(mt);
         hs_El_ElPt[4]->Fill(pt);
         hs_El_ElPhi[4]->Fill(phi);
         hs_El_ElEta[4]->Fill(eta);
         hs_El_ElMet[4]->Fill(met);
         hs_El_ElCharge[4]->Fill(charge);
         hs_El_Njets[4]->Fill(nJets);
         El_count[6]++;
       }
       // barrel + endcap
       if( (fabs(ElSCEta[i_preselected[i]]) < etaMax_B || (fabs(ElSCEta[i_preselected[i]]) > etaMin_E && fabs(ElSCEta[i_preselected[i]]) < etaMax_E) ) &&
          WPElectronID(i_preselected[i], all, WorkingPointCuts)==0) {
         hs_El_ElTransMass[5]->Fill(mt);
         hs_El_ElPt[5]->Fill(pt);
         hs_El_ElPhi[5]->Fill(phi);
         hs_El_ElEta[5]->Fill(eta);
         hs_El_ElMet[5]->Fill(met);
         hs_El_ElCharge[5]->Fill(charge);
         hs_El_Njets[5]->Fill(nJets);
         El_count[7]++;
       }
       
       // fill also the ETH ID&Iso N-1 plots at some point...
       
     }
     

     //--------------------------------------------------------------------------------------------------------------------
     // count the good leptons - pt > 10 / 25 GeV
     //--------------------------------------------------------------------------------------------------------------------
     int nLeptons10_B=0;
     int nLeptons25_B=0;
     for(int i=0; i<i_goodBarrelElectrons.size(); i++) {
       if(ElPt[i_goodBarrelElectrons[i]] > 10) nLeptons10_B++;
       if(ElPt[i_goodBarrelElectrons[i]] > 25) nLeptons25_B++;
     }
     int nLeptons10_E=0;
     int nLeptons25_E=0;
     for(int i=0; i<i_goodEndcapElectrons.size(); i++) {
       if(ElPt[i_goodEndcapElectrons[i]] > 10) nLeptons10_E++;
       if(ElPt[i_goodEndcapElectrons[i]] > 25) nLeptons25_E++;
     }
     //--------------------------------------------------------------------------------------------------------------------


//      cout<<"# vertices: "<<NVrtx<<endl;
//      for(int jvert=0; jvert<25; jvert++) {
//        if(VrtxTrkIndex[jvert][0] != -1) cout<<jvert<<". vertex: ----------------"<<endl;


//        double dphi = 0;
//        for(int jtrk=0; jtrk<25; jtrk++) {
//       if(VrtxTrkIndex[jvert][jtrk] != -1) cout<<VrtxTrkIndex[jvert][jtrk]<<", pt= "<< TrkPt[VrtxTrkIndex[jvert][jtrk]] <<", phi= "<<TrkPhi[VrtxTrkIndex[jvert][jtrk]]<<endl;
//        }
//      }


//      cout<<"here2"<<endl;

     //--------------------------------------------------------------------------------------------------------------------
     // Fill the Tree
     //--------------------------------------------------------------------------------------------------------------------

//      cout<<i_goodElectrons.size()<<endl;

     if(i_goodElectrons.size() < 1) continue;

//      cout<<"here2a"<<endl;
     

     // general info
     tRun         = Run;
     tEvent       = Event;
     tLumiSection = LumiSection;       
     
     // get the inst. lumi.
//      double avgInstLumi;
//      avgInstLumi = a.getDeliveredInstLumi(Run, LumiSection);
//      cout<<"avg. inst. lumi.: "<<avgInstLumi<<endl;


     // Electrons
     tElN = i_goodElectrons.size();
     for(int i=0; i<i_goodElectrons.size(); i++) {
     
//        if(ElIDsimpleWP80relIso[i_goodElectrons[i]] != 7) cout<<ElIDsimpleWP80relIso[i_goodElectrons[i]]<<"/"<<ElPt[i_goodElectrons[i]]<<endl;

       tElPt[i]     = ElPt[i_goodElectrons[i]];
       tElPx[i]     = ElPx[i_goodElectrons[i]];
       tElPy[i]     = ElPy[i_goodElectrons[i]];
       tElPz[i]     = ElPz[i_goodElectrons[i]];
       tElE[i]      = ElE[i_goodElectrons[i]];
       tElEta[i]    = ElEta[i_goodElectrons[i]];
       tElPhi[i]    = ElPhi[i_goodElectrons[i]];
       tElCharge[i] = ElCharge[i_goodElectrons[i]];
       tElDz[i]     = ElDzBS[i_goodElectrons[i]];
       tElPfMt[i]   = transverseMass(ElPt[i_goodElectrons[i]], PFMET,  deltaPhi(ElPhi[i_goodElectrons[i]], PFMETphi));
       tElCaloMt[i] = transverseMass(ElPt[i_goodElectrons[i]], RawMET, deltaPhi(ElPhi[i_goodElectrons[i]], RawMETphi));
       tElTcMt[i]   = transverseMass(ElPt[i_goodElectrons[i]], TCMET,  deltaPhi(ElPhi[i_goodElectrons[i]], TCMETphi));
     }

//      cout<<"here3"<<endl;
     

     // different Met
     tPfMet   = PFMET;
     tCaloMet = RawMET;
     tTcMet   = TCMET;

     // HF variables

     double SumEHF_minus  = 0.;
     double SumEHF_plus   = 0.;
     double SumEHF_Aminus  = 0.;
     double SumEHF_Aplus   = 0.;
     double SumEHF_Bminus  = 0.;
     double SumEHF_Bplus   = 0.;
     double SumETHF_minus = 0.;
     double SumETHF_plus  = 0.;
     int Ntowers_plus  = 0.;
     int Ntowers_minus = 0.;

     double SumE    = 0.;
     double minusPz = 0.;
     double plusPz  = 0.;

     double CaloTowersEMin_EB = 1.5;
     double CaloTowersEMin_EE = 2.5; // cut on ET !!!!!!!!!!!!!!!!!!
     double CaloTowersEMin_HB = 1.5;
     double CaloTowersEMin_HE = 2.;
     double CaloTowersEMin_HF = 4.;


//      cout<<"here4"<<endl;


     tCaloTowersN = NCaloTowers;
     for(int j=0; j<NCaloTowers; j++) {

       // get all calo towers
       tCaloTowersE[j] = CaloTowersE[j];
       tCaloTowersEt[j] = CaloTowersE[j] * sin(CaloTowersTheta[j]);
       tCaloTowersEta[j] = CaloTowersEta[j];
       tCaloTowersEmE[j] = CaloTowersEmE[j];
       tCaloTowersHadE[j] = CaloTowersHadE[j];
       tCaloTowersHasEB[j] = CaloTowersHasEB[j];
       tCaloTowersHasEE[j] = CaloTowersHasEE[j];
       tCaloTowersHasHB[j] = CaloTowersHasHB[j];
       tCaloTowersHasHE[j] = CaloTowersHasHE[j];
       tCaloTowersHasHF[j] = CaloTowersHasHF[j];
       tCaloTowersId[j] = CaloTowersId[j];

       // is the tower on the plus or minus side?
       int zside = (CaloTowersId[j]&0x2000)?(1):(-1);

//        int zside = 0;
//        if(CaloTowersEta[j]>0) zside = 1;
//        if(CaloTowersEta[j]<0) zside = -1;


       if(CaloTowersHasHF[j]==1 /*&& !CaloTowersHasHE[j]*/) {
         if(CaloTowersE[j] >= CaloTowersEMin_HF) {
           if(zside>0) {
             SumEHF_plus += CaloTowersE[j];
             SumETHF_plus += CaloTowersE[j] * sin(CaloTowersTheta[j]);
             Ntowers_plus++;

             if(fabs(CaloTowersEta[j]) >= 3 && fabs(CaloTowersEta[j]) <=4) {
               SumEHF_Aplus += CaloTowersE[j];
             }
             if(fabs(CaloTowersEta[j]) >= 4 && fabs(CaloTowersEta[j]) <=5) {
               SumEHF_Bplus += CaloTowersE[j];
             }

           }
           if(zside<0) {
             SumEHF_minus += CaloTowersE[j];
             SumETHF_minus += CaloTowersE[j] * sin(CaloTowersTheta[j]);
             Ntowers_minus++;

             if(fabs(CaloTowersEta[j]) >= 3 && fabs(CaloTowersEta[j]) <=4) {
               SumEHF_Aminus += CaloTowersE[j];
             }
             if(fabs(CaloTowersEta[j]) >= 4 && fabs(CaloTowersEta[j]) <=5) {
               SumEHF_Bminus += CaloTowersE[j];
             }

           }
         }
       }

       // calculate E \pm p_{z}
       if( (CaloTowersHasEB[j]==1 && CaloTowersEmE[j] > CaloTowersEMin_EB )  ||
           (CaloTowersHasEE[j]==1 && CaloTowersEmE[j] /* * sin(CaloTowersTheta[j])*/ > CaloTowersEMin_EE ) ||  // cut on ET, not on energy
           (CaloTowersHasHB[j]==1 && CaloTowersHadE[j] > CaloTowersEMin_HB ) ||
           (CaloTowersHasHE[j]==1 && CaloTowersHadE[j] > CaloTowersEMin_HE ) ||
           (CaloTowersHasHF[j]==1 && CaloTowersHadE[j] > CaloTowersEMin_HF)  ) {

         tSumE  += CaloTowersE[j];
         tSumPz += CaloTowersE[j] * cos(CaloTowersTheta[j]);

         SumE += CaloTowersE[j];
         if(zside > 0) plusPz  += CaloTowersE[j] * cos(CaloTowersTheta[j]);
         if(zside < 0) minusPz += CaloTowersE[j] * cos(CaloTowersTheta[j]);
         
       }
         
     }

     tSumEHF_minus  = SumEHF_minus;
     tSumEHF_plus   = SumEHF_plus;
     tSumEHF_Aminus  = SumEHF_Aminus;
     tSumEHF_Aplus   = SumEHF_Aplus;
     tSumEHF_Bminus  = SumEHF_Bminus;
     tSumEHF_Bplus   = SumEHF_Bplus;
     tSumETHF_minus = SumETHF_minus;
     tSumETHF_plus  = SumETHF_plus;
     tNtowers_plus  = Ntowers_plus;
     tNtowers_minus = Ntowers_minus;

     tEminusPz = SumE + minusPz;
     tEplusPz  = SumE + plusPz;


//      double E = SumEHF_plus + SumEHE_plus + SumEHB_plus + SumEEE_plus + SumEEB_plus + SumEHF_minus + SumEHE_minus + SumEHB_minus + SumEEE_minus + SumEEB_minus;
//      double minusPz = SumPzHF_minus + SumPzHE_minus + SumPzHB_minus + SumPzEE_minus + SumPzEB_minus;
//      double plusPz  = SumPzHF_plus  + SumPzHE_plus  + SumPzHB_plus  + SumPzEE_plus  + SumPzEB_plus;
//      tEminusPz = E + minusPz;
//      tEplusPz  = E + plusPz;

     // jets
     tJetN   = i_jets.size();
     for(int j=0; j<i_jets.size(); j++) {
       tJetPt[j]  = PFJPt[i_jets[j]];
       tJetEta[j] = PFJEta[i_jets[j]];
       tJetPhi[j] = PFJPhi[i_jets[j]];
     }

     // \eta_{max} and \eta_{min}
     double EtaMax = -9.99;
     double EtaMin =  9.99;
     for(int j=0; j < NPf; j++) {
       if(PfEta[j] > EtaMax) EtaMax = PfEta[j];
       if(PfEta[j] < EtaMin) EtaMin = PfEta[j];
       // keep more PFlow variables
       tPfE[j]   = PfE[j];
       tPfEt[j]  = PfEt[j];
       tPfEta[j] = PfEta[j];
       tPfPt[j]  = PfPt[j];
     }
     tEtaMax=EtaMax;
     tEtaMin=EtaMin;
     tPfN = NPf;


     // count the vertices
//      int nVert = 0;
//      for(int j=0; j<NVrtx; j++) {
//        if(VrtxNdof[j] > 4) nVert++;

//      }
     tVertexN = NVrtx;
     for(int j=0; j<NVrtx; j++) {
       tVertexNdof[j] = VrtxNdof[j];
       tVertexChi2[j] = VrtxChi2[j];
       tVertexX[j] = VrtxX[j];
       tVertexY[j] = VrtxY[j];
       tVertexZ[j] = VrtxZ[j];
       //hhh
       for(int jj=0; jj<400; jj++) {
         tVertexTrackIndex[j][jj] = VrtxTrkIndex[j][jj];
//       cout<<tVertexTrackIndex[j][jj]<<endl;
//       tTVertexTrackWeight[j][jj] = VrtxTrkWeight[j][jj];
       }
     }
     

     // count the tracks
//      int nTrck = 0;
//      tTracksN = nTrck;


      cout<<"Filling pile up"<<endl;

     // pile-up
     tPUnumInteractions = PUnumInteractions;
     for(int z=0; z<PUnumInteractions; z++) {
       tPUzPositions[z] = PUzPositions[z];
     }
     cout<<"done!"<<endl;

     tBeamSpotX = Beamspotx;
     tBeamSpotY = Beamspoty;
     tBeamSpotZ = Beamspotz;

     tTracksN = NTracks;
     for(int tr=0; tr<NTracks; tr++) {
       tTrackPt[tr] = TrkPt[tr];
       tTrackPx[tr] = TrkPx[tr];
       tTrackPy[tr] = TrkPy[tr];
       tTrackPz[tr] = TrkPz[tr];
       tTrackDzBs[tr] = Trkdzbs[tr];
       tTrackEta[tr] = TrkEta[tr];
       tTrackPhi[tr] = TrkPhi[tr];
//        tTracktdsz[t] = Trkdsz[t];
//        tTrackdxy[t] = Trkdxy[t];

     }


     //--------------------------------------------------------------------------------------------------------------------


     // fill the histograms for the barrel only
     for(int i=0; i<i_goodBarrelElectrons.size(); i++) {

       // calculate the variables to fill
       double deltaPhiElMet=deltaPhi(ElPhi[i_goodBarrelElectrons[i]], PFMETphi);
       double mt = transverseMass(ElPt[i_goodBarrelElectrons[i]] * PtCorrection, met, deltaPhiElMet);
       int charge = ElCharge[i_goodBarrelElectrons[i]];
       double pt = ElPt[i_goodBarrelElectrons[i]] * PtCorrection;
       double phi = ElPhi[i_goodBarrelElectrons[i]];
       double eta = ElEta[i_goodBarrelElectrons[i]];


       // all barrel electrons after ID&Iso with pt > 10 GeV
       hs_W_ElTransMass[0]->Fill(mt);
       hs_W_ElPt[0]->Fill(pt);
       hs_W_ElPhi[0]->Fill(phi);
       hs_W_ElEta[0]->Fill(eta);
       hs_W_ElMet[0]->Fill(met);
       hs_W_ElCharge[0]->Fill(charge);
       hs_W_Njets[0]->Fill(nJets);
       
       // cut on the pt of the electrons
       if(ElPt[i_goodBarrelElectrons[i]] > 25) {

         // all barrel electrons after ID&Iso with pt > 25 GeV
         debug4++;
         hs_W_ElTransMass[1]->Fill(mt);
         hs_W_ElPt[1]->Fill(pt);
         hs_W_ElPhi[1]->Fill(phi);
         hs_W_ElEta[1]->Fill(eta);
         hs_W_ElMet[1]->Fill(met);
         hs_W_ElCharge[1]->Fill(charge);
         hs_W_Njets[1]->Fill(nJets);
         
         //--------------------------------------------------------------------------------------------------------------------
         //--------------------------------------------------------------------------------------------------------------------
         // THIS IS THE W BOSON ANALYSIS */*
         //--------------------------------------------------------------------------------------------------------------------
         //--------------------------------------------------------------------------------------------------------------------
         
         
         // exactly one lepton
         if(nLeptons25_B == 1) {
           hs_W_ElTransMass[2]->Fill(mt);
           hs_W_ElPt[2]->Fill(pt);
           hs_W_ElPhi[2]->Fill(phi);
           hs_W_ElEta[2]->Fill(eta);
           hs_W_ElMet[2]->Fill(met);
           hs_W_ElCharge[2]->Fill(charge);
           hs_W_Njets[2]->Fill(nJets);

           // MET cut
           if(met > 30) {
             hs_W_ElTransMass[3]->Fill(mt);
             hs_W_ElPt[3]->Fill(pt);
             hs_W_ElPhi[3]->Fill(phi);
             hs_W_ElEta[3]->Fill(eta);
             hs_W_ElMet[3]->Fill(met);
             hs_W_ElCharge[3]->Fill(charge);
             hs_W_Njets[3]->Fill(nJets);

             // transverse mass cut
             if(mt > 60) {
               hs_W_ElTransMass[4]->Fill(mt);
               hs_W_ElPt[4]->Fill(pt);
               hs_W_ElPhi[4]->Fill(phi);
               hs_W_ElEta[4]->Fill(eta);
               hs_W_ElMet[4]->Fill(met);
               hs_W_ElCharge[4]->Fill(charge);
               hs_W_Njets[4]->Fill(nJets);
               
               // zero jet bin
               if(nJets==0) {
                 hs_W_ElTransMass[5]->Fill(mt);
                 hs_W_ElPt[5]->Fill(pt);
                 hs_W_ElPhi[5]->Fill(phi);
                 hs_W_ElEta[5]->Fill(eta);
                 hs_W_ElMet[5]->Fill(met);
                 hs_W_ElCharge[5]->Fill(charge);
                 hs_W_Njets[5]->Fill(nJets);     
                 

//               // DIFFRACTIVE
//               if(isDiffractive) {
                   
//                 diff_run.push_back(Run);
//                 diff_event.push_back(Event);
//                 diff_ls.push_back(LumiSection);
//                 diff_sumehf_plus.push_back(SumEHF_plus);
//                 diff_sumehf_minus.push_back(SumEHF_minus);
//                 diff_sumethf_plus.push_back(SumETHF_plus);
//                 diff_sumethf_minus.push_back(SumETHF_minus);
//                 diff_ntowershf_plus.push_back(NTowersHF_plus);
//                 diff_ntowershf_minus.push_back(NTowersHF_minus);
//               }
                 

               // W minus
               if(charge < 0) {
                 hs_W_ElTransMass[6]->Fill(mt);
                 hs_W_ElPt[6]->Fill(pt);
                 hs_W_ElPhi[6]->Fill(phi);
                 hs_W_ElEta[6]->Fill(eta);
                 hs_W_ElMet[6]->Fill(met);
                 hs_W_ElCharge[6]->Fill(charge);
                 hs_W_Njets[6]->Fill(nJets);
               }
               // W plus
               else if(charge > 0) {
                 hs_W_ElTransMass[7]->Fill(mt);
                 hs_W_ElPt[7]->Fill(pt);
                 hs_W_ElPhi[7]->Fill(phi);
                 hs_W_ElEta[7]->Fill(eta);
                 hs_W_ElMet[7]->Fill(met);
                 hs_W_ElCharge[7]->Fill(charge);
                 hs_W_Njets[7]->Fill(nJets);
               }
               }
               // 1 jet bin
               if(nJets==1) {
                 // get the angle between the jet and the elctron
                 double deltPhi = deltaPhi(phi, JPhi[i_jets[0]]);
//               cout<<"delta phi: "<<deltPhi<<endl;
                 hs_W_ElDeltaPhi[8]->Fill(deltPhi);

                 hs_W_ElTransMass[8]->Fill(mt);
                 hs_W_ElPt[8]->Fill(pt);
                 hs_W_ElPhi[8]->Fill(phi);
                 hs_W_ElEta[8]->Fill(eta);
                 hs_W_ElMet[8]->Fill(met);
                 hs_W_ElCharge[8]->Fill(charge);
                 hs_W_Njets[8]->Fill(nJets);
                 
                 if(deltPhi < 2.793) {
                   hs_W_ElTransMass[9]->Fill(mt);
                   hs_W_ElPt[9]->Fill(pt);
                   hs_W_ElPhi[9]->Fill(phi);
                   hs_W_ElEta[9]->Fill(eta);
                   hs_W_ElMet[9]->Fill(met);
                   hs_W_ElCharge[9]->Fill(charge);
                   hs_W_Njets[9]->Fill(nJets);           
                   
                   
                   // get the W pt spectrum
                   hs_W_ElDeltaPhi[9]->Fill(deltPhi);
                   hs_W_jetPt[9]->Fill(JPt[i_jets[0]]);

                   // W minus
                   if(charge < 0) {
                     hs_W_ElTransMass[10]->Fill(mt);
                     hs_W_ElPt[10]->Fill(pt);
                     hs_W_ElPhi[10]->Fill(phi);
                     hs_W_ElEta[10]->Fill(eta);
                     hs_W_ElMet[10]->Fill(met);
                     hs_W_ElCharge[10]->Fill(charge);
                     hs_W_Njets[10]->Fill(nJets);                
                     
                     
                   }
                   // W plus
                   else if(charge > 0) {
                     hs_W_ElTransMass[11]->Fill(mt);
                     hs_W_ElPt[11]->Fill(pt);
                     hs_W_ElPhi[11]->Fill(phi);
                     hs_W_ElEta[11]->Fill(eta);
                     hs_W_ElMet[11]->Fill(met);
                     hs_W_ElCharge[11]->Fill(charge);
                     hs_W_Njets[11]->Fill(nJets);                
                   }
                 }
               }
               // 2 jet bin
               if(nJets>=2) {
                 hs_W_ElTransMass[12]->Fill(mt);
                 hs_W_ElPt[12]->Fill(pt);
                 hs_W_ElPhi[12]->Fill(phi);
                 hs_W_ElEta[12]->Fill(eta);
                 hs_W_ElMet[12]->Fill(met);
                 hs_W_ElCharge[12]->Fill(charge);
                 hs_W_Njets[12]->Fill(nJets);
                 
                 // W minus
                 if(charge < 0) {
                   hs_W_ElTransMass[13]->Fill(mt);
                   hs_W_ElPt[13]->Fill(pt);
                   hs_W_ElPhi[13]->Fill(phi);
                   hs_W_ElEta[13]->Fill(eta);
                   hs_W_ElMet[13]->Fill(met);
                   hs_W_ElCharge[13]->Fill(charge);
                   hs_W_Njets[13]->Fill(nJets);          
                   
                   
                 }
                 // W plus
                 else if(charge > 0) {
                   hs_W_ElTransMass[14]->Fill(mt);
                   hs_W_ElPt[14]->Fill(pt);
                   hs_W_ElPhi[14]->Fill(phi);
                   hs_W_ElEta[14]->Fill(eta);
                   hs_W_ElMet[14]->Fill(met);
                   hs_W_ElCharge[14]->Fill(charge);
                   hs_W_Njets[14]->Fill(nJets);          
                 }
               }


             }
           } // MET cut
           
           // make N-1 plots for 0/1 jet bin (already selected == 1 electron)
           if(nJets==0) {
             hs_W_ElTransMass[15]->Fill(mt);
             hs_W_ElPt[15]->Fill(pt);
             hs_W_ElPhi[15]->Fill(phi);
             hs_W_ElEta[15]->Fill(eta);
             hs_W_ElMet[15]->Fill(met);
             hs_W_ElCharge[15]->Fill(charge);
             hs_W_Njets[15]->Fill(nJets);                
             
             if(met > 30) {
               hs_W_ElTransMass[16]->Fill(mt);
               hs_W_ElPt[16]->Fill(pt);
               hs_W_ElPhi[16]->Fill(phi);
               hs_W_ElEta[16]->Fill(eta);
               hs_W_ElMet[16]->Fill(met);
               hs_W_ElCharge[16]->Fill(charge);
               hs_W_Njets[16]->Fill(nJets);
             }
             if(mt > 60) {
               hs_W_ElTransMass[17]->Fill(mt);
               hs_W_ElPt[17]->Fill(pt);
               hs_W_ElPhi[17]->Fill(phi);
               hs_W_ElEta[17]->Fill(eta);
               hs_W_ElMet[17]->Fill(met);
               hs_W_ElCharge[17]->Fill(charge);
               hs_W_Njets[17]->Fill(nJets);
             }
           }
           if(nJets==1) {
             hs_W_ElTransMass[18]->Fill(mt);
             hs_W_ElPt[18]->Fill(pt);
             hs_W_ElPhi[18]->Fill(phi);
             hs_W_ElEta[18]->Fill(eta);
             hs_W_ElMet[18]->Fill(met);
             hs_W_ElCharge[18]->Fill(charge);
             hs_W_Njets[18]->Fill(nJets);

             if(met > 30) {
             hs_W_ElTransMass[19]->Fill(mt);
             hs_W_ElPt[19]->Fill(pt);
             hs_W_ElPhi[19]->Fill(phi);
             hs_W_ElEta[19]->Fill(eta);
             hs_W_ElMet[19]->Fill(met);
             hs_W_ElCharge[19]->Fill(charge);
             hs_W_Njets[19]->Fill(nJets);
             }
             if(mt > 60) {
             hs_W_ElTransMass[20]->Fill(mt);
             hs_W_ElPt[20]->Fill(pt);
             hs_W_ElPhi[20]->Fill(phi);
             hs_W_ElEta[20]->Fill(eta);
             hs_W_ElMet[20]->Fill(met);
             hs_W_ElCharge[20]->Fill(charge);
             hs_W_Njets[20]->Fill(nJets);
             }
           }
         }  // 1 electron
       } // pt > 25 GeV

     } // loop over good barrel electrons


       //--------------------------------------------------------------------------------------------------------------------
       //--------------------------------------------------------------------------------------------------------------------
       // THIS IS THE Z BOSON ANALYSIS */*
       //--------------------------------------------------------------------------------------------------------------------
       //--------------------------------------------------------------------------------------------------------------------
       
       // 2 electrons in the barrel
     if(nLeptons25_B >= 2) {

         // calculate the invariant mass with the 2 highest pt leptons
         double invMass;
         TLorentzVector e1, e2;
         int index1=0;
         int index2=0;
         double pt1=0;
         double pt2=0;
         for(int i=0; i<i_goodBarrelElectrons.size(); i++) {
           if(ElPt[i_goodBarrelElectrons[i]] > pt2 && ElPt[i_goodBarrelElectrons[i]] > pt1) {index2 = index1; pt2 = pt1; index1 = i_goodBarrelElectrons[i]; pt1 = ElPt[i_goodBarrelElectrons[i]];}
           if(ElPt[i_goodBarrelElectrons[i]] > pt2 && ElPt[i_goodBarrelElectrons[i]] < pt1) {index2 =  i_goodBarrelElectrons[i]; pt2 = ElPt[i_goodBarrelElectrons[i]];}
         }
  
         e1.SetPtEtaPhiE(ElPt[index1] * PtCorrection ,ElEta[index1],ElPhi[index1],ElE[index1]);
         e2.SetPtEtaPhiE(ElPt[index2] * PtCorrection ,ElEta[index2],ElPhi[index2],ElE[index2]);
         invMass = (e1+e2).M();  

         // Z mass cut flag
         bool Zmass=false;
         if(invMass > 83.7 && invMass < 98.7) Zmass=true;


         // these variables are only defined if there really are jets...

         double leadingJetPt = 0.;
         double leadingJetEta= 0.;
         double leadingJetPhi= 0.;
         double summedJetPt  = 0.;
         double summedJetEta = 0.;
         double summedJetPhi = 0.;
         double recoilPt     = 0.;
         double recoilEta    = 0.;
         double recoilPhi    = 0.;

         TLorentzVector jsum;
         for(int ji=0; ji<i_jets.size(); ji++) {
           TLorentzVector vec;
           vec.SetPxPyPzE(JPx[i_jets[ji]], JPy[i_jets[ji]], JPz[i_jets[ji]], JE[i_jets[ji]]);
           //find the leading jet
           if(vec.Pt() > leadingJetPt) {
             leadingJetPt  = vec.Pt();
             leadingJetEta = vec.Eta();
             leadingJetPhi = vec.Phi();
           }
           // sum all jets
           jsum=jsum+vec;
         }
         summedJetPt = jsum.Pt();
         summedJetEta = jsum.Eta();
         summedJetPhi = jsum.Phi();
         // add also the Z
         jsum+=e1;
         jsum+=e2;
         recoilPt = jsum.Pt();
         recoilEta = jsum.Eta();
         recoilPhi = jsum.Phi();
 
         h_Z_invMass[0]->Fill(invMass);
         h_Z_invMass2[0]->Fill(invMass);
         h_Z_ptAll[0]->Fill(ElPt[index1]);
         h_Z_ptAll[0]->Fill(ElPt[index2]);
         h_Z_etaAll[0]->Fill(ElEta[index1]);
         h_Z_etaAll[0]->Fill(ElEta[index2]);
         h_Z_phiAll[0]->Fill(ElPhi[index1]);
         h_Z_phiAll[0]->Fill(ElPhi[index2]);
         h_Z_deltaPhi[0]->Fill(deltaPhi(ElPhi[index1], ElPhi[index2]));
         if(Zmass) {
           h_Z_ptDiLepton[0]->Fill((e1+e2).Pt());
           h_Z_ptDiLepton2[0]->Fill((e1+e2).Pt());
           h_Z_etaDiLepton[0]->Fill((e1+e2).Eta());
           h_Z_phiDiLepton[0]->Fill((e1+e2).Phi());
           h_Z_metPt[0]->Fill(met);
           h_Z_metPhi[0]->Fill(metPhi);

//         // diffractive Zs
//         if(isDiffractive) {
//           diff_run.push_back(Run);
//           diff_event.push_back(Event);
//           diff_ls.push_back(LumiSection);
//           diff_sumehf_plus.push_back(SumEHF_plus);
//           diff_sumehf_minus.push_back(SumEHF_minus);
//           diff_sumethf_plus.push_back(SumETHF_plus);
//           diff_sumethf_minus.push_back(SumETHF_minus);
//           diff_ntowershf_plus.push_back(NTowersHF_plus);
//           diff_ntowershf_minus.push_back(NTowersHF_minus);
//         }
           
           if(nJets == 0) {
             sel_ZOjet_run.push_back(Run);
             sel_ZOjet_event.push_back(Event);
             sel_ZOjet_ls.push_back(LumiSection);
             
           }


         }

         // now variables with jet activity
         
         h_Z_nJets[0]->Fill(nJets);
         if(nJets > 0 && Zmass) {
           h_Z_leadingJetPt[0]->Fill(leadingJetPt);
           h_Z_leadingJetEta[0]->Fill(leadingJetEta);
           h_Z_leadingJetPhi[0]->Fill(leadingJetPhi);
           h_Z_summedJetPt[0]->Fill(summedJetPt);
           h_Z_summedJetEta[0]->Fill(summedJetEta);
           h_Z_summedJetPhi[0]->Fill(summedJetPhi);
           
           h_Z_leadingJetDeltaPt[0]->Fill(leadingJetPt - (e1+e2).Pt());
           h_Z_leadingJetDeltaEta[0]->Fill(leadingJetEta - (e1+e2).Eta());
           h_Z_leadingJetDeltaPhi[0]->Fill(deltaPhi(leadingJetPhi, (e1+e2).Phi()));
           h_Z_summedJetDeltaPt[0]->Fill(summedJetPt - (e1+e2).Pt());
           h_Z_summedJetDeltaEta[0]->Fill(summedJetEta - (e1+e2).Eta());
           h_Z_summedJetDeltaPhi[0]->Fill(deltaPhi(summedJetPhi, (e1+e2).Phi()));
           
           h_Z_recoilPt[0]->Fill(recoilPt);
           h_Z_recoilEta[0]->Fill(recoilEta);
           h_Z_recoilPhi[0]->Fill(recoilPhi);
         }
         
         
     } // 2 electrons
     
       
     // ******************************************************************************
     // ******************************************************************************


//      for(int i=0; i< tElN; i++) {
//        cout<<tElPt[i]<<endl;

//      }


     subTree->Fill();


       ////////////////////////////////////////////////////////////////////
       ////////////////////// User's code ends here ///////////////////////
       ////////////////////////////////////////////////////////////////////
       
    } // End loop over events


 END:
   cerr<<"Program ended after "<<sampleSize<<" events ..."<<endl;
   

   for(int b=0;b<MAXWcounter;b++) {
     hc_W_counter->SetBinContent(b+1,Wcounter[b]);
   }
   for(int b=0;b<counterSize;b++) {
     hc_El_counter->SetBinContent(b+1,counter[b]);
   }
   

   fclose(out);
   fclose(out_Z);

   // Summary
   cout<<"*****************************************"<<endl;
   cout<<"* SUMMARY                               *"<<endl;
   cout<<"*****************************************"<<endl;
   cout<<"# Events:            "<<El_count[0]<<endl;
   cout<<"HLT passed:          "<<El_count[1]<<endl;
   cout<<"*****************************************"<<endl;
   cout<<"Electrons"<<endl;
   cout<<"*****************************************"<<endl;
   cout<<"Preselected e in EB: "<<El_count[2]<<endl;
   cout<<"Preselected e in EE: "<<El_count[3]<<endl;
   cout<<"Presele. e in EE+EB: "<<El_count[4]<<endl;
   cout<<"*****************************************"<<endl;
   cout<<"IDed & Isol. in EB:  "<<El_count[5]<<endl;
   cout<<"IDed & Isol. in EE:  "<<El_count[6]<<endl;
   cout<<"IDed & Iso. in EB+EE:"<<El_count[7]<<endl;
   cout<<"*****************************************"<<endl;
   cout<<"Ws"<<endl;


   cout<<"*****************************************"<<endl;


   for(int ci=0; ci<15; ci++) {
     cout<<ci<<": "<<counter[ci]<<endl;
   }

   if(diff_run.size()>0) {
     cout<<"DIFFRACTION STUFF!"<<endl;
     cout<<"Run / Event / LS / E HF+ / E HF- / ET HF+ / ET hf- / #tow HF+ / #tow HF- "<<endl;
     for(int i=0; i<diff_run.size(); i++) {
       cout<< diff_run[i] << " / " <<diff_event[i]<<" / "<<diff_ls[i] << " / "
           << diff_sumehf_plus[i]    << " / " << diff_sumehf_minus[i] << " / "
           << diff_sumethf_plus[i]   << " / " << diff_sumethf_minus[i] << " / "
           << diff_ntowershf_plus[i] << " / " << diff_ntowershf_minus[i]
           <<endl;
     } 
   }

   if(sel_ZOjet_run.size() > 0) cout<<"Zs in barrel & 0 jet bin:"<<endl;
   for(int i=0; i<sel_ZOjet_run.size(); i++) {
     cout<<sel_ZOjet_run[i]<<" / "<<sel_ZOjet_event[i]<<" / "<<sel_ZOjet_ls[i]<<endl;
   }

   delete [] counter;
   counter=NULL;
   delete [] Wcounter;
   Wcounter=NULL;


   subTree->Write();
   subTree->Print();


   //////////write histos 
   h_jetPt->Write();

   
   // sum barrel and endcap into general plots
   for(int a=34; a<MAXE; a++) {
     hs_El_ElPt[a]->Add(hs_El_ElPt[a-7]);
     hs_El_ElPt[a]->Add(hs_El_ElPt[a-14]);
     hs_El_ElEta[a]->Add(hs_El_ElEta[a-7]);
     hs_El_ElEta[a]->Add(hs_El_ElEta[a-14]);
     hs_El_ElPhi[a]->Add(hs_El_ElPhi[a-7]);
     hs_El_ElPhi[a]->Add(hs_El_ElPhi[a-14]);
     hs_El_ElCharge[a]->Add(hs_El_ElCharge[a-7]);
     hs_El_ElCharge[a]->Add(hs_El_ElCharge[a-14]);
     hs_El_ElDeltaPhi[a]->Add(hs_El_ElDeltaPhi[a-7]);
     hs_El_ElDeltaPhi[a]->Add(hs_El_ElDeltaPhi[a-14]);
     hs_El_ElDeltaEta[a]->Add(hs_El_ElDeltaEta[a-7]);
     hs_El_ElDeltaEta[a]->Add(hs_El_ElDeltaEta[a-14]);
     hs_El_ElTransMass[a]->Add(hs_El_ElTransMass[a-7]);
     hs_El_ElTransMass[a]->Add(hs_El_ElTransMass[a-14]);
     hs_El_ElHoverE[a]->Add(hs_El_ElHoverE[a-7]);
     hs_El_ElHoverE[a]->Add(hs_El_ElHoverE[a-14]);
     hs_El_ElSigmaIetaIeta[a]->Add(hs_El_ElSigmaIetaIeta[a-7]);
     hs_El_ElSigmaIetaIeta[a]->Add(hs_El_ElSigmaIetaIeta[a-14]);
     hs_El_ElNumberOfMissingInnerHits[a]->Add(hs_El_ElNumberOfMissingInnerHits[a-7]);
     hs_El_ElNumberOfMissingInnerHits[a]->Add(hs_El_ElNumberOfMissingInnerHits[a-14]);
     hs_El_ElCombIso[a]->Add(hs_El_ElCombIso[a-7]);
     hs_El_ElCombIso[a]->Add(hs_El_ElCombIso[a-14]);
   }


   // write the plots to the *.root file
   
   h_numberOfEvents->Write();
   h_HLTInfo->Write();


   if(writeIdPlots==true) {
     for(int i=0;i<MAXE;i++) {
       hs_El_ElPt[i]->Write();
       hs_El_ElEta[i]->Write();
       hs_El_ElPhi[i]->Write();
       hs_El_ElDeltaPhi[i]->Write();
       hs_El_ElMet[i]->Write();
       hs_El_ElTransMass[i]->Write();
       hs_El_ElCharge[i]->Write();
       hs_El_Njets[i]->Write();
       hs_El_ElDeltaPhi[i]->Write();
       hs_El_ElDeltaEta[i]->Write();
       hs_El_ElOneOverEminusOneOverP[i]->Write();
       hs_El_ElTransMass[i]->Write();
       hs_El_Njets[i]->Write();
       hs_El_ElEthIso[i]->Write();
       hs_El_ElTrkIso[i]->Write();
       hs_El_ElEcalIso[i]->Write();
       hs_El_ElHcalIso[i]->Write();
       hs_El_ElCombIso[i]->Write();
       hs_El_ElHoverE[i]->Write();
       hs_El_ElSigmaIetaIeta[i]->Write();
       hs_El_ElNumberOfMissingInnerHits[i]->Write();
     }
     hc_El_counter->Write();
   }


   if(writeN1Plots==true) {
     for(int i=0; i<MAXN; i++) {
       hn1_El_ElDeltaPhi[i]->Write();
       hn1_El_ElDeltaEta[i]->Write();
       hn1_El_ElOneOverEminusOneOverP[i]->Write();

//        hn1_El_ElPt[i]->Write();
//        hn1_El_ElEta[i]->Write();
//        hn1_El_ElPhi[i]->Write();
//        hn1_El_ElMet[i]->Write();
//        hn1_El_ElTransMass[i]->Write();
//        hn1_El_ElCharge[i]->Write();
//        hn1_El_Njets[i]->Write();
//        hn1_El_ElTransMass[i]->Write();
//        hn1_El_Njets[i]->Write();
       hn1_El_ElEthIso[i]->Write();
//        hn1_El_ElTrkIso[i]->Write();
//        hn1_El_ElEcalIso[i]->Write();
//        hn1_El_ElHcalIso[i]->Write();
       hn1_El_ElCombIso[i]->Write();
       hn1_El_ElHoverE[i]->Write();
       hn1_El_ElSigmaIetaIeta[i]->Write();
       hn1_El_ElNumberOfMissingInnerHits[i]->Write();
     }
   }
   if(writeWenuPlots) {
     for(int i=0;i<MAXS;i++) {
       hs_W_ElPt[i]->Write();
       hs_W_ElEta[i]->Write();
       hs_W_ElPhi[i]->Write();
       hs_W_ElDeltaPhi[i]->Write();
       hs_W_ElMet[i]->Write();
       hs_W_ElTransMass[i]->Write();
       hs_W_ElCharge[i]->Write();
       hs_W_Njets[i]->Write();
       hs_W_jetPt[i]->Write();
     }
     hc_W_counter->Write();
   }

   for(int i=0;i<MAXZ;i++) {
     h_Z_invMass[i]->Write();
     h_Z_invMass2[i]->Write();
     h_Z_ptAll[i]->Write();
     h_Z_etaAll[i]->Write();
     h_Z_phiAll[i]->Write();
     h_Z_deltaPhi[i]->Write();
     h_Z_ptDiLepton[i]->Write();
     h_Z_ptDiLepton2[i]->Write();
     h_Z_etaDiLepton[i]->Write();
     h_Z_phiDiLepton[i]->Write();

     h_Z_leadingJetPt[i]->Write();
     h_Z_summedJetPt[i]->Write();
     h_Z_leadingJetEta[i]->Write();
     h_Z_summedJetEta[i]->Write();
     h_Z_leadingJetPhi[i]->Write();
     h_Z_summedJetPhi[i]->Write();
     h_Z_nJets[i]->Write();
     h_Z_leadingJetDeltaPt[i]->Write();
     h_Z_leadingJetDeltaEta[i]->Write();
     h_Z_leadingJetDeltaPhi[i]->Write();
     h_Z_summedJetDeltaPt[i]->Write();
     h_Z_summedJetDeltaEta[i]->Write();
     h_Z_summedJetDeltaPhi[i]->Write();
     h_Z_recoilPt[i]->Write();
     h_Z_recoilEta[i]->Write();
     h_Z_recoilPhi[i]->Write();
     h_Z_metPt[i]->Write();
     h_Z_metPhi[i]->Write();

   }
     std::cout << "analysisClass::Loop() ends" <<std::endl;   
}


Revision:	1.6
Committed:	Tue Dec 21 15:25:39 2010 UTC (14 years, 4 months ago) by jueugste
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.5:	+205 -26 lines
Log Message:	Added generated PU info