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.18 2007/11/30 17:24:50 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 "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");
  
  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

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

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

    }   

    // MET
    reco::Particle::LorentzVector Muon_p(0,0,0,0);
    for (int i=0; i<looseMuonCands->size();i++){
      Muon_p=Muon_p+(*looseMuonCands)[i].p4();
    }
    reco::Particle::LorentzVector MET=computeMET(iEvent, iSetup)-Muon_p;  // substract p from muons from MET 
    MET_energy=MET.energy();
    MET_pt=MET.pt();
    MET_eta=MET.eta();
    

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

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

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

      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(W,Z)
      cout<<"theWlepton->phi()="<<theWlepton->phi() <<endl;
      cout<<"theZCandidate->phi()="<<theZCandidate->phi() <<endl;
      cout<<"dphi()="<<theWlepton->phi()-theZCandidate->phi() <<endl;
      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)  { /* it's a valid electron */ 
        cout << "lepton is an electron \n" << endl;
        elWTree->Fill();
      }
    }

  }
  
  //  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;
  //  for(CandidateCollection::const_iterator jet = jetCands->begin();
  for(CaloJetCollection::const_iterator jet = jetCands->begin();
      jet != jetCands->end(); jet++) {
    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();

}


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

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

  leptonProperties["Wmu"] ->CreateBranches(wzTree, "Wmu_");
  leptonProperties["Zmu1"] ->CreateBranches(wzTree, "Zmul_");
  leptonProperties["Zmu2"] ->CreateBranches(wzTree, "Zmu2_");

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

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


}


reco::Particle::LorentzVector WZAnalyzer::computeMET(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  
  using namespace edm;
  using namespace reco;
  using namespace std;
  
  Handle<CaloMETCollection> mets;
  iEvent.getByLabel("met", mets);
  cout << "# METs=" << mets->size() << endl;
  for (CaloMETCollection::const_iterator met = mets->begin();
       met!= mets->end(); met++) {
    cout << "MET:  E = " << met->energy()
         << "\t Pt = " << met->pt()
         <<endl;
  }
  return mets->begin()->p4();
  
}

////////////////////////
// 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);
          cout << "dr by hand = " << dR_byhand <<" i="<< i<<" j="<<j <<endl;
          double dR=dR_byhand;
          
          if (dR<0.15){
            n1++;
            matching_Infos.push_back(MatchingInfo(dR,i,j));
          }
          
        }
    } 

  cout<<"table for matching reco&gen;"<< " pdgid="<<pdgid <<endl;
   for (int i=0; i<matching_Infos.size();i++)
     {
       cout <<"dr="<<matching_Infos[i].deltaR << " i="<<matching_Infos[i].genMuon <<" j="<<matching_Infos[i].recoMuon <<endl;
     }

  sort(matching_Infos.begin(),matching_Infos.end(),betterMatch);  // sorting (dR,i,j)   

  cout<<"sorted table;" <<endl;

   for (int i=0; i<matching_Infos.size();i++)
     {
       cout <<"dr="<<matching_Infos[i].deltaR << " i="<<matching_Infos[i].genMuon <<" j="<<matching_Infos[i].recoMuon <<endl;
     }

  cout<<"MC.size()="<<genparticles.size() <<"  reco->size()="<<cands->size() <<endl;
   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;
         }
       }
     }
   }
   cout << "final table:"<<endl;
   for (int i=0; i<matching_Infos.size();i++)
     {
       cout <<"dr="<<matching_Infos[i].deltaR << " i="<<matching_Infos[i].genMuon <<" j="<<matching_Infos[i].recoMuon <<endl;
     }

   // Now put result in vector of pairs

  
   //   vector<  pair<const GenParticleCandidate*, 
   //     const reco::Candidate *> >  * leptonMatching = 0;

//    if (pdgid == 11) {
//      leptonMatching = &electronMatching;
//    } else if (pdgid == 13) {
//      leptonMatching = &muonMatching;
//    }

//   if (leptonMatching == 0) return;
//   leptonMatching->clear();

   leptonMatching[pdgid].clear();

//    class matchedParticles{
//    private:
//      pair<const GenParticleCandidate*, const reco::Candidate *> pair;
//      float delta_R;
//    public:
//      matchedParticles();
//      //     matchedParticles( pair<const GenParticleCandidate*, const reco::Candidate *> pair2, float delta_R2){ pair=pair2;delta_R2=delta_R;};
//      matchedParticles( const GenParticleCandidate* genPar, const reco::Candidate * recoPar, float delta_R2) {pair.first= genPar;pair.second= recoPar ;delta_R=delta_R2;}
//    };

   //   matchedParticles::matchedParticles();
   
   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);
     
     //     leptonMatching[pdgid].push_back(pair);
   }

   
}


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

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

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 (int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
        const reco::Candidate * nesto=iter->second[j].pair.second;
        if (overlap(*nesto,*daughter)){
          cout <<" found matched genparticle from vector<class<..>>" <<endl;
          //      return matched_genParticle=iter->second[j].pair.first;
          matched_genParticle=iter->second[j].pair.first;
          //      delta_r=iter->second[j].delta_R;
        }
      }
    }
  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 (int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
        const reco::Candidate * nesto=iter->second[j].pair.second;
        if (overlap(*nesto,*daughter)){
          //      cout <<" found matched genparticle from vector<class<..>>" <<endl;
          //      return matched_genParticle=iter->second[j].pair.first;
          delta_r=iter->second[j].delta_R;
        }
      }
    }
  return delta_r;
}


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

  //collections
  cout<<" pocetak WZAnalyzer::collectMCsummary" <<endl;
  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 ) {
  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;
                        }

                //      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);
                }
        }
  }
  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.19
Committed:	Tue Dec 4 12:38:22 2007 UTC (17 years, 5 months ago) by senka
Content type:	text/plain
Branch:	MAIN
Changes since 1.18:	+16 -22 lines
Log Message:	segmentation fault solved