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 #include "TRandom3.h"
00005 
00006 
00007 class ZMuPtScaleAnalyzer : public edm::EDAnalyzer {
00008 public:
00009   ZMuPtScaleAnalyzer(const edm::ParameterSet& pset); 
00010 private:
00011   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00012   virtual void endJob();
00013   edm::InputTag   gen_;
00014   unsigned int nbinsMass_, nbinsPt_, nbinsAng_;
00015   double massMax_, ptMax_, angMax_;
00016   double  accPtMin_,accMassMin_,accMassMax_, accMassMinDen_,accMassMaxDen_, accEtaMin_, accEtaMax_, ptScale_;
00017   TH1F *h_nZ_, *h_mZ_, *h_ptZ_, *h_phiZ_, *h_thetaZ_, *h_etaZ_, *h_rapidityZ_;
00018   TH1F *h_mZMC_, *h_ptZMC_, *h_phiZMC_, *h_thetaZMC_, *h_etaZMC_, *h_rapidityZMC_;
00019   TH1F *hardpt, *softpt, * hardeta, *softeta;
00020   unsigned int nAcc_,nAccPtScaleP_,nAccPtScaleN_, nAccPtScaleSmearedFlat_ , nAccPtScaleSmearedGaus_,  nBothMuHasZHasGrandMa_;
00021  int  muPdgStatus_;
00022 };
00023 
00024 #include "DataFormats/Candidate/interface/Candidate.h"
00025 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00026 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00027 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00028 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00029 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00030 
00031 #include "HepMC/WeightContainer.h"
00032 #include "HepMC/GenEvent.h"
00033 #include "FWCore/Framework/interface/Event.h"
00034 #include "DataFormats/Common/interface/Handle.h"
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 #include "FWCore/ServiceRegistry/interface/Service.h"
00037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00038 #include <cmath>
00039 #include <iostream>
00040 
00041 using namespace std;
00042 using namespace reco;
00043 using namespace edm;
00044 
00045 ZMuPtScaleAnalyzer::ZMuPtScaleAnalyzer(const ParameterSet& pset) :
00046   gen_(pset.getParameter<InputTag>("genParticles")),
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   accMassMinDen_(pset.getUntrackedParameter<double>("accMassMinDen")),
00057   accMassMaxDen_(pset.getUntrackedParameter<double>("accMassMaxDen")),
00058   accEtaMin_(pset.getUntrackedParameter<double>("accEtaMin")),
00059   accEtaMax_(pset.getUntrackedParameter<double>("accEtaMax")),
00060   ptScale_(pset.getUntrackedParameter<double>("ptScale")),
00061   muPdgStatus_(pset.getUntrackedParameter<int>("muPdgStatus")){ 
00062 
00063   cout << ">>> Z Histogrammer constructor" << endl;
00064   Service<TFileService> fs;
00065   TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
00066  
00067   h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_,  0, massMax_);
00068   h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00069  hardpt = ZMCHisto.make<TH1F>("hardpt", "hard muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00070  softpt = ZMCHisto.make<TH1F>("softpt", "soft muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
00071 
00072  
00073   h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_,  -angMax_, angMax_);
00074   h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_,  0, angMax_);
00075   h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_,  -angMax_, angMax_);
00076   h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC y", nbinsAng_,  -angMax_, angMax_);
00077 
00078   hardeta = ZMCHisto.make<TH1F>("hard muon eta", "hard muon #eta", nbinsAng_,  -angMax_, angMax_);
00079   softeta = ZMCHisto.make<TH1F>("soft muon eta", "soft muon #eta", nbinsAng_,  -angMax_, angMax_);
00080   nAcc_=0;
00081   nAccPtScaleP_=0;
00082   nAccPtScaleN_=0;
00083   nAccPtScaleSmearedFlat_=0;
00084   nAccPtScaleSmearedGaus_=0;
00085   nBothMuHasZHasGrandMa_=0;}
00086 
00087 
00088 
00089 void ZMuPtScaleAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) { 
00090   cout << ">>> Z HistogrammerZLONLOHistogrammer.cc analyze" << endl;
00091 
00092   Handle<GenParticleCollection> gen;
00093   Handle<double> weights;
00094 
00095   event.getByLabel(gen_, gen);
00096 
00097 
00098  
00099 
00100    // get weight and fill it to histogram
00101   
00102 
00103 
00104 
00105   std::vector<GenParticle> muons;
00106 
00107   double mZGen = -100;
00108 
00109     for(unsigned int i = 0; i < gen->size(); ++i){ 
00110       const GenParticle & muMC  = (*gen)[i];
00111       // filliScaledPng only muons coming form Z
00112       if (abs(muMC.pdgId())==13 &&  muMC.status()==muPdgStatus_  && muMC.numberOfMothers()>0) {   
00113         
00114           cout << "I'm getting a muon \n" 
00115                << "with " << "muMC.numberOfMothers()  " <<  muMC.numberOfMothers() << "\n the first mother has pdgId " << muMC.mother()->pdgId() 
00116                << "with " << "muMC.mother()->numberOfMothers()  " <<  muMC.mother()->numberOfMothers()<< "\n the first grandma has pdgId " << muMC.mother()->mother()->pdgId()<<endl;
00117                cout << "with  muMC.eta() " <<  muMC.eta()<<endl;
00118                muons.push_back(muMC);
00119       }
00120           // introducing here the gen mass cut......................
00121                /*
00122    if (muPdgStatus_ ==1) {
00123             mZGen = muMC.mother()->mother()->mass(); 
00124             if (muMC.mother()->mother()->pdgId() ==23  && mZGen>accMassMinDen_ && mZGen<accMassMaxDen_ ) muons.push_back(muMC);}
00125         }
00126           if (muPdgStatus_ ==3) {
00127              mZGen = muMC.mother()->mass(); 
00128             if (muMC.mother()->pdgId() ==23  && mZGen>accMassMinDen_ && mZGen<accMassMaxDen_ ) muons.push_back(muMC);}
00129         }
00130 */
00131       
00132       
00133       const GenParticle & zMC  = (*gen)[i];
00134             if (zMC.pdgId()==23 &&  zMC.status()==3 &&zMC.numberOfDaughters()>1  ) {
00135         mZGen = zMC.mass();
00136         cout << "I'm selecting a Z MC with mass " << mZGen << endl; 
00137           if(mZGen>accMassMinDen_ && mZGen<accMassMaxDen_ ) h_mZMC_->Fill(mZGen); 
00138       
00139       
00140       
00141       }
00142       }
00143     
00144    
00145   cout << "finally I selected " << muons.size() << " muons" << endl;
00146 
00147  
00148 
00149 
00150 
00151 // if there are at least two muons, 
00152    // calculate invarant mass of first two and fill it into histogram
00153   
00154 
00155 
00156  double inv_mass = 0.0;
00157    double Zpt_ = 0.0;
00158    double Zeta_ = 0.0;
00159    double Ztheta_ = 0.0;
00160    double Zphi_ = 0.0;
00161    double Zrapidity_ = 0.0;
00162   
00163    if(muons.size()>1) {
00164 
00165 
00166     if (muons[0].mother()->mother()->pdgId()==23 && muons[1].mother()->mother()->pdgId()==23) nBothMuHasZHasGrandMa_ ++;
00167      math::XYZTLorentzVector tot_momentum(muons[0].p4());
00168      math::XYZTLorentzVector mom2(muons[1].p4());
00169      tot_momentum += mom2;
00170      inv_mass = sqrt(tot_momentum.mass2());
00171      Zpt_=tot_momentum.pt();
00172      Zeta_ = tot_momentum.eta();
00173      Ztheta_ = tot_momentum.theta();
00174      Zphi_ = tot_momentum.phi();
00175     Zrapidity_ = tot_momentum.Rapidity();
00176 
00177 
00178      
00179 
00180       double weight_sign = 1.    ;
00181      
00182       //h_mZMC_->Fill(inv_mass); 
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 
00197 
00198      if(pt1>pt2) {
00199        hardpt->Fill(pt1,weight_sign);
00200        softpt->Fill(pt2,weight_sign);
00201        hardeta->Fill(eta1,weight_sign);
00202        softeta->Fill(eta2,weight_sign);
00203      } else {
00204        hardpt->Fill(pt2,weight_sign);
00205        softpt->Fill(pt1,weight_sign);       
00206        hardeta->Fill(eta2,weight_sign);
00207        softeta->Fill(eta1,weight_sign);
00208      }
00209    
00210  //evaluating the geometric acceptance  
00211    if ( pt1 >= accPtMin_  && pt2 >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) nAcc_++;
00212    
00213            
00214    cout << "pt1" <<  pt1 << endl;
00215    
00216      // scaling the muon pt
00217    double  pt1ScaledP = pt1* ( 1. + ptScale_);
00218   cout << "pt1 ScaledP of " <<  ( 1. + ptScale_) << endl;
00219    cout << "pt1ScaledP" <<  pt1ScaledP << endl;
00220 
00221    double  pt2ScaledP = pt2 * ( 1. + ptScale_);
00222   
00223 
00224    //evaluating the geometric acceptance  
00225    if ( pt1ScaledP >= accPtMin_  && pt2ScaledP >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) 
00226 nAccPtScaleP_++;
00227    
00228 
00229 
00230   
00231      // scaling the muon pt
00232    double  pt1ScaledN = pt1* ( 1. - ptScale_);
00233    double  pt2ScaledN = pt2 * ( 1. - ptScale_);
00234   
00235 
00236    //evaluating the geometric acceptance  
00237    if ( pt1ScaledN >= accPtMin_  && pt2ScaledN >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) 
00238 nAccPtScaleN_++;
00239         
00240 
00241   // scaling the muon pt
00242    TRandom3 f;
00243    f.SetSeed(123456789);
00244    double  pt1SmearedFlat = pt1* ( 1. + ptScale_ * f.Uniform() );
00245    double  pt2SmearedFlat = pt2 * ( 1. + ptScale_ * f.Uniform() ) ;
00246   
00247 
00248    //evaluating the geometric acceptance  
00249    if ( pt1SmearedFlat >= accPtMin_  && pt2SmearedFlat >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) 
00250  nAccPtScaleSmearedFlat_++;
00251           
00252 
00253 // scaling the muon pt
00254    TRandom3 ff;
00255    ff.SetSeed(123456789);
00256    double  pt1SmearedGaus = pt1* ( 1. + ptScale_ * f.Gaus() );
00257    double  pt2SmearedGaus = pt2 * ( 1. + ptScale_ * f.Gaus() ) ;
00258   
00259 
00260    //evaluating the geometric acceptance  
00261    if ( pt1SmearedGaus >= accPtMin_  && pt2SmearedGaus >= accPtMin_ &&  fabs(eta1)>= accEtaMin_  && fabs(eta2) >= accEtaMin_ && fabs(eta1)<= accEtaMax_  && fabs(eta2) <= accEtaMax_ && inv_mass>= accMassMin_ && inv_mass<= accMassMax_) 
00262  nAccPtScaleSmearedGaus_++;
00263           
00264 
00265 
00266   
00267    }}
00268 
00269 
00270   
00271 
00272 void ZMuPtScaleAnalyzer::endJob() {
00273   cout << " number of events accepted :" << nAcc_ << endl;
00274   cout << " number of total events :" << h_mZMC_->GetEntries()  << endl;
00275    cout << " number of cases in which BothMuHasZHasGrandMa :" << nBothMuHasZHasGrandMa_  << endl;
00276   cout << " number of events pt scaled positively accepted :" << nAccPtScaleP_ << endl;
00277 
00278   cout << " number of events pt scaled negatively accepted :" << nAccPtScaleN_ << endl;
00279 
00280 cout << " number of events pt scaled smeared flattely accepted :" << nAccPtScaleSmearedFlat_ << endl;
00281 
00282 cout << " number of events pt scaled smeared gaussianely accepted :" << nAccPtScaleSmearedGaus_ << endl;
00283 
00284 
00285   double eff = (double)nAcc_ / (double) h_mZMC_->GetEntries();
00286   double err = sqrt( eff * (1. - eff) / (double) h_mZMC_->GetEntries() );
00287   cout << " geometric acceptance: " << eff << "+/-" << err << endl; 
00288 
00289  double effScaledP = (double)nAccPtScaleP_ / (double) h_mZMC_->GetEntries();
00290   double errScaledP = sqrt( effScaledP * (1. - effScaledP) / (double) h_mZMC_->GetEntries() );
00291   cout << " geometric acceptance when pt muon is positively scaled: " << effScaledP << "+/-" << errScaledP << endl; 
00292 
00293  double effScaledN = (double)nAccPtScaleN_ / (double) h_mZMC_->GetEntries();
00294   double errScaledN = sqrt( effScaledN * (1. - effScaledN) / (double) h_mZMC_->GetEntries() );
00295   cout << " geometric acceptance when pt muon is negatively scaled: " << effScaledN << "+/-" << errScaledN << endl; 
00296 
00297  double effSmearedFlat = (double) nAccPtScaleSmearedFlat_ / (double) h_mZMC_->GetEntries();
00298   double errSmearedFlat = sqrt( effSmearedFlat * (1. - effSmearedFlat) / (double) h_mZMC_->GetEntries() );
00299   cout << " geometric acceptance when pt muon is scaled with a flat smaering: " << effSmearedFlat << "+/-" << errSmearedFlat << endl; 
00300 
00301  double effSmearedGaus = (double) nAccPtScaleSmearedGaus_ / (double) h_mZMC_->GetEntries();
00302   double errSmearedGaus = sqrt( effSmearedGaus * (1. - effSmearedGaus) / (double) h_mZMC_->GetEntries() );
00303   cout << " geometric acceptance when pt muon is scaled with a gaussian smearing: " << effSmearedGaus << "+/-" << errSmearedGaus << endl; 
00304 
00305 
00306 }
00307 #include "FWCore/Framework/interface/MakerMacros.h"
00308 
00309 DEFINE_FWK_MODULE(ZMuPtScaleAnalyzer);
00310
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/ElectroWeakAnalysis/ZMuMu/plugins/ZMuPtScaleAnalyzer.cc
