WZAnalysis/src/WZAnalyzer.cc

// -*- C++ -*-
//
// Package:    WZAnalyzer
// Class:      WZAnalyzer
// 
/**\class WZAnalyzer WZAnalyzer.cc Vuko/WZAnalysis/src/WZAnalyzer.cc

 Description: <one line class summary>

 Implementation:
     <Notes on implementation>
*/
//
// Original Author:  Vuko Brigljevic
//         Created:  Fri Nov  9 11:07:27 CET 2007
// $Id: WZAnalyzer.cc,v 1.55 2008/04/22 22:14:59 kkaadze Exp $
//
//


// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "Vuko/WZAnalysis/interface/WZAnalyzer.h"
#include "DataFormats/Candidate/interface/OverlapChecker.h"

#include "Vuko/WZAnalysis/interface/ElectronProperties.h"
#include "Vuko/WZAnalysis/interface/MuonProperties.h"
#include "Vuko/WZAnalysis/interface/AnglesUtil.h"
#include "DataFormats/Common/interface/TriggerResults.h"

//--- muon AOD:
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/EgammaCandidates/interface/Electron.h"
#include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuIsoDeposit.h"
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
#include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectron.h"

#include "DataFormats/METReco/interface/CaloMET.h"


#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TTree.h" 
//#include "TVector3.h"

#include <map>

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
WZAnalyzer::WZAnalyzer(const edm::ParameterSet& iConfig)

{
  //now do what ever initialization is needed

  fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile"); 
  theMCTruthCollection         = iConfig.getParameter<edm::InputTag>("genParticles");

  //  theMCTruthCollection         = iConfig.getParameter<edm::InputTag>("genParticleCollection");

  theLooseMuonCollection       = iConfig.getParameter<edm::InputTag>("LooseMuons");
  theGoodMuonCollection        = iConfig.getParameter<edm::InputTag>("GoodMuons");
  theLooseElectronCollection   = iConfig.getParameter<edm::InputTag>("LooseElectrons");
  theAllElectronCollection    = iConfig.getParameter<edm::InputTag>("AllElectrons");
  theTightLeptonCollection     = iConfig.getParameter<edm::InputTag>("TightLeptons");
  //  theMediumElectronCollection  = iConfig.getParameter<edm::InputTag>("LooseElectrons");
  //  theTightElectronCollection   = iConfig.getParameter<edm::InputTag>("TightElectrons");
  // Z's
  theZeeCollection             = iConfig.getParameter<edm::InputTag>("ZtoEE");
  theZmumuCollection           = iConfig.getParameter<edm::InputTag>("ZtoMuMu");
  theZllCollection             = iConfig.getParameter<edm::InputTag>("ZtoLL");

  theJetCollection             = iConfig.getParameter<edm::InputTag>("Jets");

  theWeightCollection          = iConfig.getParameter<edm::InputTag>("EventWeight");
  alpgenIDTag                  =  iConfig.getParameter<edm::InputTag>("AlpgenIDTag");  

  getAlpgenID = false;
  getAlpgenID                  = iConfig.getParameter<bool>("getAlpgenID");

  storeHLTresults=false;
  storeHLTresults              = iConfig.getParameter<bool>("storeHLTresults");

  getEventWeight = false;
  getEventWeight               = iConfig.getParameter<bool>("getEventWeight");

  triggerResultsTag            = iConfig.getParameter<edm::InputTag>("TriggerResults");
  firstTriggerResult = true;

  findBackJets = false;
  findBackJets                 = iConfig.getParameter<bool>("findBackJets"); 
}


WZAnalyzer::~WZAnalyzer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
WZAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;
   using namespace reco;
   using namespace std;

   ///////////////////////////////////////
   //
   /// EVENT INITIALISATION
   ///

   //read run number and event ID number
   RunNumber = iEvent.id().run();
   EventID   = iEvent.id().event();

   // Reset a few things
   for (map<string,wzana::LeptonProperties* >::iterator iter = leptonProperties.begin();
        iter!=leptonProperties.end();   iter++) 
     {
       iter->second->ResetValues();
     }


   const Candidate * theZCandidate = 0;
   const Candidate * theWlepton    = 0;


   ////////////////////////////////////////
   // Get lists

   // Z->mumu
   Handle<CandidateCollection> zmumuCands;
   iEvent.getByLabel( theZmumuCollection , zmumuCands);


   // Z->ee
   Handle<CandidateCollection> zeeCands;
   iEvent.getByLabel( theZeeCollection , zeeCands);

   // Z->ee
   Handle<CandidateCollection> zllCands;
   iEvent.getByLabel( theZllCollection , zllCands);

   // loose electrons
   Handle<CandidateCollection> looseElectronCands;
   iEvent.getByLabel( theLooseElectronCollection , looseElectronCands);

   // All electrons
   Handle<CandidateCollection> allElectronCands;
   iEvent.getByLabel( theAllElectronCollection , allElectronCands);

   // loose muons
   Handle<CandidateCollection> looseMuonCands;
   iEvent.getByLabel( theLooseMuonCollection , looseMuonCands);

   // good muons
   Handle<CandidateCollection> goodMuonCands;
   iEvent.getByLabel( theGoodMuonCollection , goodMuonCands);


   // tight leptons
   Handle<CandidateCollection> tightLeptonCands;
   iEvent.getByLabel( theTightLeptonCollection , tightLeptonCands);

   // jets

   Handle<CaloJetCollection> jetCands;
   iEvent.getByLabel( theJetCollection , jetCands);


   // event weights (CSA07 soups...)   
   eventWeight = -1;

   Handle< double> weightHandle;
   Handle<int> genProcessID;
   Handle<int> alpgenProcessID;
   Handle<double> genEventScale;

   if (getEventWeight) {

     iEvent.getByLabel (theWeightCollection, weightHandle);
     if (weightHandle.isValid()) {
       eventWeight = * weightHandle;
     }

     iEvent.getByLabel( "genEventProcID", genProcessID );
       if (genProcessID.isValid()) {
         processID = *genProcessID;
     }

     if (processID==4 && getAlpgenID) {
             iEvent.getByLabel( alpgenIDTag , alpgenProcessID );
             if (alpgenProcessID.isValid()) {
                 alpgenID = *alpgenProcessID;
             }
     } else alpgenID=0;


     iEvent.getByLabel( "genEventScale", genEventScale );
       if (genEventScale.isValid()) {
         eventScale = *genEventScale;
     }
   }
 
   
    // get hold of TriggerResults
   Handle<TriggerResults> HLTR;  //moja promjena: ovo prebaceno 2 reda nize!
   if (storeHLTresults) {
     //      Handle<TriggerResults> HLTR;
     iEvent.getByLabel(triggerResultsTag,HLTR);
     if (firstTriggerResult) {
       firstTriggerResult = false;
       trigNames.init((edm::TriggerResults&)*HLTR);
       numTriggers = trigNames.size();
     }
   }   
   

   // MET corrected (P_muons, P_muons_calo_deposit)
   
   reco::Particle::LorentzVector MET=computeMET(iEvent, iSetup, looseMuonCands); 
   //     MET_energy=MET.energy();
   MET_energy=MET.energy(); // this is Et!!
   MET_pt=MET.pt();
   MET_px=MET.px();
   MET_py=MET.py();
   MET_eta=MET.eta();
   MET_phi=MET.phi();
   
   //Correct MET due to electron energy scale corrections
   reco::Particle::LorentzVector MET_SCCorr = correctMETforSCCorrection(MET, looseElectronCands);
   MET_Et_SCCorr = MET_SCCorr.energy();
   MET_pt_SCCorr = MET_SCCorr.pt();
   MET_phi_SCCorr = MET_SCCorr.phi();
   
   reco::Particle::LorentzVector MET_ElCorr = correctMETforElCorrection(MET, looseElectronCands);
   MET_Et_ElCorr = MET_ElCorr.energy();
   MET_pt_ElCorr = MET_ElCorr.pt();
   MET_phi_ElCorr = MET_ElCorr.phi();

   
   // selected muons

