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.28 2008/01/15 12:45:45 senka 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 "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 <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");
  theGoodElectronCollection    = iConfig.getParameter<edm::InputTag>("GoodElectrons");
  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");
  
  storeHLTresults=false;
  storeHLTresults              = iConfig.getParameter<bool>("storeHLTresults");

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

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


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
   ///

   // 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;


//    Tau.clear();
//    StableMuons.clear();
//    StableElectrons.clear();


   ////////////////////////////////////////
   // 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);

   // good electrons
   Handle<CandidateCollection> goodElectronCands;
   iEvent.getByLabel( theGoodElectronCollection , goodElectronCands);

   // 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;
   //CandidateCollection> jetCands;
   iEvent.getByLabel( theJetCollection , jetCands);


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

   Handle< double> weightHandle;
   if (getEventWeight) {
     iEvent.getByLabel (theWeightCollection, weightHandle);
     if (weightHandle.isValid()) {
       eventWeight = * weightHandle;
     }
   }


    // 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;
        triggerNames= HLTR->getTriggerNames();
        numTriggers = triggerNames.size();
      }

    }   


    // MET corrected (P_muons, P_muons_calo_deposit)

     reco::Particle::LorentzVector MET=computeMET(iEvent, iSetup, looseMuonCands); 
     MET_energy=MET.energy();
     MET_pt=MET.pt();
     MET_eta=MET.eta();
   

   // selected muons

   cout << "\t # loose mu : " <<  looseMuonCands->size() 
        << "\t # tight mu : " <<  goodMuonCands->size() 
        << "\t # loose e : " <<   looseElectronCands->size() 
        << "\t # tight e : " <<   goodElectronCands->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();


   //
   // 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))) ) {
      cout << "Reject Z candidate: built with overlapping leptons -  flavour : "
           << z1->daughter(0)->pdgId() << endl;
      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)) {
          cout << "Overlap between two Z of flavour "  
               << z1->daughter(0)->pdgId() << "\t"
               << z2->daughter(0)->pdgId() << "\t"
               << endl;
        } else {
          cout << "NON OVERLAPPING Zs " << endl;
          n++;
        }
      }
    }
  }


  zFlavour = 0;
  wlFlavour = 0;

  dPhi_WlZ_reco=-15.;    
  dPhi_WZ_reco=-15.;


  //---------------
//   cout<<"??????????????????? before MatchedGenParticle(&(*looseMuonCands)[0]);" <<endl;
//   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) ) {
//        cout << "found overlap: pt = " << l->pt()
//             << "\t eta = " << l->eta() 
//             << "\t phi = " << l->phi() 
//             << "\t compared to : pt " << theZCandidate->daughter(i)->pt()
//             << "\t eta = " << theZCandidate->daughter(i)->eta() 
//             << "\t phi = " << theZCandidate->daughter(i)->phi() 
//             << endl;
          noverlaps++;
        }
      }
