Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.4 2009/09/03 07:08:33 ceballos Exp $

#include "MitPhysics/Mods/interface/HKFactorProducer.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitAna/DataTree/interface/MCParticleCol.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include <TH1D.h>
#include <TH2D.h>
#include <TParameter.h>

using namespace mithep;

ClassImp(mithep::HKFactorProducer)

//--------------------------------------------------------------------------------------------------
HKFactorProducer::HKFactorProducer(const char *name, const char *title) : 
  BaseMod(name,title),
  fProcessID(102),
  fInputFileName("setme"),
  fMCBosonsName(ModNames::gkMCBosonsName),
  fMCEvInfoName(Names::gkMCEvtInfoBrn),
  fPt_histo(0),
  fMCEventInfo(0)
{
  // Constructor
}

//--------------------------------------------------------------------------------------------------
HKFactorProducer::~HKFactorProducer()
{
  // Destructor

  delete fPt_histo;
}

//--------------------------------------------------------------------------------------------------
void HKFactorProducer::Process()
{
  // Process entries of the tree.

  // get the bosons
  MCParticleCol *GenBosons = GetObjThisEvt<MCParticleCol>(fMCBosonsName);

  LoadBranch(fMCEvInfoName);

  Double_t theWeight = 1.0;

  // only accept the exact process id
  if (fProcessID == fMCEventInfo->ProcessId()) { 

    Double_t ptH = -1.0;
    for (UInt_t i=0; i<GenBosons->GetEntries(); ++i) {
      if(GenBosons->At(i)->PdgId() == MCParticle::kH) {
        ptH = GenBosons->At(i)->Pt();
        break;
      }
    }
  
    if(ptH >= 0) {
      // calculate bin size
      Double_t binsize = fPt_histo->GetXaxis()->GetXmax()/fPt_histo->GetNbinsX();
      // get bin
      Int_t bin = Int_t((ptH/binsize)) + 1;
      // overflow protection: use last entry
      if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX();
      theWeight = fPt_histo->GetBinContent(bin);

      if (0) {
        cout << "Bin Size: " << binsize << ", Higgs Pt: " << ptH
             << ", Bin: "<< bin  << ", KFactor: "<< fPt_histo->GetBinContent(bin) << endl;
      }

      if (GetFillHist()) {
        hDHKFactor[0]->Fill(0.5);
        hDHKFactor[1]->Fill(TMath::Min(ptH,499.999));
        hDHKFactor[2]->Fill(TMath::Min(ptH,499.999),theWeight);
        hDHKFactor[3]->Fill(TMath::Min(theWeight,9.999));
      }
    }
  }
  // process id distribution
  if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));

  TParameter<Double_t> *NNLOWeight = new TParameter<Double_t>("NNLOWeight", theWeight);
  AddObjThisEvt(NNLOWeight);
}