//    cout << "\t # loose mu : " <<  looseMuonCands->size() 
//      << "\t # tight mu : " <<  goodMuonCands->size() 
//      << "\t # loose e : " <<   looseElectronCands->size() 
//      << "\t # tight e : " <<   allElectronCands->size() 
//      << "\t # tight l : " <<   tightLeptonCands->size() 
//      << "\t # Z->mumu : " <<  zmumuCands->size() 
//      << "\t # Z->ee : " <<  zeeCands->size() 
//      << "\t # Z->ll : " <<  zllCands->size() 
//      << endl;

   N_looseMuons=0;
   N_looseElectrons=0;
   N_looseMuons = looseMuonCands->size();
   N_looseElectrons = looseElectronCands->size();
   NCaloJets = jetCands->size();


   //
   // Find best Z (in Z->ll collection (merged Z->ee & Z->mumu)
   //
   ::OverlapChecker overlap;


  float dzmass_min = 100.;

  n=1;  // number of non overlaping Zs
  for(CandidateCollection::const_iterator z1 = zllCands->begin(); 
      z1 != zllCands->end(); ++ z1 ) {

    // Check that two leptons used to build the Z are not overlapping
    if (overlap(*(z1->daughter(0)),*(z1->daughter(1))) ) {
      continue;
    }

    float dzmass = fabs( z1->mass() - 91.188);
    if (dzmass < dzmass_min) {
      theZCandidate = &(*z1);
      dzmass_min = dzmass;
    }
    //    
    //Get back Electrons from Z and W
    for(CandidateCollection::const_iterator z2 = z1;
        z2 != zllCands->end(); ++ z2 ) {
      if (z1 != z2) {
        if (overlap(*z1,*z2)) {
        } else {
          n++;
        }
      }
    }
  }


  zFlavour = 0;
  wlFlavour = 0;
  zMass = -10;
  zPt  = -10;
  zEta = -10;
  zPhi = -10;
  nWlElCand = 0;
  nWlMuCand = 0;

  dPhi_WlZ_reco=-15.;    
  dPhi_WZ_reco=-15.;
  
  dPhi_Wl_MET = -15.; 
  dPhi_Wl_MET_SCCorr = -15.;
  dPhi_Wl_MET_ElCorr = -15.;
  wlCharge =-10;
  wlEta  = -10;
  wlPhi = -10;
  wlTransMass = -10;
  wlTransMass_SCCorr = -10;
  wlTransMass_ElCorr = -10;
  
  //---------------