//       cout << "Z Flavour : " << zFlavour 
//         << "\t # overlaps for daughter " << i << " : " << noverlaps 
//         << "\t (out of " << looseLeptonCands->size() << " possible) "
//         << endl;
    }

    matching(iEvent,looseMuonCands, 13);
    matching(iEvent,looseElectronCands, 11);
    
    if (zFlavour == 11) {
      //      cout << " filling Z lepton electron properties \n";
      //      if (theZCandidate->daughter(0)) cout << " found Z daughter \n";
      //--------------   
      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) {
      //--------------   
      //      cout << " filling Z lepton muon properties \n"; 
      //      if (theZCandidate->daughter(0)) cout << "found Z daughter \n";
      const reco::Candidate * mcpart = MatchedGenParticle(theZCandidate->daughter(0));
      if (mcpart) {
        //      cout << "GenPt = genParticle->pt()="<< mcpart->pt() <<"\n";     
      } else {
        cout << "NO MATCHED GEN PARTICLE \n";
      }
      //    cout << "dR="<< dR(theZCandidate->daughter(0))<<"\n";
      //--------------   
      leptonProperties["Zmu1"]->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 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();

      cout << "found W lepton \n";

      if (abs(wlFlavour) == 11) {
        leptonProperties["Wel"]->FillProperties(theWlepton, iEvent, iSetup, MatchedGenParticle(theWlepton),dR(theWlepton));
      } else if (abs(wlFlavour) == 13) {
        //      cout << "found W lepton: it is a muon \n";
        leptonProperties["Wmu"]->FillProperties(theWlepton, iEvent, iSetup, MatchedGenParticle(theWlepton),dR(theWlepton));
      }

      dPhi_WlZ_reco=acos(cos(theWlepton->phi()-theZCandidate->phi()));    
      dPhi_WZ_reco=acos(cos((theWlepton->p4()+MET).phi()-theZCandidate->phi()));

//       const reco::Muon * muon =      dynamic_cast<const reco::Muon *>(&(*theWlepton));
//       if(muon != 0)  { /* it's a valid muon */ 
//              cout << "lepton is a muon \n" << endl;
//       }
//       const reco::PixelMatchGsfElectron * electron = dynamic_cast<const reco::PixelMatchGsfElectron *>(&(*theWlepton));
//       if (electron != 0) {
//              cout << "lepton is an electron \n" << endl;
//       }
    }

  }
  
  //  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;

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

  //    cout << "trigName" << triggerNames[i] << "wasRUN:" << HLTR->wasrun(i) << "error " << HLTR->error(i)
  //            << " accept" << HLTR->accept(i) << "\n";

        if (!HLTR->wasrun(i)) {
                //cout << "trigger not run? "<< triggerNames[i] << "\n";
                continue;
        }
        if (HLTR->error(i) ) {
                //cout << "error trigger: "<< triggerNames[i] <<"\n";
                continue;
        }
        if (HLTR->accept(i)) {
                //if (triggerNames[i].compare("mcpath")!=0) 
                //TRaccept=true;
        }else {
                continue;
        }

        if (triggerNames[i].compare("HLT1Electron")==0) triggerBitmask |= 1;
        if (triggerNames[i].compare("HLT1ElectronRelaxed")==0) triggerBitmask |= 2;
        if (triggerNames[i].compare("HLT2Electron")==0) triggerBitmask |= 4;
        if (triggerNames[i].compare("HLT2ElectronRelaxed")==0) triggerBitmask |= 8;

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

        if (triggerNames[i].compare("HLTXElectronMuon")==0) triggerBitmask |= 256;
        if (triggerNames[i].compare("HLTXElectronMuonRelaxed")==0) triggerBitmask |= 512;
     }
   }
  

  // 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();
  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->daughter(ct)->pt()>0 
                    && (abs(theZCandidate->daughter(ct)->pdgId())==13 
                        || abs(theZCandidate->daughter(ct)->pdgId())==11) ) { 

                        dPhi[ct]=fabs(theZCandidate->daughter(ct)->phi() - jet->phi());
                        dEta[ct]=fabs(theZCandidate->daughter(ct)->eta() - jet->eta());
                }
        }

        if (theWlepton->pt()>0 && (abs(theWlepton->pdgId())==13 || abs(theWlepton->pdgId())==11 )) {
                dPhi[2]=fabs(theWlepton->phi() - jet->phi());
                dEta[2]=fabs(theWlepton->eta() - jet->eta());
        }

        for (int k=0;k<3;k++) {
                if (matchedLept[k] || dPhi[k]<0.0 || dEta[k]<0.0) continue;

                while (dPhi[k]>6.28318531) dPhi[k]-= 6.28318531;
                if (dPhi[k]>3.14159265) dPhi[k]= 6.28318531-dPhi[k];

                float deltaR=sqrt(dPhi[k]*dPhi[k]+dEta[k]*dEta[k]);
                if (deltaR<0.15) {
                        skipJet=true;
                        matchedLept[k]=true;
                        break;
                }
        }
        if (skipJet) continue;
   }

    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]);

  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;
        }
        leptonProperties["electron"]->FillProperties(&(*el), 
                                                     iEvent, iSetup, 
                                                     MatchedGenParticle(&(*el)),dR(&(*el)));
        electronTree->Fill();

      }
    }
  }


}


// ------------ 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");

  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["Wel"]->CreateBranches(wzTree, "WEl_");
  leptonProperties["ZEl1"]->CreateBranches(wzTree, "ZEl1_");
  leptonProperties["ZEl2"]->CreateBranches(wzTree, "ZEl2_");

  leptonProperties["Wmu"] ->CreateBranches(wzTree, "Wmu_");
  leptonProperties["Zmu1"] ->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");

  }

  // MET branch
  wzTree->Branch("MET_energy",&MET_energy,"MET_energy/F");
  wzTree->Branch("MET_pt",&MET_pt,"MET_pt/F");
  wzTree->Branch("MET_eta",&MET_eta,"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");

  // # 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");   

}

