Inheritance diagram for JetCorrectorDemo:

Public Types
typedef reco::Particle::LorentzVector	LorentzVector
Public Member Functions
	JetCorrectorDemo (const edm::ParameterSet &)
	~JetCorrectorDemo ()
Private Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
virtual void	endJob ()
Private Attributes
edm::Service< TFileService >	fs
bool	mDebug
double	mEtaMax
double	mEtaMin
TH2F *	mJECvsEta
TH2F *	mJECvsPt
std::string	mJetCorService
int	mNGraphPoints
int	mNHistoPoints
std::string	mPayloadName
double	mPtMax
double	mPtMin
TRandom *	mRandom
TGraph *	mUncCorPt [100]
std::string	mUncertaintyFile
std::string	mUncertaintyTag
TGraph *	mUncEta [100]
bool	mUseCondDB
std::vector< double >	mVEta
TGraphErrors *	mVGraphCorPt [100]
TGraphErrors *	mVGraphEta [100]
TGraphErrors *	mVGraphPt [100]
std::vector< double >	mVPt
double	veta [100][1000]
double	vex_eta [100][1000]
double	vex_pt [100][1000]
double	vjec_eta [100][1000]
double	vjec_pt [100][1000]
double	vjecUnc_eta [100][1000]
double	vjecUnc_pt [100][1000]
double	vpt [100][1000]
double	vptcor [100][1000]
double	vUnc_eta [100][1000]
double	vUnc_pt [100][1000]

Detailed Description

Definition at line 28 of file JetCorrectorDemo.cc.

Member Typedef Documentation

typedef reco::Particle::LorentzVector JetCorrectorDemo::LorentzVector

Definition at line 32 of file JetCorrectorDemo.cc.

Constructor & Destructor Documentation

JetCorrectorDemo::JetCorrectorDemo ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 56 of file JetCorrectorDemo.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), mDebug, mEtaMax, mEtaMin, mJetCorService, mNGraphPoints, mNHistoPoints, mPayloadName, mPtMax, mPtMin, mUncertaintyFile, mUncertaintyTag, mUseCondDB, mVEta, and mVPt.

{
  mJetCorService     = iConfig.getParameter<std::string>          ("JetCorrectionService");
  mPayloadName       = iConfig.getParameter<std::string>          ("PayloadName");
  mUncertaintyTag    = iConfig.getParameter<std::string>          ("UncertaintyTag");
  mUncertaintyFile   = iConfig.getParameter<std::string>          ("UncertaintyFile");
  mNHistoPoints      = iConfig.getParameter<int>                  ("NHistoPoints");
  mNGraphPoints      = iConfig.getParameter<int>                  ("NGraphPoints");
  mEtaMin            = iConfig.getParameter<double>               ("EtaMin");
  mEtaMax            = iConfig.getParameter<double>               ("EtaMax");
  mPtMin             = iConfig.getParameter<double>               ("PtMin");
  mPtMax             = iConfig.getParameter<double>               ("PtMax");
  mVEta              = iConfig.getParameter<std::vector<double> > ("VEta");
  mVPt               = iConfig.getParameter<std::vector<double> > ("VPt");
  mDebug             = iConfig.getUntrackedParameter<bool>        ("Debug",false);
  mUseCondDB         = iConfig.getUntrackedParameter<bool>        ("UseCondDB",false);
}

JetCorrectorDemo::~JetCorrectorDemo ( )

Definition at line 74 of file JetCorrectorDemo.cc.

{
  
}

Member Function Documentation