//   matching(iEvent,looseMuonCands, 13);
//   matching(iEvent,looseElectronCands, 11);
//   MatchedGenParticle(&(*looseMuonCands)[0]);
  //-------------------

  int nwl_candidates=0;
  if (theZCandidate) {


    zMass      = theZCandidate->mass();
    zFlavour  =  abs(theZCandidate->daughter(0)->pdgId());
    zPt       = theZCandidate->pt();
    zEta     = theZCandidate->eta();
    zPhi      = theZCandidate->phi();

    // Sanity check: can we identify the Z daughters with leptons from the loose lists

    Handle<CandidateCollection> looseLeptonCands;
    if (zFlavour == 11) {
      looseLeptonCands = looseElectronCands;      
    } else if (zFlavour == 13) {
      looseLeptonCands = looseMuonCands;      
    }
    for (int i=0; i<2; i++) {
      int noverlaps=0;
      for(CandidateCollection::const_iterator l = looseLeptonCands->begin();
          l != looseLeptonCands->end(); l++) {
        if (overlap(*(theZCandidate->daughter(i)), *l) ) {
          noverlaps++;
        }
      }
    }

    matching(iEvent,looseMuonCands, 13);
    matching(iEvent,looseElectronCands, 11);
    
    if (zFlavour == 11) {
      //--------------   
      leptonProperties["ZEl1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup, 
                                               MatchedGenParticle(theZCandidate->daughter(0)),dR(theZCandidate->daughter(0)));
      leptonProperties["ZEl2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup, 
                                               MatchedGenParticle(theZCandidate->daughter(1)),dR(theZCandidate->daughter(1)));
    } else if (zFlavour == 13) {
      //--------------   
      const reco::Candidate * mcpart = MatchedGenParticle(theZCandidate->daughter(0));
      //--------------   
      leptonProperties["Zmu1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup, 
                                               MatchedGenParticle(theZCandidate->daughter(0)),dR(theZCandidate->daughter(0)));
      leptonProperties["Zmul"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup, 
                                               MatchedGenParticle(theZCandidate->daughter(0)),dR(theZCandidate->daughter(0)));
      leptonProperties["Zmu2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup, 
                                               MatchedGenParticle(theZCandidate->daughter(1)),dR(theZCandidate->daughter(1)));

    }

    float max_pt = 0.;

    /// Find W lepton


    // Now find lepton that will be associated to W
    //   among leptons not used for the Z reconstruction
    for(CandidateCollection::const_iterator lepton = tightLeptonCands->begin();
        lepton != tightLeptonCands->end(); lepton++) {


      if ( overlap(*theZCandidate, *lepton) )  continue; // Ignore if lepton used for the Z

      nwl_candidates++;
      if (abs(lepton->pdgId()) == 11) {
        nWlElCand++;
      } else if (abs(lepton->pdgId()) == 13) {
        nWlMuCand++;
      }

      // If more than one candidate: choose leading pt
      if (lepton->pt() > max_pt) {
        theWlepton = &(*lepton);
        max_pt     = lepton->pt();
      }

      if(lepton->hasMasterClone()) {
        //      cout << "SHOUT: TIGHT LEPTON IS SHALLOW COPY !!! (SHOULD NOT BE!) \n";
      }
      //
    }
    if (theWlepton) {
      wlFlavour = theWlepton->pdgId();
      wlCharge  = theWlepton->charge();
      wlPt      = theWlepton->pt();
      wlEta     = theWlepton->eta();
      wlPhi     = theWlepton->phi();
      wlTransMass= kinem::TransvMass(theWlepton->pt(),MET.pt(),
                                     theWlepton->px(), MET.px(),
                                     theWlepton->py(),MET.py());

      if (abs(wlFlavour) == 11) {
        leptonProperties["Wel"]->FillProperties(theWlepton, iEvent, iSetup, MatchedGenParticle(theWlepton),dR(theWlepton));
      } else if (abs(wlFlavour) == 13) {
        leptonProperties["Wmu"]->FillProperties(theWlepton, iEvent, iSetup, MatchedGenParticle(theWlepton),dR(theWlepton));
      }

      dPhi_WlZ_reco=kinem::delta_phi(theWlepton->phi(),theZCandidate->phi());    
      dPhi_WZ_reco=kinem::delta_phi((theWlepton->p4()+MET).phi(),theZCandidate->phi());

      dPhi_Wl_MET = kinem::delta_phi(theWlepton->phi(),MET.phi());
      dPhi_Wl_MET_SCCorr = kinem::delta_phi(theWlepton->phi(),MET_SCCorr.phi());
      wlTransMass_SCCorr= kinem::TransvMass(theWlepton->pt(),MET_SCCorr.pt(),
                                     theWlepton->px(), MET_SCCorr.px(),
                                     theWlepton->py(),MET_SCCorr.py());
      dPhi_Wl_MET_ElCorr = kinem::delta_phi(theWlepton->phi(),MET_ElCorr.phi());
      wlTransMass_ElCorr= kinem::TransvMass(theWlepton->pt(),MET_ElCorr.pt(),
                                     theWlepton->px(), MET_ElCorr.px(),
                                     theWlepton->py(),MET_ElCorr.py());
    }
    
  }
  
  //  Handle <reco::GenParticleCandidateCollection> mccands;
  Handle<CandidateCollection> mccands;
  iEvent.getByLabel( theMCTruthCollection,mccands );
  
  collectMCsummary(mccands);
  
  if (storeHLTresults) {
    
    //def: 1:single electron 2:single relaxed 3:double electron 4: double relaxed 5,6,7,8: muon e+mu: 9
    
    triggerBitmask=0;
    triggerBitmaskQCD=0;

    for (size_t i=0; i<numTriggers; ++i) {

        if (!HLTR->wasrun(i)) {
                continue;
        }
        if (HLTR->error(i) ) {
                continue;
        }
        if (HLTR->accept(i)) {
                //if (triggerNames[i].compare("mcpath")!=0) 
                //TRaccept=true;
        }else {
                continue;
        }
        if (trigNames.triggerName(i).compare("HLT1Electron")==0) triggerBitmask |= 1;
        if (trigNames.triggerName(i).compare("HLT1ElectronRelaxed")==0) triggerBitmask |= 2;
        if (trigNames.triggerName(i).compare("HLT2Electron")==0) triggerBitmask |= 4;
        if (trigNames.triggerName(i).compare("HLT2ElectronRelaxed")==0) triggerBitmask |= 8;

        if (trigNames.triggerName(i).compare("HLT1MuonIso")==0) triggerBitmask |= 16;
        if (trigNames.triggerName(i).compare("HLT1MuonNonIso")==0) triggerBitmask |= 32;
        if (trigNames.triggerName(i).compare("HLT2MuonNonIso")==0) triggerBitmask |= 64;
        if (trigNames.triggerName(i).compare("HLT2MuonZ")==0) triggerBitmask |= 128;

        if (trigNames.triggerName(i).compare("HLTXElectronMuon")==0) triggerBitmask |= 256;
        if (trigNames.triggerName(i).compare("HLTXElectronMuonRelaxed")==0) triggerBitmask |= 512;

        if (trigNames.triggerName(i).compare("CandHLT2MuonIso")==0) triggerBitmask |= 1024;

        //for QCD control sample
        if (trigNames.triggerName(i).compare("HLT1Jet")==0) triggerBitmaskQCD        |= 1;
        if (trigNames.triggerName(i).compare("HLT2Jet")==0) triggerBitmaskQCD        |= 2;
        if (trigNames.triggerName(i).compare("HLT3Jet")==0) triggerBitmaskQCD        |= 4;
        if (trigNames.triggerName(i).compare("HLT4Jet")==0) triggerBitmaskQCD        |= 8;
        if (trigNames.triggerName(i).compare("HLT2JetAco")==0) triggerBitmaskQCD     |= 16;
        if (trigNames.triggerName(i).compare("HLTXMuonJets")==0) triggerBitmaskQCD   |= 32;
        if (trigNames.triggerName(i).compare("HLT1MET")==0) triggerBitmaskQCD         |= 64;
        if (trigNames.triggerName(i).compare("HLT1Jet1METAco")==0) triggerBitmaskQCD |= 128;
        if (trigNames.triggerName(i).compare("HLT1Jet1MET")==0) triggerBitmaskQCD    |= 256;
        if (trigNames.triggerName(i).compare("HLT2Jet1MET")==0) triggerBitmaskQCD    |= 512;
        if (trigNames.triggerName(i).compare("HLT3Jet1MET")==0) triggerBitmaskQCD    |= 1024;
        if (trigNames.triggerName(i).compare("HLT4Jet1MET")==0) triggerBitmaskQCD    |= 2048;
        if (trigNames.triggerName(i).compare("HLT1MET1HT")==0) triggerBitmaskQCD     |= 4096;
        if (trigNames.triggerName(i).compare("HLTS2JetAco")==0) triggerBitmaskQCD       |= 8192;

     }
   }
  
  // Jet properties

  float max_et = 0.;
  float max2_et = 0.;
  const Candidate * leadingJet   = 0;
  const Candidate * secondJet    = 0;

  bool matchedLept[3]={false,false,false};

  //  for(CandidateCollection::const_iterator jet = jetCands->begin();
  NbJets = 0;
  NbJetsWithLeptons = 0;

  for(CaloJetCollection::const_iterator jet = jetCands->begin();
      jet != jetCands->end(); jet++) {

        //ignore jets which are within dR 0.15 to a WZ lepton
   {

     float dPhi[3]={-1.0,-1.0,-1.0};
     float dEta[3]={-1.0,-1.0,-1.0};
     bool skipJet=false;

     for (int ct=0;ct<2;ct++) {
             if (theZCandidate)
                     if (theZCandidate->daughter(ct)->pt()>0 
                         && (abs(theZCandidate->daughter(ct)->pdgId())==13 
                             || abs(theZCandidate->daughter(ct)->pdgId())==11) ) {
                       dPhi[ct]=kinem::delta_phi(theZCandidate->daughter(ct)->phi(),jet->phi());
                       dEta[ct]=fabs(theZCandidate->daughter(ct)->eta() - jet->eta());
                     }
     }
     if (theWlepton)
             if (theWlepton->pt()>0 && (abs(theWlepton->pdgId())==13 || abs(theWlepton->pdgId())==11 )) {
               dPhi[2]=kinem::delta_phi(theWlepton->phi(),jet->phi());
               dEta[2]=fabs(theWlepton->eta() - jet->eta());
             }
     for (int k=0;k<3;k++) {
             if (matchedLept[k] || dEta[k]<0.0) continue;

                float deltaR=sqrt(dPhi[k]*dPhi[k]+dEta[k]*dEta[k]);
                if (deltaR<0.15) {
                        skipJet=true;
                        matchedLept[k]=true;
                        break;
                }
        }
        if(jet->pt() > 10 && fabs(jet->eta())<2.5 )NbJetsWithLeptons++;

        if (skipJet || fabs(jet->eta()) > 3.) continue;

        if(jet->pt() > 10 && fabs(jet->eta())<2.5 )NbJets++;
   }

   if (jet->et() > max_et ) {
     if (leadingJet) {
       secondJet = leadingJet;
       max2_et   = max_et;
     }
     leadingJet = &(*jet);
     max_et     = jet->et();
     
   } else {
     if (jet->et() > max2_et ) {
       secondJet = &(*jet);
       max2_et   = jet->et();
     }
   }
  }

  if (leadingJet) {
          jet1Et  = leadingJet->et();
          jet1Phi = leadingJet->phi();
          jet1Eta = leadingJet->eta();

  } else {
          jet1Et  = -10;
          jet1Phi = 0.;
          jet1Eta = 0.;
  }
  if (secondJet) {
          jet2Et  = secondJet->et();
          jet2Phi = secondJet->phi();
          jet2Eta = secondJet->eta();
  } else {
          jet2Et  = -10;
          jet2Phi = 0.;
          jet2Eta = 0.;

  }

  // matching:

  matching(iEvent,looseMuonCands, 13);
  matching(iEvent,looseElectronCands, 11);

  getMother(&(*mccands)[1]);

  numberOfElectronTrees  = 0;

  //count number of electron treee entries which will be added
  for(CandidateCollection::const_iterator el = looseElectronCands->begin(); 
      el != looseElectronCands->end(); ++ el ) {
     if (el->pt()>10 && fabs(el->eta())<2.5) {
                numberOfElectronTrees++;
     }
  }

  wzTree->Fill();


  // Fill electrons in a tree
  for(CandidateCollection::const_iterator el = looseElectronCands->begin(); 
      el != looseElectronCands->end(); ++ el ) {

    electronUse = 0;

    if (theZCandidate) {

      if (overlap(*theZCandidate, *el) ) {
        electronUse = 23;
      } else if (theWlepton) {
        if (overlap(*theWlepton, *el) ) {
          electronUse = 24;
        }
      }
    }
    if (el->pt()>10 && fabs(el->eta())<2.5) {
      leptonProperties["electron"]->FillProperties(&(*el), 
                                                   iEvent, iSetup, 
                                                   MatchedGenParticle(&(*el)),dR(&(*el)));

      if (findBackJets) {

        backJet_found=0;
        backJet_Et=0;
        backJet_Eta=0;
        backJet_Phi=0;
        backJet_maxEInEmTowers=0;
        backJet_maxEInHadTowers=0;
        backJet_energyFractionHadronic=0;
        backJet_emEnergyFraction=0;
        backJet_hadEnergyInHB=0;
        backJet_hadEnergyInHO=0;
        backJet_hadEnergyInHE=0;
        backJet_hadEnergyInHF=0;
        backJet_emEnergyInEB=0;
        backJet_emEnergyInEE=0;
        backJet_emEnergyInHF=0;
        backJet_towersArea=0;
        backJet_n90=0;
        backJet_n60=0;
        backJet_NbTracks = 0;
        backJet_match=0;
        backJet_matchPt=0;
        backJet_match_dR=-1.0;

        float oldDeltaPhi = -1.0;

        edm::Handle<reco::TrackCollection> BarrelKFTracks;
        iEvent.getByLabel("ctfWithMaterialTracks",BarrelKFTracks);

        for (CaloJetCollection::const_iterator jet = jetCands->begin();
                        jet != jetCands->end(); jet++) {

                int NbTracksInJets = 0;

                for( reco::TrackCollection::const_iterator KFT = BarrelKFTracks->begin();KFT !=BarrelKFTracks->end(); ++ KFT ){

                  double DPhi = kinem::delta_phi(KFT->phi(),jet->phi());
                
                        double dR=sqrt( (KFT->eta() - jet->eta())*(KFT->eta() - jet->eta())   + DPhi*DPhi );

                        if (KFT->pt()>2.0 && dR<0.3) NbTracksInJets++;

                        //cout << "pass track eta:" << KFT->eta() << "  pt " << KFT->pt() << "  dr " << dR << endl;

                }

                if (NbTracksInJets<=3) continue;

                float dPhi = fabs(kinem::delta_phi(jet->phi(),el->phi()) - kinem::PI);

                if (dPhi<0.7 && (dPhi<oldDeltaPhi || oldDeltaPhi<0.0)) {
                        oldDeltaPhi=dPhi;
 
                        backJet_found=1;
                        backJet_Et=jet->et();
                        backJet_Eta=jet->eta();
                        backJet_Phi=jet->phi();
                        backJet_maxEInEmTowers=jet->maxEInEmTowers();
                        backJet_maxEInHadTowers=jet->maxEInHadTowers();
                        backJet_energyFractionHadronic=jet->energyFractionHadronic();
                        backJet_emEnergyFraction=jet->emEnergyFraction();
                        backJet_hadEnergyInHB=jet->hadEnergyInHB();
                        backJet_hadEnergyInHO=jet->hadEnergyInHO();
                        backJet_hadEnergyInHE=jet->hadEnergyInHE();
                        backJet_hadEnergyInHF=jet->hadEnergyInHF();
                        backJet_emEnergyInEB=jet->emEnergyInEB();
                        backJet_emEnergyInEE=jet->emEnergyInEE();
                        backJet_emEnergyInHF=jet->emEnergyInHF();
                        backJet_towersArea=jet->towersArea();
                        backJet_n90=jet->n90();
                        backJet_n60=jet->n60();
                        //fill variables
                        backJet_NbTracks = NbTracksInJets;


                        //find monte carlo electrons close to the jet

                        float eDR = -1.0;
                }
        }
      }

      electronTree->Fill();
      
    }
    
  }

  EventCounter++;
}