// ------------ 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");

  // 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("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");

  // MET properties

  // Jet properties
  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");

  //
  electronTree->Branch("used",   &electronUse, "used/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;
  
  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++) {
//     cout << "MET:  E = " << met->energy()
//       << "\t Pt = " << met->pt()
//       <<endl;
  }

      // sum of p of muons
    reco::Particle::LorentzVector Muon_p(0,0,0,0);   // sum of muons 4-vectors
    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
    double px=0;
    double py=0;
    double pz=0;    
    double Et=0;

    for (unsigned int i=0; i<MuonCands->size();i++){
      Muon_p=Muon_p+(*MuonCands)[i].p4();
      px=(*MuonCands)[i].px();
      py=(*MuonCands)[i].py();
      pz=(*MuonCands)[i].pz();
      Et=(*MuonCands)[i].et();
      Muon_p_4v=Muon_p_4v+reco::Particle::LorentzVector(px,py,0,Et);
    }
    //    cout<< "####### sum Muon_p="<< Muon_p<<endl;
    //    cout<< "####### sum Muon_p correct="<< Muon_p_4v<<endl;
      // sum of p of muons left in calorimeters
    reco::Particle::LorentzVector Muon_p_cal(0,0,0,0); // sum of muons depositions 4-vectors
    reco::Particle::LorentzVector Muon_p_cal_4v(0,0,0,0);  // sum of muons depositions 4-vectors in (px,py,0,Et), since MET looks like that
    edm::Handle<reco::MuIsoDepositAssociationMap> ecalIso;
    iEvent.getByLabel("muGlobalIsoDepositCalByAssociatorTowers","ecal", ecalIso);
    edm::Handle<reco::MuIsoDepositAssociationMap> hcalIso;
    iEvent.getByLabel("muGlobalIsoDepositCalByAssociatorTowers","hcal", hcalIso);
    edm::Handle<reco::MuIsoDepositAssociationMap> hocalIso;
    iEvent.getByLabel("muGlobalIsoDepositCalByAssociatorTowers","ho", hocalIso);

    float Et_ecal; 
    float Et_hcal; 
    float Et_ho; 
    double eta_ecal; 
    double eta_hcal; 
    double eta_ho; 
    double phi_ecal; 
    double phi_hcal; 
    double phi_ho; 
    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;

    for (unsigned int i=0; i<MuonCands->size();i++){
      Et_ecal=-2.;  
      Et_hcal=-2.;    
      Et_ho=-2.;  
      eta_ecal=-2.;  
      eta_hcal=-2.;    
      eta_ho=-2.;  
      phi_ecal=-2.;  
      phi_hcal=-2.;    
      phi_ho=-2.;  
      const reco::Muon* muons=dynamic_cast<const reco::Muon *> (&(*MuonCands)[i]);
      const reco::MuIsoDeposit& depEcal     = (*ecalIso)[muons->combinedMuon()];
      const reco::MuIsoDeposit& depHcal     = (*hcalIso)[muons->combinedMuon()];
      const reco::MuIsoDeposit& depHOcal   = (*hocalIso)[muons->combinedMuon()];
      Et_ecal=depEcal.muonEnergy();  
      Et_hcal=depHcal.muonEnergy();    
      Et_ho=depHOcal.muonEnergy();  
      eta_ecal=depEcal.eta();  
      eta_hcal=depHcal.eta();    
      eta_ho=depHOcal.eta(); 
      phi_ecal=depEcal.phi();  
      phi_hcal=depHcal.phi();    
      phi_ho=depHOcal.phi(); 
      tg_theta_pola_ecal=TMath::Exp(-eta_ecal);
      tg_theta_pola_hcal=TMath::Exp(-eta_hcal);      
      tg_theta_pola_ho=TMath::Exp(-eta_ho);
      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);
      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_ecal);
      px_hcal=pt_hcal*TMath::Cos(phi_hcal);
      px_ho=pt_ho*TMath::Cos(phi_ho);
      py_ecal=pt_ecal*TMath::Sin(phi_ecal);
      py_hcal=pt_hcal*TMath::Sin(phi_hcal);
      py_ho=pt_ho*TMath::Sin(phi_ho);

      Muon_p_cal=Muon_p_cal + reco::Particle::LorentzVector(px_ecal+px_hcal+px_ho,py_ecal+py_hcal+py_ho,pz_ecal+pz_hcal+pz_ho,E_ecal+E_hcal+E_ho);
      Muon_p_cal_4v=Muon_p_cal_4v + reco::Particle::LorentzVector(px_ecal+px_hcal+px_ho,py_ecal+py_hcal+py_ho,0,Et_ecal+Et_hcal+Et_ho);
    }

    reco::Particle::LorentzVector MET=mets->begin()->p4()-Muon_p_4v-Muon_p_cal_4v;  // substract p from muons and p of calo muon deposit from MET 
    //    cout <<" MET no corrections"<<mets->begin()->p4() << endl;
    //    cout <<" MET with corrections"<< MET << endl;

  return MET;
  
}

