Reference Manual

Go to the documentation of this file.
00001 #include "FWCore/Framework/interface/EDAnalyzer.h"
00002 #include "FWCore/Utilities/interface/InputTag.h"
00003 #include "TH1.h"
00004 
00005 class ZLONLOHistogrammer : public edm::EDAnalyzer {
00006 public:
00007   ZLONLOHistogrammer(const edm::ParameterSet& pset); 
00008 private:
00009   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00010   virtual void endJob();
00011   edm::InputTag   gen_, weights_;
00012   unsigned int nbinsMass_, nbinsPt_, nbinsAng_;
00013   double massMax_, ptMax_, angMax_;
00014   double  accPtMin_,accMassMin_,accMassMax_, accEtaMin_, accEtaMax_;
00015   TH1F *h_nZ_, *h_mZ_, *h_ptZ_, *h_phiZ_, *h_thetaZ_, *h_etaZ_, *h_rapidityZ_;
00016   TH1F *h_mZMC_, *h_ptZMC_, *h_phiZMC_, *h_thetaZMC_, *h_etaZMC_, *h_rapidityZMC_;
00017   TH1F *hardpt, *softpt, * hardeta, *softeta;
00018   TH1F * h_weight_histo; 
00019   bool isMCatNLO_; 
00020   unsigned int nAcc_, nBothMuHasZHasGrandMa_;
00021 };
00022 
00023 #include "DataFormats/Candidate/interface/Candidate.h"
00024 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00025 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00026 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00027 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00028 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00029 
00030 #include "HepMC/WeightContainer.h"
00031 #include "HepMC/GenEvent.h"
00032 #include "FWCore/Framework/interface/Event.h"
00033 #include "DataFormats/Common/interface/Handle.h"
00034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00035 #include "FWCore/ServiceRegistry/interface/Service.h"
00036 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00037 #include <cmath>
00038 #include <iostream>
00039 
00040 using namespace std;
00041 using namespace reco;
00042 using namespace edm;
00043 
00044 ZLONLOHistogrammer::ZLONLOHistogrammer(const ParameterSet& pset) :
00045   gen_(pset.getParameter<InputTag>("genParticles")),
00046   weights_(pset.getParameter<InputTag>("weights")),
00047   nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
00048   nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
00049   nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
00050   massMax_(pset.getUntrackedParameter<double>("massMax")),
00051   ptMax_(pset.getUntrackedParameter<double>("ptMax")),
00052   angMax_(pset.getUntrackedParameter<double>("angMax")), 
00053   accPtMin_(pset.getUntrackedParameter<double>("accPtMin")), 
00054   accMassMin_(pset.getUntrackedParameter<double>("accMassMin")),
00055   accMassMax_(pset.getUntrackedParameter<double>("accMassMax")),
00056   accEtaMin_(pset.getUntrackedParameter<double>("accEtaMin")),
00057   accEtaMax_(pset.getUntrackedParameter<double>("accEtaMax")),
00058   isMCatNLO_(pset.getUntrackedParameter<bool>("isMCatNLO"))  { 
00059   cout << ">>> Z Histogrammer constructor" << endl;
00060   Service<TFileService> fs;
00061   TFileDirectory ZHisto = fs->mkdir( "ZRecHisto" );
00062   TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
00063   h_nZ_ = ZHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
00064   h_mZ_ = ZHisto.make<TH1F>("ZMass", "Z mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
00065   h_ptZ_ = ZHisto.make<TH1F>("ZPt", "Z p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00066   h_phiZ_ = ZHisto.make<TH1F>("ZPhi", "Z #phi", nbinsAng_,  -angMax_, angMax_);
00067   h_thetaZ_ = ZHisto.make<TH1F>("Ztheta", "Z #theta", nbinsAng_,  0, angMax_);
00068   h_etaZ_ = ZHisto.make<TH1F>("ZEta", "Z #eta", nbinsAng_,  -angMax_, angMax_);
00069   h_rapidityZ_ = ZHisto.make<TH1F>("ZRapidity", "Z rapidity", nbinsAng_,  -angMax_, angMax_);
00070   h_weight_histo  = ZMCHisto.make<TH1F>("weight_histo","weight_histo",20,-10,10);       
00071 
00072   h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
00073   h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00074  hardpt = ZMCHisto.make<TH1F>("hardpt", "hard muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00075  softpt = ZMCHisto.make<TH1F>("softpt", "soft muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00076 
00077  
00078   h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_,  -angMax_, angMax_);
00079   h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_,  0, angMax_);
00080   h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_,  -angMax_, angMax_);
00081   h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC y", nbinsAng_,  -angMax_, angMax_);
00082 
00083   hardeta = ZMCHisto.make<TH1F>("hard muon eta", "hard muon #eta", nbinsAng_,  -angMax_, angMax_);
00084   softeta = ZMCHisto.make<TH1F>("soft muon eta", "soft muon #eta", nbinsAng_,  -angMax_, angMax_);
00085   nAcc_=0;
00086   nBothMuHasZHasGrandMa_=0;
00087   
00088 }
00089 
00090 void ZLONLOHistogrammer::analyze(const edm::Event& event, const edm::EventSetup& setup) { 
00091   cout << ">>> Z Histogrammer analyze" << endl;
00092 
00093   Handle<GenParticleCollection> gen;
00094   Handle<double> weights;
00095 
00096   event.getByLabel(gen_, gen);
00097  event.getByLabel(weights_, weights );
00098 
00099  
00100 
00101    // get weight and fill it to histogram
00102   
00103   double weight = * weights;
00104   if(!weight) weight=1.;
00105   h_weight_histo->Fill(weight); 
00106   
00107    if(isMCatNLO_) {
00108     weight > 0 ?  weight=1. : weight=-1.;
00109   }
00110 
00111   
00112 
00113   std::vector<GenParticle> muons;
00114   if (!isMCatNLO_){
00115     // LO....
00116     for(unsigned int i = 0; i < gen->size(); ++i){ 
00117       const GenParticle & muMC  = (*gen)[i];
00118       // filling only muons coming form Z
00119       if (abs(muMC.pdgId())==13 &&  muMC.status()==1  && muMC.numberOfMothers()>0) {   
00120         if (muMC.mother()->numberOfMothers()> 0 ){
00121           cout << "I'm getting a muon \n" 
00122                << "with " << "muMC.numberOfMothers()  " <<  muMC.numberOfMothers() << "\n the first mother has pdgId " << muMC.mother()->pdgId() 
00123                << "with " << "muMC.mother()->numberOfMothers()  " <<  muMC.mother()->numberOfMothers()<< "\n the first grandma has pdgId " << muMC.mother()->mother()->pdgId()<<endl;
00124           if (muMC.mother()->mother()->pdgId() ==23 ) muons.push_back(muMC);
00125         }
00126       }
00127     }
00128   } else {
00129     // NLO
00130     for(unsigned int i = 0; i < gen->size(); ++i){ 
00131       const GenParticle & muMC  = (*gen)[i];
00132       if (abs(muMC.pdgId())==13 &&  muMC.status()==1  && muMC.numberOfMothers()>0) {                                                         if (muMC.mother()->numberOfMothers()> 0 ){
00133         cout << "I'm getting a muon \n" 
00134              << "with " << "muMC.numberOfMothers()  " <<  muMC.numberOfMothers() << "\n the first mother has pdgId " << muMC.mother()->pdgId() 
00135              << "with " << "muMC.mother()->numberOfMothers()  " <<  muMC.mother()->numberOfMothers()<< "\n the first grandma has pdgId " << muMC.mother()->mother()->pdgId()<<endl;
00136         // filling withoput requiring that the grandma is a Z...... sometimes the grandma are still muons, otherwise those are fake muons, but the first two are always the desired muons.... 
00137         muons.push_back(muMC);
00138       }   
00139       }
00140     }
00141   }
00142 
00143   cout << "finally I selected " << muons.size() << " muons" << endl;
00144 
00145  
00146 
00147 
00148 
00149 // if there are at least two muons, 
00150    // calculate invarant mass of first two and fill it into histogram
00151   
00152   //if  isMCatNLO_......
00153 
00154 
00155 
00156 
00157  double inv_mass = 0.0;
00158    double Zpt_ = 0.0;
00159    double Zeta_ = 0.0;
00160    double Ztheta_ = 0.0;
00161    double Zphi_ = 0.0;
00162    double Zrapidity_ = 0.0;
00163   
00164    if(muons.size()>1) {
00165      if (muons[0].mother()->mother()->pdgId()==23 && muons[1].mother()->mother()->pdgId()==23) nBothMuHasZHasGrandMa_ ++;
00166      math::XYZTLorentzVector tot_momentum(muons[0].p4());
00167      math::XYZTLorentzVector mom2(muons[1].p4());
00168      tot_momentum += mom2;
00169      inv_mass = sqrt(tot_momentum.mass2());
00170      Zpt_=tot_momentum.pt();
00171      Zeta_ = tot_momentum.eta();
00172      Ztheta_ = tot_momentum.theta();
00173      Zphi_ = tot_momentum.phi();
00174     Zrapidity_ = tot_momentum.Rapidity();
00175 
00176 
00177      
00178      // IMPORTANT: use the weight of the event ...
00179      
00180       double weight_sign = (weight > 0) ? 1. : -1.;
00181       //double weight_sign = 1.    ;
00182       h_mZMC_->Fill(inv_mass,weight_sign);
00183       h_ptZMC_->Fill(Zpt_,weight_sign);
00184       h_etaZMC_->Fill(Zeta_,weight_sign);
00185       h_thetaZMC_->Fill(Ztheta_,weight_sign);
00186       h_phiZMC_->Fill(Zphi_,weight_sign);
00187       h_rapidityZMC_-> Fill (Zrapidity_,weight_sign );    
00188 
00189       double pt1 = muons[0].pt();
00190       double pt2 = muons[1].pt();
00191       double eta1 = muons[0].eta();
00192       double eta2 = muons[1].eta();
00193       
00194 
00195 
00196      if(pt1>pt2) {
00197        hardpt->Fill(pt1,weight_sign);
00198        softpt->Fill(pt2,weight_sign);
00199        hardeta->Fill(eta1,weight_sign);
00200        softeta->Fill(eta2,weight_sign);
00201      } else {
00202        hardpt->Fill(pt2,weight_sign);
00203        softpt->Fill(pt1,weight_sign);       
00204        hardeta->Fill(eta2,weight_sign);
00205        softeta->Fill(eta1,weight_sign);
00206      }
00207    
00208 
00209    //evaluating the geometric acceptance  
00210    if ( pt1 >= accPtMin_  && pt2 >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) nAcc_++;
00211    
00212           
00213    }
00214 
00215 } 
00216   
00217 
00218 void ZLONLOHistogrammer::endJob() {
00219   cout << " number of events accepted :" << nAcc_ << endl;
00220   cout << " number of total events :" << h_mZMC_->GetEntries()  << endl;
00221    cout << " number of cases in which BothMuHasZHasGrandMa :" << nBothMuHasZHasGrandMa_  << endl;
00222   double eff = (double)nAcc_ / (double) h_mZMC_->GetEntries();
00223   double err = sqrt( eff * (1. - eff) / (double) h_mZMC_->GetEntries() );
00224   cout << " geometric acceptance: " << eff << "+/-" << err << endl; 
00225 }
00226 #include "FWCore/Framework/interface/MakerMacros.h"
00227 
00228 DEFINE_FWK_MODULE(ZLONLOHistogrammer);
00229
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/ElectroWeakAnalysis/ZMuMu/plugins/ZLONLOHistogrammer.cc