// ------------ method called once each job just before starting event loop  ------------
void 
WZAnalyzer::beginJob(const edm::EventSetup&)
{

  using namespace wzana;

  //  theFile   = new TFile( "wz.root", "RECREATE" ) ;
  theFile   = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;


  wzTree = new TTree("WZTree","WZ informations");

  electronTree = new TTree("ElTree","electron informations");

//  jetTree = new TTree("JetTree","jet deltaR information");

  leptonProperties["electron"] = new ElectronProperties();

  leptonProperties["Wel"]  = new ElectronProperties();
  leptonProperties["ZEl1"] = new ElectronProperties();
  leptonProperties["ZEl2"] = new ElectronProperties();

  leptonProperties["Wmu"]  = new MuonProperties();
  leptonProperties["Zmu1"] = new MuonProperties();
  leptonProperties["Zmu2"] = new MuonProperties();

  leptonProperties["Zmul"] = new MuonProperties();       // SAME COMMENT AS BELOW: TEMPORARY FIX

  leptonProperties["Wel"]->CreateBranches(wzTree, "WEl_");
  leptonProperties["ZEl1"]->CreateBranches(wzTree, "ZEl1_");
  leptonProperties["ZEl2"]->CreateBranches(wzTree, "ZEl2_");

  leptonProperties["Wmu"]->CreateBranches(wzTree, "Wmu_");

  leptonProperties["Zmu1"]->CreateBranches(wzTree, "Zmu1_");
  // Vuko's NOTE: there was type in the key "Zmul" instead of "Zmu1"
  // ("L" instead of "one")
  // I am fixing it but in order not to break everybody's root scripts
  // I am adding the next line so we have both until all trees
  // are contain "Zmu1"
  // when that is the case, the next line can be dropped
  leptonProperties["Zmul"]->CreateBranches(wzTree, "Zmul_");
  leptonProperties["Zmu2"]->CreateBranches(wzTree, "Zmu2_");
  leptonProperties["electron"]->CreateBranches(electronTree);

  initialiseTree();
  
  //prepare HLT TriggerResults branch
  if (storeHLTresults) {
        wzTree->Branch("WZ_hltBitmask",&triggerBitmask,"WZ_hltBitmask/I");
        wzTree->Branch("QCD_hltBitmask",&triggerBitmaskQCD,"QCD_hltBitmask/I");

  }

  // MET branch
  wzTree->Branch("MET_Et",&MET_energy,"MET_Et/F");
  wzTree->Branch("MET_pt",&MET_pt,"MET_pt/F");
  wzTree->Branch("MET_eta",&MET_eta,"MET_eta/F");
  wzTree->Branch("MET_phi",&MET_phi,"MET_eta/F");
  wzTree->Branch("dPhi_WlZ_reco",&dPhi_WlZ_reco,"dPhi_WlZ_reco/F");
  wzTree->Branch("dPhi_WZ_reco",&dPhi_WZ_reco,"dPhi_WZ_reco/F");
  wzTree->Branch("dPhi_Wl_MET",&dPhi_Wl_MET,"dPhi_Wl_MET/F");

  // MET branch after corrections
  wzTree->Branch("MET_Et_SCCorr",&MET_Et_SCCorr,"MET_Et_SCCorr/F");
  wzTree->Branch("MET_pt_SCCorr",&MET_pt_SCCorr,"MET_pt_SCCorr/F");
  wzTree->Branch("MET_phi_SCCorr",&MET_phi_SCCorr,"MET_phi_SCCorr/F");
  wzTree->Branch("dPhi_Wl_MET_SCCorr",&dPhi_Wl_MET_SCCorr,"dPhi_Wl_MET_SCCorr/F");

  wzTree->Branch("MET_Et_ElCorr",&MET_Et_ElCorr,"MET_Et_SCCorr/F");
  wzTree->Branch("MET_pt_ElCorr",&MET_pt_ElCorr,"MET_pt_SCCorr/F");
  wzTree->Branch("MET_phi_ElCorr",&MET_phi_ElCorr,"MET_phi_SCCorr/F");
  wzTree->Branch("dPhi_Wl_MET_ElCorr",&dPhi_Wl_MET_ElCorr,"dPhi_Wl_MET_ElCorr/F");

  wzTree->Branch("WlTransMass",     &wlTransMass, "WlTransMass/F");
  wzTree->Branch("WlTransMass_SCCorr", &wlTransMass_SCCorr,  "WlTransMass_SCCorr/F");
  wzTree->Branch("WlTransMass_ElCorr",     &wlTransMass_ElCorr,      "WlTransMass_ElCorr/F");


  // # of non overlaping Zs
  wzTree->Branch("N_nonOverlaping_Z",&n,"N_nonOverlaping_Z/I");   
  wzTree->Branch("N_looseMuons",&N_looseMuons,"N_looseMuons/I");   
  wzTree->Branch("N_looseElectrons",&N_looseElectrons,"N_looseElectrons/I");   


  EventCounter=0;
}

// ------------ method called once each job just after ending the event loop  ------------
void 
WZAnalyzer::endJob() {

  theFile->Write();
  theFile->Close();

}


