CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/ElectroWeakAnalysis/ZMuMu/plugins/zPdfUnc.cc

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 zPdfUnc : public edm::EDAnalyzer {
00006 public:
00007   zPdfUnc(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_, pdfweights_;
00012   unsigned int pdfmember_ , nbinsMass_, nbinsPt_, nbinsAng_ ;
00013   double massMax_, ptMax_, angMax_;
00014   double  accPtMin_,accMassMin_,accMassMax_, accEtaMin_, accEtaMax_,  accMassMinDenominator_ ;
00015   TH1F *h_nZ_;
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   double nAcc_, nAccReW_,  nBothMuHasZHasGrandMa_;
00021   std::string filename_;
00022 
00023 };
00024 
00025 #include "DataFormats/Candidate/interface/Candidate.h"
00026 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00027 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00028 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00029 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00030 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00031 
00032 #include "HepMC/WeightContainer.h"
00033 #include "HepMC/GenEvent.h"
00034 #include "FWCore/Framework/interface/Event.h"
00035 #include "DataFormats/Common/interface/Handle.h"
00036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00037 #include "FWCore/ServiceRegistry/interface/Service.h"
00038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00039 #include <cmath>
00040 #include <iostream>
00041 #include <fstream>
00042 
00043 using namespace std;
00044 using namespace reco;
00045 using namespace edm;
00046 
00047 zPdfUnc::zPdfUnc(const ParameterSet& pset) :
00048   gen_(pset.getParameter<InputTag>("genParticles")),
00049   pdfweights_(pset.getParameter<InputTag>("pdfweights")),
00050   pdfmember_(pset.getUntrackedParameter<unsigned int>("pdfmember")),
00051   nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
00052   nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
00053   nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
00054   massMax_(pset.getUntrackedParameter<double>("massMax")),
00055   ptMax_(pset.getUntrackedParameter<double>("ptMax")),
00056   angMax_(pset.getUntrackedParameter<double>("angMax")), 
00057   accPtMin_(pset.getUntrackedParameter<double>("accPtMin")), 
00058   accMassMin_(pset.getUntrackedParameter<double>("accMassMin")),
00059   accMassMax_(pset.getUntrackedParameter<double>("accMassMax")),
00060   accEtaMin_(pset.getUntrackedParameter<double>("accEtaMin")),
00061   accEtaMax_(pset.getUntrackedParameter<double>("accEtaMax")),
00062   accMassMinDenominator_(pset.getUntrackedParameter<double>("accMassMinDenominator")),
00063   isMCatNLO_(pset.getUntrackedParameter<bool>("isMCatNLO")),
00064   filename_(pset.getUntrackedParameter<std::string>("outfilename")) { 
00065   cout << ">>> Z Histogrammer constructor" << endl;
00066   Service<TFileService> fs;
00067   
00068   TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
00069   h_nZ_ = ZMCHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
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   nAccReW_ =0;
00087   nBothMuHasZHasGrandMa_=0;
00088   
00089 }
00090 
00091 void zPdfUnc::analyze(const edm::Event& event, const edm::EventSetup& setup) { 
00092   cout << ">>> Z Histogrammer analyze" << endl;
00093 
00094   Handle<GenParticleCollection> gen;
00095   Handle< vector<double> > pdfweights;
00096 
00097   event.getByLabel(gen_, gen);
00098  event.getByLabel(pdfweights_, pdfweights );
00099 
00100  
00101 
00102    // get weight and fill it to histogram
00103  std::vector<double> weights = (*pdfweights);
00104 
00105   unsigned int nmembers = weights.size();
00106 
00107 
00108   double weight =  weights[pdfmember_];
00109   
00110   h_weight_histo->Fill(weight); 
00111   
00112   cout << "found nmember: " << nmembers << endl;
00113   cout << "taken the member n " << pdfmember_ << " weight= " << weight <<   endl;
00114 
00115   
00116 
00117   std::vector<GenParticle> muons;
00118 
00119     
00120     for(unsigned int i = 0; i < gen->size(); ++i){ 
00121       const GenParticle & muMC  = (*gen)[i];
00122       // filling only muons coming form Z
00123       if (abs(muMC.pdgId())==13 &&  muMC.status()==1  && muMC.numberOfMothers()>0) {   
00124         if (muMC.mother()->numberOfMothers()> 0 ){
00125           cout << "I'm getting a muon \n" 
00126                << "with " << "muMC.numberOfMothers()  " <<  muMC.numberOfMothers() << "\n the first mother has pdgId " << muMC.mother()->pdgId() 
00127                << "with " << "muMC.mother()->numberOfMothers()  " <<  muMC.mother()->numberOfMothers()<< "\n the first grandma has pdgId " << muMC.mother()->mother()->pdgId()<<endl;
00128           if (muMC.mother()->mother()->pdgId() ==23 ) muons.push_back(muMC);
00129         }
00130       }
00131     }
00132 
00133   cout << "finally I selected " << muons.size() << " muons" << endl;
00134  
00135 // if there are at least two muons, 
00136    // calculate invarant mass of first two and fill it into histogram
00137 
00138   double inv_mass = 0.0;
00139    double Zpt_ = 0.0;
00140    double Zeta_ = 0.0;
00141    double Ztheta_ = 0.0;
00142    double Zphi_ = 0.0;
00143    double Zrapidity_ = 0.0;
00144   
00145    if(muons.size()>1) {
00146      if (muons[0].mother()->mother()->pdgId()==23 && muons[1].mother()->mother()->pdgId()==23) nBothMuHasZHasGrandMa_ ++;
00147      math::XYZTLorentzVector tot_momentum(muons[0].p4());
00148      math::XYZTLorentzVector mom2(muons[1].p4());
00149      tot_momentum += mom2;
00150      inv_mass = sqrt(tot_momentum.mass2());
00151 
00152    Zpt_=tot_momentum.pt();
00153      Zeta_ = tot_momentum.eta();
00154      Ztheta_ = tot_momentum.theta();
00155      Zphi_ = tot_momentum.phi();
00156     Zrapidity_ = tot_momentum.Rapidity();
00157 
00158 
00159      
00160      // IMPORTANT: use the weight of the event ...
00161      
00162     double weight_sign = weight;
00163       //double weight_sign = 1.    ;
00164     // fill the denominator numbers only if mass>accMassMinDenominator and go on in that case
00165     if (inv_mass > accMassMinDenominator_ ) {  
00166       h_mZMC_->Fill(inv_mass,weight_sign);
00167       h_ptZMC_->Fill(Zpt_,weight_sign);
00168       h_etaZMC_->Fill(Zeta_,weight_sign);
00169       h_thetaZMC_->Fill(Ztheta_,weight_sign);
00170       h_phiZMC_->Fill(Zphi_,weight_sign);
00171       h_rapidityZMC_-> Fill (Zrapidity_,weight_sign );    
00172     
00173       double pt1 = muons[0].pt();
00174       double pt2 = muons[1].pt();
00175       double eta1 = muons[0].eta();
00176       double eta2 = muons[1].eta();
00177       if(pt1>pt2) {
00178         hardpt->Fill(pt1,weight_sign);
00179         softpt->Fill(pt2,weight_sign);
00180         hardeta->Fill(eta1,weight_sign);
00181         softeta->Fill(eta2,weight_sign);
00182       } else {
00183         hardpt->Fill(pt2,weight_sign);
00184         softpt->Fill(pt1,weight_sign);       
00185         hardeta->Fill(eta2,weight_sign);
00186         softeta->Fill(eta1,weight_sign);
00187       }
00188       
00189       //evaluating the geometric acceptance  
00190       if ( pt1 >= accPtMin_  && pt2 >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_)  {
00191         nAcc_ ++;
00192         nAccReW_ +=weight; 
00193       }  
00194     }
00195    }
00196 } 
00197   
00198 
00199 void zPdfUnc::endJob() {
00200   cout << " number of events accepted :" << nAcc_ << endl;
00201   cout << " number of events accepted reweigthed :" << nAccReW_ << endl;
00202   double nev = h_mZMC_->GetEntries();
00203   double nev_weigthed = h_mZMC_->Integral( ) ;
00204   cout << " number of total events :" << nev   << endl;
00205     cout << " number of total weighted events :" << nev_weigthed << endl;
00206  cout << " number of cases in which BothMuHasZHasGrandMa :" << nBothMuHasZHasGrandMa_  << endl;
00207   double eff = (double)nAcc_ / (double) h_mZMC_->GetEntries();
00208  double eff_rew = (double)nAccReW_ / (double) h_mZMC_->Integral(0,nbinsMass_ +1 );
00209   double err = sqrt( eff * (1. - eff) / (double) h_mZMC_->GetEntries() );
00210   double err_rew = sqrt( eff_rew * (1. - eff_rew) / (double) h_mZMC_->Integral(0,nbinsMass_ +1 ) );
00211   cout << " geometric acceptance: " << eff << "+/-" << err << endl; 
00212   cout << " geometric acceptance reweighted: " << eff_rew << "+/-" << err_rew << endl; 
00213 
00214   ofstream myfile;
00215   myfile.open(filename_.c_str(), std::ios::app);  
00216   myfile<< eff << " "<< eff_rew << " " << nev << " " << nev_weigthed << endl ;
00217   myfile.close(); 
00218 }
00219 #include "FWCore/Framework/interface/MakerMacros.h"
00220 
00221 DEFINE_FWK_MODULE(zPdfUnc);
00222