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.7 2007/11/23 18:01:55 vuko 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/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
#include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectron.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

  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");
  
  storeHLTresults=false;
  storeHLTresults              = iConfig.getParameter<bool>("storeHLTresults");

  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;


#ifdef THIS_IS_AN_EVENT_EXAMPLE
   Handle<ExampleData> pIn;
   iEvent.getByLabel("example",pIn);
#endif
   
#ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
   ESHandle<SetupData> pSetup;
   iSetup.get<SetupRecord>().get(pSetup);
#endif

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

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

   
    // get hold of TriggerResults
    Handle<TriggerResults> HLTR;
    if (storeHLTresults) {
      iEvent.getByType(HLTR);
      if (firstTriggerResult) {
        firstTriggerResult = false;
        triggerNames= HLTR->getTriggerNames();
        numTriggers = triggerNames.size();
      }

    }   

   // Reset a few things

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

   // selected muons


   // 

   int nzee   = zeeCands->size();
   int nzmumu = zmumuCands->size();

   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;

   Handle<CandidateCollection>  zCands[2] = {zeeCands,zmumuCands};

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


  float dzmass_min = 100.;

  for(CandidateCollection::const_iterator z1 = zllCands->begin(); 
      z1 != zllCands->end(); ++ z1 ) {
    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;
        }
      }
    }
  }


  zFlavour = 0;
  wlFlavour = 0;

  /// Find W lepton


  int nwl_candidates=0;
  if (theZCandidate) {


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

    if (zFlavour == 11) {
      leptonProperties["ZEl1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup);
      leptonProperties["ZEl2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup);
    } else if (zFlavour == 13) {
      leptonProperties["Zmu1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup);
      leptonProperties["Zmu2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup);
    }

    float max_pt = 0.;


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

      int id = lepton->pdgId();
      cout << "Tight lepton: " << id << endl;

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


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


      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)  { /* it's a valid electron */ 
        cout << "lepton is an electron \n" << endl;
        elWTree->Fill();
      }
    }

  }
  
  Handle <reco::GenParticleCandidateCollection> mccands;
  iEvent.getByLabel( "genParticleCandidateProducer",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;
     }
   }
  
  wzTree->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" ) ;

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

  elWTree = new TTree("WElTree","W electron informations");

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

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

  }
  
  
}

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

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

}