void WZAnalyzer::initialiseTree() {
  // Event properties
  wzTree->Branch("weight",    &eventWeight,  "weight/F");
  wzTree->Branch("processID", &processID, "processID/I");
  wzTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
  wzTree->Branch("genEventScale",    &eventScale,  "genEventScale/F");
  wzTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
  wzTree->Branch("EventID", &EventID, "EventID/I");

  // Z properties  
  wzTree->Branch("Zmass",    &zMass,    "Zmass/F");
  wzTree->Branch("ZId",      &zFlavour, "Zid/I");
  wzTree->Branch("ZPt",      &zPt,      "ZPt/F");
  wzTree->Branch("ZEta",     &zEta,     "ZEta/F");
  wzTree->Branch("ZPhi",     &zPhi,     "ZPhi/F");

  // W Properties
  wzTree->Branch("WlId",     &wlFlavour, "Wlid/I");
  wzTree->Branch("WlCharge", &wlCharge,  "WlCharge/I");
  wzTree->Branch("WlPt",     &wlPt,      "WlPt/F");

  wzTree->Branch("nwlelcand",     &nWlElCand,      "nwlelcand/I");
  wzTree->Branch("nwlmucand",     &nWlMuCand,      "nwlmucand/I");

  
  wzTree->Branch("MC_leptonZ1_pt", &MCleptZ1_pt,"MC_leptonZ1_pt/F");
  wzTree->Branch("MC_leptonZ2_pt", &MCleptZ2_pt,"MC_leptonZ2_pt/F");
  wzTree->Branch("MC_leptonW_pt",  &MCleptW_pt, "MC_leptonW_pt/F");

  wzTree->Branch("MC_leptonZ1_eta", &MCleptZ1_eta,"MC_leptonZ1_eta/F");
  wzTree->Branch("MC_leptonZ2_eta", &MCleptZ2_eta,"MC_leptonZ2_eta/F");
  wzTree->Branch("MC_leptonW_eta",  &MCleptW_eta, "MC_leptonW_eta/F");

  wzTree->Branch("MC_leptonZ1_phi", &MCleptZ1_phi,"MC_leptonZ1_phi/F");
  wzTree->Branch("MC_leptonZ2_phi", &MCleptZ2_phi,"MC_leptonZ2_phi/F");
  wzTree->Branch("MC_leptonW_phi",  &MCleptW_phi, "MC_leptonW_phi/F");

  wzTree->Branch("MC_leptonZ1_pdgid", &MCleptZ1_pdgid,"MC_leptonZ1_pdgid/I");
  wzTree->Branch("MC_leptonZ2_pdgid", &MCleptZ2_pdgid,"MC_leptonZ2_pdgid/I");
  wzTree->Branch("MC_leptonW_pdgid",  &MCleptW_pdgid, "MC_leptonW_pdgid/I");

  /*
  wzTree->Branch("MC_leptonW_origin", &MCleptW_pdgid, "MC_leptonW_origin/I");
  wzTree->Branch("MC_leptonZ1_origin", &MCleptW_pdgid, "MC_leptonZ1_origin/I");
  wzTree->Branch("MC_leptonZ2_origin", &MCleptW_pdgid, "MC_leptonZ2_origin/I");
  */

  wzTree->Branch("MC_TauDecayType_fromZ1", &MC_tauDecayTypeZ1,"MC_TauDecayType_fromZ1/I");
  wzTree->Branch("MC_TauDecayType_fromZ2", &MC_tauDecayTypeZ2,"MC_TauDecayType_fromZ2/I");
  wzTree->Branch("MC_TauDecayType_fromW",  &MC_tauDecayTypeW, "MC_TauDecayType_fromW/I");

   // Jet properties
  wzTree->Branch("NbJets",     &NbJets,  "NbJets/I");
  wzTree->Branch("NCaloJets",  &NCaloJets,"NCaloJets/I");

  wzTree->Branch("Jet1Et",     &jet1Et,  "Jet1Et/F");
  wzTree->Branch("Jet1Phi",    &jet1Phi, "Jet1Phi/F");
  wzTree->Branch("Jet1Eta",    &jet1Eta, "Jet1Eta/F");

  wzTree->Branch("Jet2Et",     &jet2Et,  "Jet2Et/F");
  wzTree->Branch("Jet2Phi",    &jet2Phi, "Jet2Phi/F");
  wzTree->Branch("Jet2Eta",    &jet2Eta, "Jet2Eta/F");
  wzTree->Branch("ElectronTreesInEvent", &numberOfElectronTrees, "ElectronTreesInEvent/I");
  wzTree->Branch("EventCounter", &EventCounter, "EventCounter/I");
  //
  electronTree->Branch("used",   &electronUse, "used/I");
  electronTree->Branch("weight",    &eventWeight,  "weight/F");
  electronTree->Branch("processID", &processID, "processID/I");
  electronTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
  electronTree->Branch("genEventScale",    &eventScale,  "genEventScale/F");
  electronTree->Branch("ElectronTreesInEvent", &numberOfElectronTrees, "ElectronTreesInEvent/I");

  electronTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
  electronTree->Branch("EventID", &EventID, "EventID/I");
  electronTree->Branch("EventCounter", &EventCounter, "EventCounter/I");
  electronTree->Branch("NbJetsWithLeptons",     &NbJetsWithLeptons,  "NbJetsWithLeptons/I");
 
  if (findBackJets) {
    electronTree->Branch("MET_Et",&MET_energy,"MET_Et/F");
    electronTree->Branch("MET_pt",&MET_pt,"MET_pt/F");
    electronTree->Branch("MET_eta",&MET_eta,"MET_eta/F");
    electronTree->Branch("MET_phi",&MET_phi,"MET_phi/F");

    electronTree->Branch("backJet_found",&backJet_found,"backJet_found/I");
    electronTree->Branch("backJet_Et",&backJet_Et,"backJet_Et/F");
    electronTree->Branch("backJet_Eta",&backJet_Eta,"backJet_Eta/F");
    electronTree->Branch("backJet_Phi",&backJet_Phi,"backJet_Phi/F");
    electronTree->Branch("backJet_maxEInEmTowers",&backJet_maxEInEmTowers,"backJet_maxEInEmTowers/F");
    electronTree->Branch("backJet_maxEInHadTowers",&backJet_maxEInHadTowers,"backJet_maxEInHadTowers/F");
    electronTree->Branch("backJet_energyFractionHadronic",&backJet_energyFractionHadronic,"backJet_energyFractionHadronic/F");
    electronTree->Branch("backJet_emEnergyFraction",&backJet_emEnergyFraction,"backJet_emEnergyFraction/F");
    electronTree->Branch("backJet_hadEnergyInHB",&backJet_hadEnergyInHB,"backJet_hadEnergyInHB/F");
    electronTree->Branch("backJet_hadEnergyInHO",&backJet_hadEnergyInHO,"backJet_hadEnergyInHO/F");
    electronTree->Branch("backJet_hadEnergyInHE",&backJet_hadEnergyInHE,"backJet_hadEnergyInHE/F");
    electronTree->Branch("backJet_hadEnergyInHF",&backJet_hadEnergyInHF,"backJet_hadEnergyInHF/F");
    electronTree->Branch("backJet_emEnergyInEB",&backJet_emEnergyInEB,"backJet_emEnergyInEB/F");
    electronTree->Branch("backJet_emEnergyInEE",&backJet_emEnergyInEE,"backJet_emEnergyInEE/F");
    electronTree->Branch("backJet_emEnergyInHF",&backJet_emEnergyInHF,"backJet_emEnergyInHF/F");
    electronTree->Branch("backJet_towersArea",&backJet_towersArea,"backJet_towersArea/F");
    electronTree->Branch("backJet_n90",&backJet_n90,"backJet_n90/F");
    electronTree->Branch("backJet_n60",&backJet_n60,"backJet_n60/F");

    electronTree->Branch("backJet_NbTracks",&backJet_NbTracks,"backJet_NbTracks/I");

    electronTree->Branch("backJet_match",&backJet_match,"backJet_match/I");
    electronTree->Branch("backJet_matchPt",&backJet_matchPt,"backJet_matchPt/F");
    electronTree->Branch("backJet_match_dR",&backJet_match_dR,"backJet_match_dR/F");


  }

  electronTree->Branch("WZ_hltBitmask",&triggerBitmask,"WZ_hltBitmask/I");
  electronTree->Branch("QCD_hltBitmask",&triggerBitmaskQCD,"QCD_hltBitmask/I");
}


    //////////////////////////////////////////
    // MET & MET corrections
    // MET=MET-sum_p(muons)-sum_p(muons calos depositions)