void JetCorrectorDemo::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[private, virtual]`

Implements edm::EDAnalyzer.

Definition at line 79 of file JetCorrectorDemo.cc.

References JetCorrector::correction(), gather_cfg::cout, alignCSCRings::e, eta(), edm::FileInPath::fullPath(), edm::EventSetup::get(), JetCorrector::getJetCorrector(), JetCorrectionUncertainty::getUncertainty(), i, patTestJEC_cfi::jec, mDebug, mEtaMax, mEtaMin, mJECvsEta, mJECvsPt, mJetCorService, mNGraphPoints, mNHistoPoints, mPayloadName, mPtMax, mPtMin, mRandom, mUncertaintyFile, mUncertaintyTag, mUseCondDB, mVEta, mVPt, funct::pow(), JetCorrectionUncertainty::setJetEta(), JetCorrectionUncertainty::setJetPt(), tmp, veta, vex_eta, vex_pt, vjec_eta, vjec_pt, vjecUnc_eta, vjecUnc_pt, vpt, vptcor, vUnc_eta, and vUnc_pt.

{
  const JetCorrector* corrector = JetCorrector::getJetCorrector(mJetCorService,iSetup);
  JetCorrectionUncertainty *jecUnc(0);
  if (mUncertaintyTag != "")
    {
      if (mUseCondDB)
        {
          edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
          iSetup.get<JetCorrectionsRecord>().get(mPayloadName,JetCorParColl); 
          JetCorrectorParameters const & JetCorPar = (*JetCorParColl)[mUncertaintyTag];
          jecUnc = new JetCorrectionUncertainty(JetCorPar);
          std::cout<<"Configured Uncertainty from CondDB"<<std::endl;   
        }
      else
        {  
          edm::FileInPath fip("CondFormats/JetMETObjects/data/"+mUncertaintyFile+".txt");
          jecUnc = new JetCorrectionUncertainty(fip.fullPath()); 
        }
    }
  double jec,rawPt,corPt,eta,unc;
  TLorentzVector P4;
  double dEta = (mEtaMax-mEtaMin)/mNGraphPoints;
  for(int i=0;i<mNHistoPoints;i++)
    {
      rawPt  = mRandom->Uniform(mPtMin,mPtMax);
      eta = mRandom->Uniform(mEtaMin,mEtaMax);
      P4.SetPtEtaPhiE(rawPt,eta,0,0); 
      LorentzVector rawP4(P4.Px(),P4.Py(),P4.Pz(),P4.E());
      jec = corrector->correction(rawP4);
      mJECvsEta->Fill(eta,jec);
      mJECvsPt->Fill(rawPt,jec);
    }
  //--------- Pt Graphs ------------------
  for(unsigned ieta=0;ieta<mVEta.size();ieta++)
    {
      double rPt  = pow((3500./TMath::CosH(mVEta[ieta]))/mPtMin,1./mNGraphPoints);
      for(int i=0;i<mNGraphPoints;i++) 
        {
          rawPt  = mPtMin*pow(rPt,i);
          eta = mVEta[ieta];
          vpt[ieta][i] = rawPt;
          P4.SetPtEtaPhiE(rawPt,eta,0,0); 
          LorentzVector rawP4(P4.Px(),P4.Py(),P4.Pz(),P4.E());
          jec = corrector->correction(rawP4);// the jec is a function of the raw pt
          vjec_eta[ieta][i] = jec;
          vptcor[ieta][i] = rawPt*jec;
          unc = 0.0;
          if (mUncertaintyTag != "")
            {
              jecUnc->setJetEta(eta);
              jecUnc->setJetPt(rawPt*jec);// the uncertainty is a function of the corrected pt
              unc = jecUnc->getUncertainty(true);
            }
          vjecUnc_eta[ieta][i] = unc*jec;
          vUnc_eta[ieta][i] = unc;
          vex_eta[ieta][i] = 0.0;
          if (mDebug)
            std::cout<<rawPt<<" "<<eta<<" "<<jec<<" "<<rawPt*jec<<" "<<unc<<std::endl;
        }
    }
  //--------- Eta Graphs -------------------
  for(unsigned ipt=0;ipt<mVPt.size();ipt++)
    {
      for(int i=0;i<mNGraphPoints;i++)
        {
          eta = mEtaMin + i*dEta;
          corPt  = mVPt[ipt];
          veta[ipt][i] = eta;
          //---------- find the raw pt -----------
          double e = 1.0;
          int nLoop(0);
          rawPt = corPt; 
          while(e > 0.0001 && nLoop < 10) 
             {
               P4.SetPtEtaPhiE(rawPt,eta,0,0); 
               LorentzVector rawP4(P4.Px(),P4.Py(),P4.Pz(),P4.E()); 
               jec = corrector->correction(rawP4); 
               double tmp = rawPt * jec;
               e = fabs(tmp-corPt)/corPt;
               if (jec > 0)
                 rawPt = corPt/jec;
               nLoop++;
             } 
          //--------- calculate the jec for the rawPt --------
          P4.SetPtEtaPhiE(rawPt,eta,0,0); 
          LorentzVector rawP4(P4.Px(),P4.Py(),P4.Pz(),P4.E());
          jec = corrector->correction(rawP4);// the jec is a function of the raw pt
          vjec_pt[ipt][i] = jec;
          unc = 0.0;
          if (mUncertaintyTag != "")
            {
              jecUnc->setJetEta(eta);
              jecUnc->setJetPt(corPt);// the uncertainty is a function of the corrected pt
              unc = jecUnc->getUncertainty(true);
            }
          vjecUnc_pt[ipt][i] = unc*jec;
          vUnc_pt[ipt][i] = unc;
          vex_pt[ipt][i] = 0.0;
          if (mDebug)
            std::cout<<rawPt<<" "<<eta<<" "<<jec<<" "<<rawPt*jec<<" "<<unc<<std::endl;
        }
    }

}

void JetCorrectorDemo::beginJob ( void ) [private, virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 185 of file JetCorrectorDemo.cc.

References Exception, fs, mEtaMax, mEtaMin, mJECvsEta, mJECvsPt, mNGraphPoints, mPtMax, mPtMin, mRandom, mVEta, and mVPt.

{
  if (mNGraphPoints > 1000)
    throw  cms::Exception("JetCorrectorDemo","Too many graph points !!! Maximum is 1000 !!!");
  if (mVEta.size() > 100)
    throw  cms::Exception("JetCorrectorDemo","Too many eta values !!! Maximum is 100 !!!");
  if (mVPt.size() > 100)
    throw  cms::Exception("JetCorrectorDemo","Too many pt values !!! Maximum is 100 !!!");
  mJECvsEta = fs->make<TH2F>("JECvsEta","JECvsEta",200,mEtaMin,mEtaMax,100,0,5);
  mJECvsPt  = fs->make<TH2F>("JECvsPt","JECvsPt",200,mPtMin,mPtMax,100,0,5);
  mRandom   = new TRandom();
  mRandom->SetSeed(0);
}

void JetCorrectorDemo::endJob ( void ) [private, virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 199 of file JetCorrectorDemo.cc.

References fs, mNGraphPoints, mUncCorPt, mUncEta, mVEta, mVGraphCorPt, mVGraphEta, mVGraphPt, mVPt, mergeVDriftHistosByStation::name, veta, vex_eta, vex_pt, vjec_eta, vjec_pt, vjecUnc_eta, vjecUnc_pt, vpt, vptcor, vUnc_eta, and vUnc_pt.

{
  char name[1000];
  for(unsigned ipt=0;ipt<mVPt.size();ipt++)
    {
      mVGraphEta[ipt] = fs->make<TGraphErrors>(mNGraphPoints,veta[ipt],vjec_pt[ipt],vex_pt[ipt],vjecUnc_pt[ipt]);
      sprintf(name,"JEC_vs_Eta_CorPt%1.1f",mVPt[ipt]);
      mVGraphEta[ipt]->SetName(name);
      mUncEta[ipt] = fs->make<TGraph>(mNGraphPoints,veta[ipt],vUnc_pt[ipt]);
      sprintf(name,"UNC_vs_Eta_CorPt%1.1f",mVPt[ipt]);
      mUncEta[ipt]->SetName(name);
    }
  for(unsigned ieta=0;ieta<mVEta.size();ieta++)
    {
      mVGraphPt[ieta] = fs->make<TGraphErrors>(mNGraphPoints,vpt[ieta],vjec_eta[ieta],vex_eta[ieta],vjecUnc_eta[ieta]);
      sprintf(name,"JEC_vs_RawPt_eta%1.1f",mVEta[ieta]);
      mVGraphPt[ieta]->SetName(name);
      mVGraphCorPt[ieta] = fs->make<TGraphErrors>(mNGraphPoints,vptcor[ieta],vjec_eta[ieta],vex_eta[ieta],vjecUnc_eta[ieta]);
      sprintf(name,"JEC_vs_CorPt_eta%1.1f",mVEta[ieta]);
      mVGraphCorPt[ieta]->SetName(name);
      mUncCorPt[ieta] = fs->make<TGraph>(mNGraphPoints,vptcor[ieta],vUnc_eta[ieta]);
      sprintf(name,"UNC_vs_CorPt_eta%1.1f",mVEta[ieta]);
      mUncCorPt[ieta]->SetName(name);
    }
}