//--------------------------------------------------------------------------------------------------
void HKFactorProducer::SlaveBegin()
{
  // Book branch and histograms if wanted.

  ReqBranch(fMCEvInfoName, fMCEventInfo);

  if (!fPt_histo) {
    Info("SlaveBegin", "Using %s as input data file", fInputFileName.Data());
    fPt_histo = new  HWWKfactorList("KFactorList", fInputFileName);
  }

  if (GetFillHist()) {
    char sb[1024];
    sprintf(sb,"hDHKFactor_%d", 0);  hDHKFactor[0]  = new TH1D(sb,sb,1,0,1); 
    sprintf(sb,"hDHKFactor_%d", 1);  hDHKFactor[1]  = new TH1D(sb,sb,500,0.,500.); 
    sprintf(sb,"hDHKFactor_%d", 2);  hDHKFactor[2]  = new TH1D(sb,sb,500,0.,500.); 
    sprintf(sb,"hDHKFactor_%d", 3);  hDHKFactor[3]  = new TH1D(sb,sb,400,0.,4.); 
    sprintf(sb,"hDHKFactor_%d", 4);  hDHKFactor[4]  = new TH1D(sb,sb,1000,-0.5,999.5); 
    for(Int_t i=0; i<5; i++) AddOutput(hDHKFactor[i]);
  }
}
Revision:	1.5
Committed:	Thu Sep 3 07:21:52 2009 UTC (15 years, 8 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012, Mit_011a
Changes since 1.4:	+2 -2 lines
Log Message:	new process id plot
#	User	Rev	Content
1	ceballos	1.5	// $Id: HKFactorProducer.cc,v 1.4 2009/09/03 07:08:33 ceballos Exp $
2	ceballos	1.1
3			#include "MitPhysics/Mods/interface/HKFactorProducer.h"
4			#include "MitCommon/MathTools/interface/MathUtils.h"
5	loizides	1.3	#include "MitAna/DataTree/interface/MCParticleCol.h"
6	ceballos	1.1	#include "MitPhysics/Init/interface/ModNames.h"
7			#include <TH1D.h>
8			#include <TH2D.h>
9			#include <TParameter.h>
10
11			using namespace mithep;
12
13			ClassImp(mithep::HKFactorProducer)
14
15			//--------------------------------------------------------------------------------------------------
16			HKFactorProducer::HKFactorProducer(const char name, const char title) :
17			BaseMod(name,title),
18			fProcessID(102),
19	loizides	1.2	fInputFileName("setme"),
20	ceballos	1.1	fMCBosonsName(ModNames::gkMCBosonsName),
21	loizides	1.2	fMCEvInfoName(Names::gkMCEvtInfoBrn),
22			fPt_histo(0),
23	ceballos	1.1	fMCEventInfo(0)
24			{
25	loizides	1.2	// Constructor
26	ceballos	1.1	}
27
28	loizides	1.2	//--------------------------------------------------------------------------------------------------
29	ceballos	1.1	HKFactorProducer::~HKFactorProducer()
30			{
31	loizides	1.2	// Destructor
32
33	ceballos	1.1	delete fPt_histo;
34			}
35
36			//--------------------------------------------------------------------------------------------------
37			void HKFactorProducer::Process()
38			{
39			// Process entries of the tree.
40
41	loizides	1.2	// get the bosons
42			MCParticleCol *GenBosons = GetObjThisEvt<MCParticleCol>(fMCBosonsName);
43
44			LoadBranch(fMCEvInfoName);
45	ceballos	1.1
46			Double_t theWeight = 1.0;
47	loizides	1.2
48			// only accept the exact process id
49			if (fProcessID == fMCEventInfo->ProcessId()) {
50
51	ceballos	1.1	Double_t ptH = -1.0;
52			for (UInt_t i=0; i<GenBosons->GetEntries(); ++i) {
53			if(GenBosons->At(i)->PdgId() == MCParticle::kH) {
54			ptH = GenBosons->At(i)->Pt();
55			break;
56			}
57			}
58
59	loizides	1.2	if(ptH >= 0) {
60			// calculate bin size
61	ceballos	1.1	Double_t binsize = fPt_histo->GetXaxis()->GetXmax()/fPt_histo->GetNbinsX();
62	loizides	1.2	// get bin
63	ceballos	1.1	Int_t bin = Int_t((ptH/binsize)) + 1;
64			// overflow protection: use last entry
65			if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX();
66			theWeight = fPt_histo->GetBinContent(bin);
67
68	loizides	1.2	if (0) {
69			cout << "Bin Size: " << binsize << ", Higgs Pt: " << ptH
70			<< ", Bin: "<< bin << ", KFactor: "<< fPt_histo->GetBinContent(bin) << endl;
71	ceballos	1.1	}
72	loizides	1.2
73			if (GetFillHist()) {
74	ceballos	1.1	hDHKFactor[0]->Fill(0.5);
75			hDHKFactor[1]->Fill(TMath::Min(ptH,499.999));
76			hDHKFactor[2]->Fill(TMath::Min(ptH,499.999),theWeight);
77			hDHKFactor[3]->Fill(TMath::Min(theWeight,9.999));
78			}
79			}
80			}
81	ceballos	1.4	// process id distribution
82	ceballos	1.5	if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));
83	ceballos	1.1
84			TParameter<Double_t> *NNLOWeight = new TParameter<Double_t>("NNLOWeight", theWeight);
85			AddObjThisEvt(NNLOWeight);
86			}
87
88			//--------------------------------------------------------------------------------------------------
89			void HKFactorProducer::SlaveBegin()
90			{
91			// Book branch and histograms if wanted.
92
93	loizides	1.2	ReqBranch(fMCEvInfoName, fMCEventInfo);
94
95			if (!fPt_histo) {
96			Info("SlaveBegin", "Using %s as input data file", fInputFileName.Data());
97			fPt_histo = new HWWKfactorList("KFactorList", fInputFileName);
98			}
99	ceballos	1.1
100	loizides	1.2	if (GetFillHist()) {
101	ceballos	1.1	char sb[1024];
102			sprintf(sb,"hDHKFactor_%d", 0); hDHKFactor[0] = new TH1D(sb,sb,1,0,1);
103			sprintf(sb,"hDHKFactor_%d", 1); hDHKFactor[1] = new TH1D(sb,sb,500,0.,500.);
104			sprintf(sb,"hDHKFactor_%d", 2); hDHKFactor[2] = new TH1D(sb,sb,500,0.,500.);
105			sprintf(sb,"hDHKFactor_%d", 3); hDHKFactor[3] = new TH1D(sb,sb,400,0.,4.);
106	ceballos	1.4	sprintf(sb,"hDHKFactor_%d", 4); hDHKFactor[4] = new TH1D(sb,sb,1000,-0.5,999.5);
107			for(Int_t i=0; i<5; i++) AddOutput(hDHKFactor[i]);
108	ceballos	1.1	}
109			}