reco::Particle::LorentzVector WZAnalyzer::computeMET(const edm::Event& iEvent, const edm::EventSetup& iSetup, Handle<reco::CandidateCollection> MuonCands) {
  
  using namespace edm;
  using namespace reco;
  using namespace std;

  double px_MET=0, py_MET=0, Ex_MET=0, Ey_MET=0;
  double px,py,pz,pt,Et,Ex,Ey;
  double px_muons=0;
  double py_muons=0;
  double Ex_muons=0;
  double Ey_muons=0;
  double Et_muons=0;   
  double Et_ecal; 
  double Et_hcal; 
  double Et_ho; 
  double Et_cal;  
  double Ex_ecal; 
  double Ex_hcal; 
  double Ex_ho; 
  double Ey_ecal; 
  double Ey_hcal; 
  double Ey_ho; 
  double eta_cal; 
  double phi_cal; 
  double tg_theta_pola_ecal; 
  double tg_theta_pola_hcal; 
  double tg_theta_pola_ho; 
  double tg_theta_ecal; 
  double tg_theta_hcal; 
  double tg_theta_ho; 
  double E_ecal; 
  double E_hcal; 
  double E_ho; 
  double me=0.0005; 
  double sin_theta_ecal, sin_theta_hcal, sin_theta_ho, cos_theta_ecal,cos_theta_hcal,cos_theta_ho;
  double p_ecal, p_hcal,p_ho, pt_ecal,pt_hcal,pt_ho, px_ecal,px_hcal,px_ho, py_ecal,py_hcal,py_ho, pz_ecal,pz_hcal,pz_ho;
  double Ex_cal=0; 
  double Ey_cal=0;    
  double px_cal=0; 
  double py_cal=0;
  double Ex_kon=0;
  double Ey_kon=0;

  reco::Particle::LorentzVector Muon_p_4v(0,0,0,0);  // sum of muons 4-vectors in (px,py,0,Et), since MET looks like that
  reco::Particle::LorentzVector Muon_p_cal_4v(0,0,0,0); 
  
  Handle<CaloMETCollection> mets;
  iEvent.getByLabel("met", mets);
  //  cout << "# METs=" << mets->size() << endl;
  
  if (mets->size() != 1) {
    cout << "ALARM: # METS is not 1: !" << endl;
    return reco::Particle::LorentzVector(0.,0.,0.,0.);
  }
  
  for (CaloMETCollection::const_iterator met = mets->begin();
       met!= mets->end(); met++) {
    px_MET=met->px();
    py_MET=met->py();
    Ex_MET=met->et()*TMath::Cos(met->phi());
    Ey_MET=met->et()*TMath::Sin(met->phi());

  }
  
  // sum of p of muons
  for (unsigned int i=0; i<MuonCands->size();i++){
    px=(*MuonCands)[i].px();
    py=(*MuonCands)[i].py();
    pz=(*MuonCands)[i].pz();
    pt=(*MuonCands)[i].pt();
    Et=(*MuonCands)[i].et();
    Ex=Et*TMath::Cos((*MuonCands)[i].phi());
    Ey=Et*TMath::Sin((*MuonCands)[i].phi());
    px_muons=px_muons+px;
    py_muons=py_muons+py;
    Ex_muons=Ex_muons+Ex;
    Ey_muons=Ey_muons+Ey;
   
  }
  Et_muons=TMath::Sqrt(Ex_muons*Ex_muons+Ey_muons*Ey_muons);
  Muon_p_4v=reco::Particle::LorentzVector(px_muons,py_muons,0,Et_muons);
  
// sum of p of muons left in calorimeters  
  for (unsigned int i=0; i<MuonCands->size();i++){
    cout<<"" << endl;
    const reco::Muon* muons=dynamic_cast<const reco::Muon *> (&(*MuonCands)[i]);
    Et_ecal=muons->getCalEnergy().emS9;  
    Et_hcal=muons->getCalEnergy().hadS9;
    Et_ho=muons->getCalEnergy().hoS9;  
    eta_cal=muons->eta();  
    phi_cal=muons->phi();  
    tg_theta_pola_ecal=TMath::Exp(-eta_cal);
    tg_theta_pola_hcal=TMath::Exp(-eta_cal);      
    tg_theta_pola_ho=TMath::Exp(-eta_cal);
    tg_theta_ecal=2*tg_theta_pola_ecal/(1-tg_theta_pola_ecal*tg_theta_pola_ecal);
    tg_theta_hcal=2*tg_theta_pola_hcal/(1-tg_theta_pola_hcal*tg_theta_pola_hcal);      
    tg_theta_ho=2*tg_theta_pola_ho/(1-tg_theta_pola_ho*tg_theta_pola_ho);
    E_ecal=TMath::Abs(Et_ecal*TMath::Sqrt(1+tg_theta_ecal*tg_theta_ecal)/tg_theta_ecal);
    E_hcal=TMath::Abs(Et_hcal*TMath::Sqrt(1+tg_theta_hcal*tg_theta_hcal)/tg_theta_hcal);
    E_ho=TMath::Abs(Et_ho*TMath::Sqrt(1+tg_theta_ho*tg_theta_ho)/tg_theta_ho);
    Ex_ecal=Et_ecal*TMath::Cos(phi_cal);
    Ex_hcal=Et_hcal*TMath::Cos(phi_cal);
    Ex_ho=Et_ho*TMath::Cos(phi_cal);
    Ey_ecal=Et_ecal*TMath::Sin(phi_cal);
    Ey_hcal=Et_hcal*TMath::Sin(phi_cal);
    Ey_ho=Et_ho*TMath::Sin(phi_cal);
    sin_theta_ecal=tg_theta_ecal/TMath::Sqrt(1+tg_theta_ecal*tg_theta_ecal);
    sin_theta_hcal=tg_theta_hcal/TMath::Sqrt(1+tg_theta_hcal*tg_theta_hcal);
    sin_theta_ho=tg_theta_ho/TMath::Sqrt(1+tg_theta_ho*tg_theta_ho);
    cos_theta_ecal=TMath::Sqrt(1-sin_theta_ecal*sin_theta_ecal);
    cos_theta_hcal=TMath::Sqrt(1-sin_theta_hcal*sin_theta_hcal);
    cos_theta_ho=TMath::Sqrt(1-sin_theta_ho*sin_theta_ho);
    if (E_ecal<me){
      p_ecal=0;
    }
    else {
      p_ecal=TMath::Sqrt(E_ecal*E_ecal-me*me);
    }
    if (E_hcal<me){
      p_hcal=0;
    }
    else {
      p_hcal=TMath::Sqrt(E_hcal*E_hcal-me*me);
    }      
    if (E_ho<me){
      p_ho=0;
    }
    else {
      p_ho=TMath::Sqrt(E_ho*E_ho-me*me);
    }

    pz_ecal=p_ecal*cos_theta_ecal;
    pz_hcal=p_hcal*cos_theta_hcal;
    pz_ho=p_ho*cos_theta_ho;
    pt_ecal=TMath::Abs(p_ecal*sin_theta_ecal);
    
    pt_hcal=TMath::Abs(p_hcal*sin_theta_hcal);
    
    pt_ho=TMath::Abs(p_ho*sin_theta_ho);      
    
    px_ecal=pt_ecal*TMath::Cos(phi_cal);
    px_hcal=pt_hcal*TMath::Cos(phi_cal);
    px_ho=pt_ho*TMath::Cos(phi_cal);
    py_ecal=pt_ecal*TMath::Sin(phi_cal);
    py_hcal=pt_hcal*TMath::Sin(phi_cal);
    py_ho=pt_ho*TMath::Sin(phi_cal);      

    
    Ex_cal=Ex_cal+Ex_ecal+Ex_hcal+Ex_ho;
    Ey_cal=Ey_cal+Ey_ecal+Ey_hcal+Ey_ho;
    
    px_cal=px_cal+px_ecal+px_hcal+px_ho;
    py_cal=py_cal+py_ecal+py_hcal+py_ho;
    

  }
  Et_cal=TMath::Sqrt(Ex_cal*Ex_cal+Ey_cal*Ey_cal);
  Muon_p_cal_4v=reco::Particle::LorentzVector(px_cal,py_cal,0, Et_cal);  

  Ex_kon=Ex_MET-Ex_muons+Ex_cal;
  Ey_kon=Ey_MET-Ey_muons+Ey_cal;

  reco::Particle::LorentzVector MET= reco::Particle::LorentzVector(px_MET-px_muons+px_cal,py_MET-py_muons+py_cal,0,TMath::Sqrt(Ex_kon*Ex_kon+Ey_kon*Ey_kon));  // substract p from muons and p of calo muon deposit from MET 

  return MET;
  
}