////////////////////////
// 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};

  bool firstZlept = true;
  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;
                        }
                }
                int parentItr=0;
                bool Zchild = false;
                bool Wchild = false;
                bool Tauchild = false;
                
                reco::GenParticleCandidate* mcParticleRef = *index;

                //find W/Z mother
                while ( mcParticleRef->mother()!=0) {
                
                        if (parentItr>=2) break;
                        parentItr++;
                        mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
                        if (mcParticleRef->mother()==0) break;
                        
                        //if tau
                        if (abs((mcParticleRef)->pdgId())==15  ) {
                                Tauchild=true;
                                break;
                        }

                        if (abs((mcParticleRef)->pdgId())==23  ) {
                                Zchild=true;
                                break;
                        }
                        if (abs((mcParticleRef)->pdgId())==24  ) {
                                Wchild=true;
                                break;
                        }
                }
        
                if (maxPt >0) {
                        if (Tauchild) {                 
                                parentItr = 0;
                                while (mcParticleRef->mother()!=0) {
                                        if (parentItr>=2) break;
                                        parentItr++;
                                        mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
                                        if (mcParticleRef->mother()==0) break;

                                        if (abs((mcParticleRef)->pdgId())==23  ) {
                                                if (MCleptZ1_pt == float(mcParticleRef->pt())
                                                && MCleptZ1_eta == float(mcParticleRef->eta())
                                                && MCleptZ1_phi == float(mcParticleRef->phi()) )
                                                        MC_tauDecayTypeZ1 = i;
                                                else MC_tauDecayTypeZ2 = i;
                                                break;
                                        }
                                        if (abs((mcParticleRef)->pdgId())==24  ) {
                                                MC_tauDecayTypeW = i;
                                                break;
                                        }
                                }
                        }
                
                        if (Wchild) {
                                MCleptW_pdgid=(*index)->pdgId();
                                MCleptW_pt=(float)(*index)->pt();
                                MCleptW_eta=(float)(*index)->eta();
                                MCleptW_phi=(float)(*index)->phi();
                                if (abs(MCleptW_pdgid)==15) 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 (abs(MCleptZ1_pdgid)==15) MC_tauDecayTypeZ1 =3;
                                }
                                else  {
                                        MCleptZ2_pdgid=(*index)->pdgId();
                                        MCleptZ2_pt=(float)(*index)->pt();
                                        MCleptZ2_eta=(float)(*index)->eta();
                                        MCleptZ2_phi=(float)(*index)->phi();
                                        if (abs(MCleptZ2_pdgid)==15) MC_tauDecayTypeZ2 =3;
                                }
                        }       
                        StableLeptons[i]->erase(index);
                }
        }
  }
//   std::cout << "# Electrons end: " << StableElectrons.size()
//          << "# Muons     end: " << StableMuons.size() << std::endl
//          << "# Tau     end: " << Tau.size() << std::endl;

}
//define this as a plug-in
// DEFINE_FWK_MODULE(WZAnalyzer);
Revision:	1.29
Committed:	Mon Jan 21 13:29:26 2008 UTC (17 years, 3 months ago) by smorovic
Content type:	text/plain
Branch:	MAIN
Changes since 1.28:	+43 -2 lines
Log Message:	add jet rejection if DeltaR(jet,WZlepton)<0.15