void WZAnalyzer::initialiseTree() {

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

void WZAnalyzer::collectMCsummary(Handle <reco::GenParticleCandidateCollection> mccands) {

  //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;
  
  reco::Candidate::iterator p;
  for (p = ((reco::Candidate*)&mccands)->begin(); p != ((reco::Candidate*)&mccands)->end(); ++p ) {

        reco::Candidate* p_tmp = &(*p);
        reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) p_tmp ;
  
        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;
                        }

                //      cout << "mother: " << (mcParticleRef)->pdg_id()
                //      <<" pt="<< c_p4.pt() << " eta="<<c_p4.eta()
                //      <<" status=" <<mcParticleRef->status() <<"\n";

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

}
//define this as a plug-in
// DEFINE_FWK_MODULE(WZAnalyzer);
Revision:	1.8
Committed:	Mon Nov 26 16:10:55 2007 UTC (17 years, 5 months ago) by smorovic
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+190 -7 lines
Log Message:	monte carlo analysis
#	User	Rev	Content
1	vuko	1.1	// -- C++ --
2			//
3			// Package: WZAnalyzer
4			// Class: WZAnalyzer
5			//
6	vuko	1.2	/**\class WZAnalyzer WZAnalyzer.cc Vuko/WZAnalysis/src/WZAnalyzer.cc
7	vuko	1.1
8			Description: <one line class summary>
9
10			Implementation:
11			<Notes on implementation>
12			*/
13			//
14			// Original Author: Vuko Brigljevic
15			// Created: Fri Nov 9 11:07:27 CET 2007
16	smorovic	1.8	// $Id: WZAnalyzer.cc,v 1.7 2007/11/23 18:01:55 vuko Exp $
17	vuko	1.1	//
18			//
19
20
21			// system include files
22			#include <memory>
23
24			// user include files
25			#include "FWCore/Framework/interface/Frameworkfwd.h"
26			#include "FWCore/Framework/interface/EDAnalyzer.h"
27
28			#include "FWCore/Framework/interface/Event.h"
29			#include "FWCore/Framework/interface/MakerMacros.h"
30
31			#include "FWCore/ParameterSet/interface/ParameterSet.h"
32
33	vuko	1.2	#include "Vuko/WZAnalysis/interface/WZAnalyzer.h"
34	vuko	1.1	#include "DataFormats/Candidate/interface/OverlapChecker.h"
35
36	vuko	1.2	#include "Vuko/WZAnalysis/interface/ElectronProperties.h"
37	vuko	1.3	#include "Vuko/WZAnalysis/interface/MuonProperties.h"
38	vuko	1.1
39	smorovic	1.6	#include "DataFormats/Common/interface/TriggerResults.h"
40	vuko	1.1
41			//--- muon AOD:
42			#include "DataFormats/JetReco/interface/CaloJetCollection.h"
43			#include "DataFormats/EgammaCandidates/interface/Electron.h"
44			#include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
45			#include "DataFormats/MuonReco/interface/MuonFwd.h"
46			#include "DataFormats/MuonReco/interface/Muon.h"
47			#include "DataFormats/EgammaReco/interface/BasicCluster.h"
48			#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
49			#include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectron.h"
50
51
52			#include "TFile.h"
53			#include "TH1F.h"
54			#include "TH2F.h"
55			#include "TTree.h"
56
57			#include <map>
58
59			//
60			// constants, enums and typedefs
61			//
62
63			//
64			// static data member definitions
65			//
66
67			//
68			// constructors and destructor
69			//
70			WZAnalyzer::WZAnalyzer(const edm::ParameterSet& iConfig)
71
72			{
73			//now do what ever initialization is needed
74
75			theLooseMuonCollection = iConfig.getParameter<edm::InputTag>("LooseMuons");
76			theGoodMuonCollection = iConfig.getParameter<edm::InputTag>("GoodMuons");
77			theLooseElectronCollection = iConfig.getParameter<edm::InputTag>("LooseElectrons");
78			theGoodElectronCollection = iConfig.getParameter<edm::InputTag>("GoodElectrons");
79			theTightLeptonCollection = iConfig.getParameter<edm::InputTag>("TightLeptons");
80			// theMediumElectronCollection = iConfig.getParameter<edm::InputTag>("LooseElectrons");
81			// theTightElectronCollection = iConfig.getParameter<edm::InputTag>("TightElectrons");
82			// Z's
83			theZeeCollection = iConfig.getParameter<edm::InputTag>("ZtoEE");
84			theZmumuCollection = iConfig.getParameter<edm::InputTag>("ZtoMuMu");
85	vuko	1.3	theZllCollection = iConfig.getParameter<edm::InputTag>("ZtoLL");
86	smorovic	1.6
87			storeHLTresults=false;
88			storeHLTresults = iConfig.getParameter<bool>("storeHLTresults");
89	vuko	1.1
90	smorovic	1.6	firstTriggerResult = true;
91	vuko	1.1	}
92
93
94			WZAnalyzer::~WZAnalyzer()
95			{
96
97			// do anything here that needs to be done at desctruction time
98			// (e.g. close files, deallocate resources etc.)
99
100			}
101
102
103			//
104			// member functions
105			//
106
107			// ------------ method called to for each event ------------
108			void
109			WZAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
110			{
111			using namespace edm;
112			using namespace reco;
113			using namespace std;
114
115
116			#ifdef THIS_IS_AN_EVENT_EXAMPLE
117			Handle<ExampleData> pIn;
118			iEvent.getByLabel("example",pIn);
119			#endif
120
121			#ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
122			ESHandle<SetupData> pSetup;
123			iSetup.get<SetupRecord>().get(pSetup);
124			#endif
125
126			const Candidate * theZCandidate = 0;
127			const Candidate * theWlepton = 0;
128
129			// Z->mumu
130			Handle<CandidateCollection> zmumuCands;
131			iEvent.getByLabel( theZmumuCollection , zmumuCands);
132
133
134			// Z->ee
135			Handle<CandidateCollection> zeeCands;
136			iEvent.getByLabel( theZeeCollection , zeeCands);
137
138	vuko	1.3	// Z->ee
139			Handle<CandidateCollection> zllCands;
140			iEvent.getByLabel( theZllCollection , zllCands);
141
142	vuko	1.1	// loose electrons
143			Handle<CandidateCollection> looseElectronCands;
144			iEvent.getByLabel( theLooseElectronCollection , looseElectronCands);
145
146			// good electrons
147			Handle<CandidateCollection> goodElectronCands;
148			iEvent.getByLabel( theGoodElectronCollection , goodElectronCands);
149
150			// loose muons
151			Handle<CandidateCollection> looseMuonCands;
152			iEvent.getByLabel( theLooseMuonCollection , looseMuonCands);
153
154			// good muons
155			Handle<CandidateCollection> goodMuonCands;
156			iEvent.getByLabel( theGoodMuonCollection , goodMuonCands);
157
158
159			// tight leptons
160			Handle<CandidateCollection> tightLeptonCands;
161			iEvent.getByLabel( theTightLeptonCollection , tightLeptonCands);
162
163	smorovic	1.6
164			// get hold of TriggerResults
165			Handle<TriggerResults> HLTR;
166			if (storeHLTresults) {
167	vuko	1.7	iEvent.getByType(HLTR);
168			if (firstTriggerResult) {
169			firstTriggerResult = false;
170			triggerNames= HLTR->getTriggerNames();
171			numTriggers = triggerNames.size();
172			}
173	smorovic	1.6
174	vuko	1.7	}
175	vuko	1.1
176			// Reset a few things
177
178	vuko	1.3	for (map<string,wzana::LeptonProperties* >::iterator iter = leptonProperties.begin();
179			iter!=leptonProperties.end(); iter++)
180	vuko	1.1	{
181			iter->second->ResetValues();
182			}
183
184			// selected muons
185
186
187			//
188
189			int nzee = zeeCands->size();
190			int nzmumu = zmumuCands->size();
191
192			cout << "\t # loose mu : " << looseMuonCands->size()
193			<< "\t # tight mu : " << goodMuonCands->size()
194			<< "\t # loose e : " << looseElectronCands->size()
195			<< "\t # tight e : " << goodElectronCands->size()
196			<< "\t # tight l : " << tightLeptonCands->size()
197			<< "\t # Z->mumu : " << zmumuCands->size()
198			<< "\t # Z->ee : " << zeeCands->size()
199	vuko	1.7	<< "\t # Z->ll : " << zllCands->size()
200	vuko	1.1	<< endl;
201
202			Handle<CandidateCollection> zCands[2] = {zeeCands,zmumuCands};
203
204			//
205	vuko	1.7	// Find best Z (in Z->ll collection (merged Z->ee & Z->mumu)
206	vuko	1.1	//
207			::OverlapChecker overlap;
208
209	vuko	1.3
210			float dzmass_min = 100.;
211
212			for(CandidateCollection::const_iterator z1 = zllCands->begin();
213			z1 != zllCands->end(); ++ z1 ) {
214			float dzmass = fabs( z1->mass() - 91.188);
215			if (dzmass < dzmass_min) {
216	vuko	1.7	theZCandidate = &(*z1);
217	vuko	1.3	dzmass_min = dzmass;
218			}
219			//
220			//Get back Electrons from Z and W
221			for(CandidateCollection::const_iterator z2 = z1;
222			z2 != zllCands->end(); ++ z2 ) {
223			if (z1 != z2) {
224			if (overlap(z1,z2)) {
225	vuko	1.7	cout << "Overlap between two Z of flavour "
226			<< z1->daughter(0)->pdgId() << "\t"
227			<< z2->daughter(0)->pdgId() << "\t"
228			<< endl;
229	vuko	1.3	} else {
230			cout << "NON OVERLAPPING Zs " << endl;
231			}
232			}
233			}
234			}
235
236
237	vuko	1.1	zFlavour = 0;
238			wlFlavour = 0;
239
240	vuko	1.4	/// Find W lepton
241
242	vuko	1.3
243	vuko	1.1	int nwl_candidates=0;
244			if (theZCandidate) {
245
246
247			zMass = theZCandidate->mass();
248			zFlavour = abs(theZCandidate->daughter(0)->pdgId());
249			zPt = theZCandidate->pt();
250			zEta = theZCandidate->eta();
251			zPhi = theZCandidate->phi();
252
253	vuko	1.4	if (zFlavour == 11) {
254	vuko	1.5	leptonProperties["ZEl1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup);
255			leptonProperties["ZEl2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup);
256	vuko	1.4	} else if (zFlavour == 13) {
257	vuko	1.5	leptonProperties["Zmu1"]->FillProperties(theZCandidate->daughter(0), iEvent, iSetup);
258			leptonProperties["Zmu2"]->FillProperties(theZCandidate->daughter(1), iEvent, iSetup);
259	vuko	1.4	}
260
261	vuko	1.1	float max_pt = 0.;
262
263
264			// Now find lepton that will be associated to W
265			// among leptons not used for the Z reconstruction
266			for(CandidateCollection::const_iterator lepton = tightLeptonCands->begin();
267			lepton != tightLeptonCands->end(); lepton++) {
268
269
270			if ( overlap(theZCandidate, lepton) ) continue; // Ignore if lepton used for the Z
271
272			nwl_candidates++;
273
274			// If more than one candidate: choose leading pt
275			if (lepton->pt() > max_pt) {
276			theWlepton = &(*lepton);
277			max_pt = lepton->pt();
278			}
279
280			int id = lepton->pdgId();
281			cout << "Tight lepton: " << id << endl;
282
283			if(lepton->hasMasterClone()) {
284			cout << "SHOUT: TIGHT LEPTON IS SHALLOW COPY !!! (SHOULD NOT BE!) \n";
285			}
286			//
287			}
288			if (theWlepton) {
289			wlFlavour = theWlepton->pdgId();
290			wlCharge = theWlepton->charge();
291			wlPt = theWlepton->pt();
292			wlEta = theWlepton->eta();
293			wlPhi = theWlepton->phi();
294
295
296	vuko	1.3	if (abs(wlFlavour) == 11) {
297			leptonProperties["Wel"]->FillProperties(theWlepton, iEvent, iSetup);
298			} else if (abs(wlFlavour) == 13) {
299			leptonProperties["Wmu"]->FillProperties(theWlepton, iEvent, iSetup);
300			}
301
302
303	vuko	1.1	const reco::Muon * muon = dynamic_cast<const reco::Muon >(&(theWlepton));
304			if(muon != 0) { /* it's a valid muon */
305			cout << "lepton is a muon \n" << endl;
306			}
307			const reco::PixelMatchGsfElectron * electron = dynamic_cast<const reco::PixelMatchGsfElectron >(&(theWlepton));
308			if(electron != 0) { /* it's a valid electron */
309			cout << "lepton is an electron \n" << endl;
310			elWTree->Fill();
311			}
312			}
313
314			}
315	smorovic	1.6
316	smorovic	1.8	Handle <reco::GenParticleCandidateCollection> mccands;
317			iEvent.getByLabel( "genParticleCandidateProducer",mccands );
318	smorovic	1.6
319	smorovic	1.8	collectMCsummary(mccands);
320	smorovic	1.6
321	smorovic	1.8	if (storeHLTresults) {
322
323	smorovic	1.6	//def: 1:single electron 2:single relaxed 3:double electron 4: double relaxed 5,6,7,8: muon e+mu: 9
324	smorovic	1.8
325	smorovic	1.6	triggerBitmask=0;
326	smorovic	1.8
327	smorovic	1.6	for (size_t i=0; i<numTriggers; ++i) {
328
329			// cout << "trigName" << triggerNames[i] << "wasRUN:" << HLTR->wasrun(i) << "error " << HLTR->error(i)
330			// << " accept" << HLTR->accept(i) << "\n";
331
332			if (!HLTR->wasrun(i)) {
333			//cout << "trigger not run? "<< triggerNames[i] << "\n";
334			continue;
335			}
336			if (HLTR->error(i) ) {
337			//cout << "error trigger: "<< triggerNames[i] <<"\n";
338			continue;
339			}
340			if (HLTR->accept(i)) {
341			//if (triggerNames[i].compare("mcpath")!=0)
342			//TRaccept=true;
343			}else {
344			continue;
345			}
346
347			if (triggerNames[i].compare("HLT1Electron")==0) triggerBitmask \|= 1;
348			if (triggerNames[i].compare("HLT1ElectronRelaxed")==0) triggerBitmask \|= 2;
349			if (triggerNames[i].compare("HLT2Electron")==0) triggerBitmask \|= 4;
350			if (triggerNames[i].compare("HLT2ElectronRelaxed")==0) triggerBitmask \|= 8;
351
352			if (triggerNames[i].compare("HLT1MuonIso")==0) triggerBitmask \|= 16;
353			if (triggerNames[i].compare("HLT1MuonNonIso")==0) triggerBitmask \|= 32;
354			if (triggerNames[i].compare("HLT2MuonNonIso")==0) triggerBitmask \|= 64;
355			if (triggerNames[i].compare("HLT2MuonZ")==0) triggerBitmask \|= 128;
356
357			if (triggerNames[i].compare("HLTXElectronMuon")==0) triggerBitmask \|= 256;
358			if (triggerNames[i].compare("HLTXElectronMuonRelaxed")==0) triggerBitmask \|= 512;
359			}
360			}
361
362	vuko	1.1	wzTree->Fill();
363
364			}
365
366
367			// ------------ method called once each job just before starting event loop ------------
368			void
369			WZAnalyzer::beginJob(const edm::EventSetup&)
370			{
371
372			using namespace wzana;
373
374			theFile = new TFile( "wz.root", "RECREATE" ) ;
375
376			wzTree = new TTree("WZTree","WZ informations");
377
378			elWTree = new TTree("WElTree","W electron informations");
379
380	vuko	1.3	leptonProperties["Wel"] = new ElectronProperties();
381			leptonProperties["ZEl1"] = new ElectronProperties();
382			leptonProperties["ZEl2"] = new ElectronProperties();
383
384			leptonProperties["Wmu"] = new MuonProperties();
385			leptonProperties["Zmu1"] = new MuonProperties();
386			leptonProperties["Zmu2"] = new MuonProperties();
387	vuko	1.1
388	vuko	1.3	leptonProperties["Wel"]->CreateBranches(wzTree, "WEl_");
389			leptonProperties["Wel"]->CreateBranches(elWTree,"ElW_");
390	vuko	1.4	leptonProperties["ZEl1"]->CreateBranches(wzTree, "ZEl1_");
391			leptonProperties["ZEl2"]->CreateBranches(wzTree, "ZEl2_");
392	vuko	1.1
393			initialiseTree();
394	smorovic	1.6
395			//prepare HLT TriggerResults branch
396			if (storeHLTresults) {
397			wzTree->Branch("WZ_hltBitmask",&triggerBitmask,"WZ_hltBitmask/I");
398
399			}
400
401
402	vuko	1.1	}
403
404			// ------------ method called once each job just after ending the event loop ------------
405			void
406			WZAnalyzer::endJob() {
407
408			theFile->Write();
409			theFile->Close();
410
411			}
412
413
414			void WZAnalyzer::initialiseTree() {
415
416			// Z properties
417			wzTree->Branch("Zmass", &zMass, "Zmass/F");
418			wzTree->Branch("ZId", &zFlavour, "Zid/I");
419			wzTree->Branch("ZPt", &zPt, "ZPt/F");
420			wzTree->Branch("ZEta", &zEta, "ZEta/F");
421			wzTree->Branch("ZPhi", &zPhi, "ZPhi/F");
422
423			// W Properties
424			wzTree->Branch("WlId", &wlFlavour, "Wlid/I");
425			wzTree->Branch("WlCharge", &wlCharge, "WlCharge/I");
426			wzTree->Branch("WlPt", &wlPt, "WlPt/F");
427	smorovic	1.8
428			wzTree->Branch("MC_leptonZ1_pt", &MCleptZ1_pt,"MC_leptonZ1_pt/F");
429			wzTree->Branch("MC_leptonZ2_pt", &MCleptZ2_pt,"MC_leptonZ2_pt/F");
430			wzTree->Branch("MC_leptonW_pt", &MCleptW_pt, "MC_leptonW_pt/F");
431
432			wzTree->Branch("MC_leptonZ1_eta", &MCleptZ1_eta,"MC_leptonZ1_eta/F");
433			wzTree->Branch("MC_leptonZ2_eta", &MCleptZ2_eta,"MC_leptonZ2_eta/F");
434			wzTree->Branch("MC_leptonW_eta", &MCleptW_eta, "MC_leptonW_eta/F");
435
436			wzTree->Branch("MC_leptonZ1_phi", &MCleptZ1_phi,"MC_leptonZ1_phi/F");
437			wzTree->Branch("MC_leptonZ2_phi", &MCleptZ2_phi,"MC_leptonZ2_phi/F");
438			wzTree->Branch("MC_leptonW_phi", &MCleptW_phi, "MC_leptonW_phi/F");
439
440			wzTree->Branch("MC_leptonZ1_pdgid", &MCleptZ1_pdgid,"MC_leptonZ1_pdgid/I");
441			wzTree->Branch("MC_leptonZ2_pdgid", &MCleptZ2_pdgid,"MC_leptonZ2_pdgid/I");
442			wzTree->Branch("MC_leptonW_pdgid", &MCleptW_pdgid, "MC_leptonW_pdgid/I");
443
444			wzTree->Branch("MC_TauDecayType_fromZ1", &MC_tauDecayTypeZ1,"MC_TauDecayType_fromZ1/I");
445			wzTree->Branch("MC_TauDecayType_fromZ2", &MC_tauDecayTypeZ2,"MC_TauDecayType_fromZ2/I");
446			wzTree->Branch("MC_TauDecayType_fromW", &MC_tauDecayTypeW, "MC_TauDecayType_fromW/I");
447			}
448	vuko	1.1
449	smorovic	1.8	void WZAnalyzer::collectMCsummary(Handle <reco::GenParticleCandidateCollection> mccands) {
450
451			//collections
452			MCleptZ1_pdgid = -1;
453			MCleptZ2_pdgid = -1;
454			MCleptW_pdgid = -1;
455
456			MCleptZ1_pt = -1;
457			MCleptZ2_pt = -1;
458			MCleptW_pt = -1;
459
460			MCleptZ1_eta = -1;
461			MCleptZ2_eta = -1;
462			MCleptW_eta = -1;
463
464			MCleptZ1_phi = -1;
465			MCleptZ2_phi = -1;
466			MCleptW_phi = -1;
467
468
469
470
471			vector<reco::GenParticleCandidate*> Tau;
472			vector<reco::GenParticleCandidate*> StableMuons;
473			vector<reco::GenParticleCandidate*> StableElectrons;
474
475			reco::Candidate::iterator p;
476			for (p = ((reco::Candidate)&mccands)->begin(); p != ((reco::Candidate)&mccands)->end(); ++p ) {
477	vuko	1.1
478	smorovic	1.8	reco::Candidate* p_tmp = &(*p);
479			reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) p_tmp ;
480
481			if ( (ptr)->status() == 1
482			/&& (ptr)->momentum().perp() > 2./ ) { // stable particles with pt>2 only
483			if ( abs((ptr)->pdgId()) == 11 ) {
484			StableElectrons.push_back((ptr));
485			//cout << "electron MCT\n";
486			}
487			else if ( abs((ptr)->pdgId()) == 13 ) {
488			StableMuons.push_back((ptr)) ;
489			//cout << "muon MCT\n";
490			}
491			}
492			else if ((ptr)->status() == 2) {
493			if ( abs((ptr)->pdgId()) == 15 ) {//tau
494			Tau.push_back((ptr));
495			}
496			}
497			}
498			//std::cout << "# Electrons : " << StableElectrons.size()
499			// << "# Muons : " << StableMuons.size() << std::endl;
500			// << "# Tau : " << Tau.size() << std::endl;
501
502			vector<reco::GenParticleCandidate> StableLeptons[3] = {&StableElectrons, &StableMuons, &Tau};
503
504			bool firstZlept = true;
505			MC_tauDecayTypeZ1 = 0;
506			MC_tauDecayTypeZ2 = 0;
507			MC_tauDecayTypeW = 0;
508
509
510			for (int i=2; i>=0; i--) {
511			while (StableLeptons[i]->size() > 0) {
512			float maxPt = 0;
513			vector<reco::GenParticleCandidate*>::iterator index = StableLeptons[i]->end();
514
515			for ( vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
516			lepton != StableLeptons[i]->end(); lepton++ ) {
517
518			if ((*lepton)->pt() > maxPt) {
519			maxPt = (*lepton)->pt();
520			index = lepton;
521			}
522			}
523			int parentItr=0;
524			bool Zchild = false;
525			bool Wchild = false;
526			bool Tauchild = false;
527
528			reco::GenParticleCandidate* mcParticleRef = *index;
529
530			//find W/Z mother
531			while ( mcParticleRef->mother()!=0) {
532
533			if (parentItr>=2) break;
534			parentItr++;
535			mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
536			if (mcParticleRef->mother()==0) break;
537
538			//if tau
539			if (abs((mcParticleRef)->pdgId())==15 ) {
540			Tauchild=true;
541			break;
542			}
543
544			// cout << "mother: " << (mcParticleRef)->pdg_id()
545			// <<" pt="<< c_p4.pt() << " eta="<<c_p4.eta()
546			// <<" status=" <<mcParticleRef->status() <<"\n";
547
548			if (abs((mcParticleRef)->pdgId())==23 ) {
549			Zchild=true;
550			break;
551			}
552			if (abs((mcParticleRef)->pdgId())==24 ) {
553			Wchild=true;
554			break;
555			}
556			}
557
558			if (maxPt >0) {
559			if (Tauchild) {
560			parentItr = 0;
561			while (mcParticleRef->mother()!=0) {
562			if (parentItr>=2) break;
563			parentItr++;
564			mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
565			if (mcParticleRef->mother()==0) break;
566
567			if (abs((mcParticleRef)->pdgId())==23 ) {
568			if (MCleptZ1_pt == float(mcParticleRef->pt())
569			&& MCleptZ1_eta == float(mcParticleRef->eta())
570			&& MCleptZ1_phi == float(mcParticleRef->phi()) )
571			MC_tauDecayTypeZ1 = i;
572			else MC_tauDecayTypeZ2 = i;
573			break;
574			}
575			if (abs((mcParticleRef)->pdgId())==24 ) {
576			MC_tauDecayTypeW = i;
577			break;
578			}
579			}
580			}
581
582			if (Wchild) {
583			MCleptW_pdgid=(*index)->pdgId();
584			MCleptW_pt=(float)(*index)->pt();
585			MCleptW_eta=(float)(*index)->eta();
586			MCleptW_phi=(float)(*index)->phi();
587			if (abs(MCleptW_pdgid)==15) MC_tauDecayTypeW =3;//default to hadronic decay
588			}
589			if (Zchild) {
590			if (firstZlept) {
591			firstZlept=false;
592			MCleptZ1_pdgid=(*index)->pdgId();
593			MCleptZ1_pt=(float)(*index)->pt();
594			MCleptZ1_eta=(float)(*index)->eta();
595			MCleptZ1_phi=(float)(*index)->phi();
596			if (abs(MCleptZ1_pdgid)==15) MC_tauDecayTypeZ1 =3;
597			}
598			else {
599			MCleptZ2_pdgid=(*index)->pdgId();
600			MCleptZ2_pt=(float)(*index)->pt();
601			MCleptZ2_eta=(float)(*index)->eta();
602			MCleptZ2_phi=(float)(*index)->phi();
603			if (abs(MCleptZ2_pdgid)==15) MC_tauDecayTypeZ2 =3;
604			}
605			}
606			StableLeptons[i]->erase(index);
607			}
608			}
609			}
610	vuko	1.1
611	smorovic	1.8	}
612	vuko	1.1	//define this as a plug-in
613	vuko	1.2	// DEFINE_FWK_MODULE(WZAnalyzer);