//MET Corrections from difference between corrected and uncorretced SC
reco::Particle::LorentzVector WZAnalyzer::correctMETforSCCorrection(reco::Particle::LorentzVector MET_p_4v, Handle<reco::CandidateCollection> ElectronCands) {

  double SCRawEnergy = 0;
  double SCEnergy = 0;
  double deltaE = 0;
  double deltaET = 0;
  double deltaEx = 0;
  double deltaEy = 0;
 
  reco::Particle::LorentzVector SCCorr_p_4v(0,0,0,0);  // sum of correction fraction of ET in (px,py,0,Et)
  
  //Loop over electron candidates and obtain the difference between uncorr and corr SC
  for (unsigned int i =0; i<ElectronCands->size(); ++i) {
    const reco::PixelMatchGsfElectron * electron = 0;
    if ( (*ElectronCands)[i].hasMasterClone() ) {
      reco::CandidateBaseRef master = (*ElectronCands)[i].masterClone();
      electron = dynamic_cast<const reco::PixelMatchGsfElectron *>(&(*master));
    } else {
    electron = dynamic_cast<const reco::PixelMatchGsfElectron *>(&(*ElectronCands)[i]);
    }
    const reco::SuperCluster& ElSC = *(electron->superCluster());
    SCRawEnergy = ElSC.rawEnergy();
    SCEnergy = ElSC.energy();
    //b.c. both uncrected and corrected SC clusters are oriented the same way, difference between vectors is just magnitudinal.
    deltaE = SCEnergy - SCRawEnergy;
    deltaET = deltaE/cosh(ElSC.eta());
    deltaEx += deltaET*cos(ElSC.phi());
    deltaEy += deltaET*sin(ElSC.phi());
  }

  deltaET = 0;
  deltaET = sqrt(deltaEx*deltaEx + deltaEy*deltaEy);
  SCCorr_p_4v = reco::Particle::LorentzVector(deltaEx, deltaEy, 0, deltaET);
  //  return SCCorr_p_4v;
  return (MET_p_4v - SCCorr_p_4v);
}


//MET Corrections from difference between el energy and rawenergy of SC
reco::Particle::LorentzVector WZAnalyzer::correctMETforElCorrection(reco::Particle::LorentzVector MET_p_4v, Handle<reco::CandidateCollection> ElectronCands) {

  double px = 0;
  double py = 0;
  double pz = 0;
  double Ex = 0;
  double Ey = 0;
  double Et = 0;
  double SCRawEnergy = 0;
  double SCPt = 0;
  double SCEt = 0;
  double SCpx = 0; 
  double SCpy = 0; 
  double SCEx = 0;
  double SCEy = 0;
  
  const reco::PixelMatchGsfElectron * electron = 0;

  reco::Particle::LorentzVector El_p_4v(0,0,0,0);   //  in (px,py,0,Et)
  reco::Particle::LorentzVector SC_p_4v(0,0,0,0);

  //Loop over electron candidates and obtain the difference between uncorr and corr SC
  for (unsigned int i =0; i<ElectronCands->size(); ++i) {
    
    px += (*ElectronCands)[i].px();
    py += (*ElectronCands)[i].py();
    pz += (*ElectronCands)[i].pz();
    Ex += (*ElectronCands)[i].et()*cos((*ElectronCands)[i].phi());
    Ey += (*ElectronCands)[i].et()*sin((*ElectronCands)[i].phi());
    
    if ( (*ElectronCands)[i].hasMasterClone() ) {
      reco::CandidateBaseRef master = (*ElectronCands)[i].masterClone();
      electron = dynamic_cast<const reco::PixelMatchGsfElectron *>(&(*master));
    } else {
      electron = dynamic_cast<const reco::PixelMatchGsfElectron *>(&(*ElectronCands)[i]);
    }
    const reco::SuperCluster& ElSC = *(electron->superCluster());
    SCRawEnergy = ElSC.rawEnergy();
    SCPt = SCRawEnergy/cosh(ElSC.eta());
    SCEt = SCRawEnergy/cosh(ElSC.eta());
    SCpx += SCPt*cos(ElSC.phi());
    SCpy += SCPt*sin(ElSC.phi());
    SCEx += SCEt*cos(ElSC.phi());
    SCEy += SCEt*sin(ElSC.phi());
  }

  Et = sqrt(Ex*Ex + Ey*Ey);
  El_p_4v = reco::Particle::LorentzVector(px,py,0,Et);
  
  SCEt = sqrt(SCEx*SCEx + SCEy * SCEy);  
  SC_p_4v = reco::Particle::LorentzVector(SCpx,SCpy,0,SCEt);

  reco::Particle::LorentzVector ElCorr_p_4v = El_p_4v - SC_p_4v;
  return (MET_p_4v - ElCorr_p_4v);
}


////////////////////////
// matching
//

class MatchingInfo {
public:

  MatchingInfo(float dr, int gen, int reco) : deltaR(dr), genMuon(gen),recoMuon(reco) {}
  
  bool operator<(const MatchingInfo m) { return deltaR < m.deltaR;}
  
  float deltaR;
  int   genMuon;
  int   recoMuon;
  
};


bool betterMatch(MatchingInfo m1, MatchingInfo m2) 
{ return m1.deltaR < m2.deltaR;}


// method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:

void WZAnalyzer::matching(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, int pdgid){

  using namespace edm;
  using namespace reco;
  using namespace std;

  //-------------------------
  Handle<CandidateCollection> mccands;
  iEvent.getByLabel( theMCTruthCollection,mccands );
  
  vector <GenParticleCandidate*> genparticles;
  
  for(CandidateCollection::const_iterator p = mccands->begin();
      p != mccands->end(); p++) {
    
    //  reco::Candidate* p_tmp = &(*p);
    reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
    
    if ( (ptr)->status() == 1
         /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
      if ( abs((ptr)->pdgId()) == pdgid ) {
        genparticles.push_back((ptr));
        //cout << "electron MCT\n";
      }
    }
  }
  

  int n1=0; 
  vector<MatchingInfo> matching_Infos;
  for ( unsigned int i=0; i<genparticles.size(); i++ ) // po svim MC muonima
    {
      
      for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
        {
          
          double d_eta2=(*cands)[j].eta()-(genparticles)[i]->eta();
          double d_phi=fabs((*cands)[j].phi()-(genparticles)[i]->phi());
          double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
          double dR=dR_byhand;
          
          if (dR<0.15){
            n1++;
            matching_Infos.push_back(MatchingInfo(dR,i,j));
          }
          
        }
    } 
  
  sort(matching_Infos.begin(),matching_Infos.end(),betterMatch);  // sorting (dR,i,j)   
  
  if (genparticles.size()!=0 && cands->size()!=0){
    // Now exclude combinations using same gen or rec muons
    for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
         it != matching_Infos.end(); it++) {
      for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
           it2!=it; it2++) {
        if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
          matching_Infos.erase(it);
          it=matching_Infos.begin();
          break;
        }
      }
    }
  }
  
   // Now put result in vector of pairs

   leptonMatching[pdgid].clear();
   
   for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
        match !=matching_Infos.end(); match++) { // po tablici dr,i,j

     pair<const GenParticleCandidate*, const reco::Candidate *> pair;
     //     vector<pair,float > matchedParticles;
     pair.first = genparticles[match->genMuon];
     pair.second = & (*cands)[match->recoMuon];

     matchedParticles cestice(genparticles[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);

     leptonMatching[pdgid].push_back(cestice);

   }
  
}


// get mother of particle 
// returns mother = 1 if mother is Z boson
// returns mother = 2 if mother is W boson
// returns mother = 4 if mother is b
// returns mother = 0 else 

WZAnalyzer::LeptonOrigin WZAnalyzer::getMother(const reco::Candidate* genParticle){

  int pdg_id=genParticle->pdgId();
  //  cout << "particle id : " << pdg_id << endl;

  while (genParticle = genParticle->mother() ) {  //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
    //    cout << "Going up: ";
    //    cout << "Mother id : " << genParticle->pdgId() << endl;
    if (abs(genParticle->pdgId())!=abs(pdg_id)) {
      //      cout<< " good mother " <<endl;
      if (abs(genParticle->pdgId())==23){ // mother is Z
        return zboson;
      } 
      if (abs(genParticle->pdgId())==24) { // mother is W
        return wboson;
      }
      //       if (abs(genParticle->pdgId())==23 || abs(genParticle->pdgId())==24) { // mother is W or Z
      //        WZ_matched=1; 
      //        mother=3;
      //       } 
      if ((((abs(genParticle->pdgId()))/100)%10 ==5)||(((abs(genParticle->pdgId()))/1000)%10 ==5  )) {  // mother is b
        return bdecay;
      }    
      if ((((abs(genParticle->pdgId()))/100)%10 ==4)||(((abs(genParticle->pdgId()))/1000)%10 ==4  )) {  // mother is c
        return cdecay;
      }    
      return other;
    }
  }

  return other; 
}

