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.44 2008/02/13 18:49:38 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");
  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;
}


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;


//    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;
   Handle<int> genProcessID;
   Handle<int> alpgenProcessID;
   Handle<double> genEventScale;

   //cout << "get weight "<<getEventWeight << endl;
   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();
   //   cout << " MET_energy(Et)="<< MET_energy<< " MET_pt="<< MET_pt<<" MET_Ex="<< MET_energy*TMath::Cos(MET.phi())<<" MET_Ey="<< MET_energy*TMath::Sin(MET.phi()) <<" MET_px="<< MET_px <<" MET_py="<< MET_py <<endl;
   
   
   // 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;
  zMass = -10;
  zPt  = -10;
  zEta = -10;
  zPhi = -10;

  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 (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;
     }
   }
  
  // 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]=fabs(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]=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(jet->pt() > 10 && fabs(jet->eta())<2.5 )NbJetsWithLeptons++;
        if (skipJet) 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)));
      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_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("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");
  wzTree->Branch("processID", &processID, "processID/I");
  wzTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
  wzTree->Branch("genEventScale",    &eventScale,  "genEventScale/I");
  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("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("NbJets",     &NbJets,  "NbJets/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");
  //
  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/I");
  electronTree->Branch("ElectronTreesInEvent", &numberOfElectronTrees, "ElectronTreesInEvent/I");

  electronTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
  electronTree->Branch("EventID", &EventID, "EventID/I");
  electronTree->Branch("NbJetsWithLeptons",     &NbJetsWithLeptons,  "NbJetsWithLeptons/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;
  
}

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


  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);
        }
  }
//   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.45
Committed:	Tue Feb 19 11:50:47 2008 UTC (17 years, 2 months ago) by beaucero
Content type:	text/plain
Branch:	MAIN
CVS Tags:	sm_21_Feb_08, SB_19Feb08
Changes since 1.44:	+3 -3 lines
Log Message:	Comment cout for MET