const reco::Candidate * WZAnalyzer::MatchedGenParticle(const reco::Candidate * daughter){
  ::OverlapChecker overlap;
  matched_genParticle=0;

  for (map< int, 
         std::vector<  matchedParticles >   > ::iterator iter = leptonMatching.begin();
       iter!=leptonMatching.end();   iter++) // entire map
    {
      for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
        const reco::Candidate * nesto=iter->second[j].pair.second;
        if (overlap(*nesto,*daughter)){
          matched_genParticle=iter->second[j].pair.first;
        }
      }
    }
  return matched_genParticle;
}


float WZAnalyzer::dR(const reco::Candidate * daughter){

  ::OverlapChecker overlap;

  float delta_r = -10;
  for (map< int, 
         std::vector<  matchedParticles >   > ::iterator iter = leptonMatching.begin();
       iter!=leptonMatching.end();   iter++) // entire map
    {
      for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
        const reco::Candidate * nesto=iter->second[j].pair.second;
        if (overlap(*nesto,*daughter)){
          delta_r=iter->second[j].delta_R;
        }
      }
    }
  return delta_r;
}


void WZAnalyzer::collectMCsummary(Handle <reco::CandidateCollection> mccands) {
 
  using namespace reco;

  //collections

  MCleptZ1_pdgid = -1;
  MCleptZ2_pdgid = -1;
  MCleptW_pdgid = -1;

  MCleptZ1_pt = -1;
  MCleptZ2_pt = -1;
  MCleptW_pt = -1;
  
  MCleptZ1_eta = -1;
  MCleptZ2_eta = -1;
  MCleptW_eta = -1;
  
  MCleptZ1_phi = -1;
  MCleptZ2_phi = -1;
  MCleptW_phi = -1;
   
  vector<reco::GenParticleCandidate*> Tau;
  vector<reco::GenParticleCandidate*> StableMuons;
  vector<reco::GenParticleCandidate*> StableElectrons;

  
  for(CandidateCollection::const_iterator p = mccands->begin();
      p != mccands->end(); p++) {

    //  reco::Candidate* p_tmp = &(*p);
    reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
  
        if ( (ptr)->status() == 1
                /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
                if ( abs((ptr)->pdgId()) == 11 ) {

                        StableElectrons.push_back((ptr));
                        //cout << "electron MCT\n";
                }
                else if ( abs((ptr)->pdgId()) == 13 ) {

                        StableMuons.push_back((ptr)) ;
                        //cout << "muon MCT\n";
                }
        }
        else if ((ptr)->status() == 2) {
                if ( abs((ptr)->pdgId()) == 15 ) {//tau
                        Tau.push_back((ptr));
                }
        }
  }

//   std::cout << "# Electrons : " << StableElectrons.size()
//          << "# Muons     : " << StableMuons.size() << std::endl
//          << "# Tau     : " << Tau.size() << std::endl;


  vector<reco::GenParticleCandidate*> * StableLeptons[3] = {&StableElectrons, &StableMuons, &Tau};

  vector<reco::GenParticleCandidate*> TauChildLeptons;

  
   //erase electrons and muons from tau from GenparticleCandidate Vector and put to Tau Vector
   for (int i=1; i>=0; i--) {
       for (   vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
               lepton != StableLeptons[i]->end();lepton++) {

                reco::GenParticleCandidate* mcParticleRef = *lepton;

               int parentItr=0;
                 while ( mcParticleRef->mother()!=0 && parentItr<2) {

                         parentItr++;

                         mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
                         if (mcParticleRef==0) break;

                         //if tau
                         if (abs((mcParticleRef)->pdgId())==15 ) {
                               //put into collection
                               TauChildLeptons.push_back(*lepton);
                               StableLeptons[i]->erase(lepton);
                               lepton--;
                               break;
                         }
               }
       }
   }


  bool firstZlept = true;
//   int MC_tauDecayTypeZ1 = 0;
//   int MC_tauDecayTypeZ2 = 0;
//   int MC_tauDecayTypeW  = 0;
  MC_tauDecayTypeZ1 = 0;
  MC_tauDecayTypeZ2 = 0;
  MC_tauDecayTypeW  = 0;


  for (int i=2; i>=0; i--) {
        while (StableLeptons[i]->size() > 0) {
                float maxPt = 0;
                vector<reco::GenParticleCandidate*>::iterator  index = StableLeptons[i]->end();

                for (   vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
                        lepton != StableLeptons[i]->end(); lepton++ ) {

                        if ((*lepton)->pt() > maxPt) {
                                maxPt = (*lepton)->pt();
                                index = lepton;
                        }
                }
                bool Zchild = false;
                bool Wchild = false;
                bool isTau = false;

        
                reco::GenParticleCandidate* mcParticleRef = *index;

                if (abs((*index)->pdgId()) == 15) isTau = true;

                //find W/Z mother
                int parentItr=0;
                while ( mcParticleRef->mother()!=0 && parentItr<2) {

                        if (parentItr>=2) break;
                        parentItr++;
                        mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
                        if (mcParticleRef==0) break;

                        if (abs((mcParticleRef)->pdgId())==23  ) {
                                Zchild=true;
                                break;
                        }
                        if (abs((mcParticleRef)->pdgId())==24  ) {
                                Wchild=true;
                                break;
                        }
                }
        
                
                if (maxPt >0) {
                        int * MCtauDecayTypePtr = 0;
        
                        if (Wchild) {
                                MCleptW_pdgid=(*index)->pdgId();
                                MCleptW_pt=(float)(*index)->pt();
                                MCleptW_eta=(float)(*index)->eta();
                                MCleptW_phi=(float)(*index)->phi();
                               if (isTau) {
                                       MCtauDecayTypePtr = &MC_tauDecayTypeW;
                                       MC_tauDecayTypeW =3;//default to hadronic decay
                               }

                        }
                        if (Zchild) {
                                if (firstZlept) {
                                        firstZlept=false;
                                        MCleptZ1_pdgid=(*index)->pdgId();
                                        MCleptZ1_pt=(float)(*index)->pt();
                                        MCleptZ1_eta=(float)(*index)->eta();
                                        MCleptZ1_phi=(float)(*index)->phi();
                                        if (isTau) {
                                                MCtauDecayTypePtr = &MC_tauDecayTypeZ1;
                                                MC_tauDecayTypeZ1 =3;
                                        }

                                }
                                else  {
                                        MCleptZ2_pdgid=(*index)->pdgId();
                                        MCleptZ2_pt=(float)(*index)->pt();
                                        MCleptZ2_eta=(float)(*index)->eta();
                                        MCleptZ2_phi=(float)(*index)->phi();
                                       if (isTau) {
                                               MCtauDecayTypePtr = &MC_tauDecayTypeZ2;
                                               MC_tauDecayTypeZ2 =3;
                                       }

                                }
                        }
                       //check type of tau decay
                       if (MCtauDecayTypePtr) {
                               for( vector<reco::GenParticleCandidate*>::iterator p = TauChildLeptons.begin();p != TauChildLeptons.end(); p++) {
                                       int pitr = 0;
                                       reco::GenParticleCandidate* mcParticleRef = *p;
                                       while (mcParticleRef->mother() && pitr<2) {
                                               pitr++;
                                               mcParticleRef =(reco::GenParticleCandidate*) mcParticleRef->mother();
                                               if (mcParticleRef==0) break;
                                               if (mcParticleRef == *index) {
                                                       if (abs((*p)->pdgId()) == 11) *MCtauDecayTypePtr = 1;
                                                       if (abs((*p)->pdgId()) == 13) *MCtauDecayTypePtr = 2;
                                               }
                                       }
                                }
                        }
                }
                StableLeptons[i]->erase(index);
        }
  }
}
//define this as a plug-in
// DEFINE_FWK_MODULE(WZAnalyzer);

Revision:	1.56
Committed:	Wed Apr 23 14:23:27 2008 UTC (17 years ago) by beaucero
Content type:	text/plain
Branch:	MAIN
CVS Tags:	SB_23Apr08
Changes since 1.55:	+42 -156 lines
Log Message:	Add Transverse Mass Calculation + DeltaPhi