Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.6 2009/12/06 14:59:28 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),
  fIsData(kFALSE),
  fPt_histo(0),
  fMCEventInfo(0)
{
  // Constructor
}

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

  delete fPt_histo;
}

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

  Double_t theWeight = 1.0;

  if(fIsData == kFALSE){
    // get the bosons
    MCParticleCol *GenBosons = GetObjThisEvt<MCParticleCol>(fMCBosonsName);

    LoadBranch(fMCEvInfoName);

    // 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::Max(TMath::Min(theWeight,3.999),-3.999));
        }
      }
    } else {
      theWeight = fMCEventInfo->Weight();
      hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
      theWeight = theWeight/fabs(theWeight);
      hDHKFactor[0]->Fill(0.5,theWeight);
    }
    // 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.

  if(fIsData == kFALSE){
    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,-4.0,4.0); 
    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.7
Committed:	Sat Mar 13 20:50:32 2010 UTC (15 years, 1 month ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_014e, Mit_014d, Mit_014c, Mit_014b, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1
Changes since 1.6:	+8 -3 lines
Log Message:	adding more info
#	User	Rev	Content
1	ceballos	1.7	// $Id: HKFactorProducer.cc,v 1.6 2009/12/06 14:59:28 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	ceballos	1.6	fIsData(kFALSE),
23	loizides	1.2	fPt_histo(0),
24	ceballos	1.1	fMCEventInfo(0)
25			{
26	loizides	1.2	// Constructor
27	ceballos	1.1	}
28
29	loizides	1.2	//--------------------------------------------------------------------------------------------------
30	ceballos	1.1	HKFactorProducer::~HKFactorProducer()
31			{
32	loizides	1.2	// Destructor
33
34	ceballos	1.1	delete fPt_histo;
35			}
36
37			//--------------------------------------------------------------------------------------------------
38			void HKFactorProducer::Process()
39			{
40			// Process entries of the tree.
41
42			Double_t theWeight = 1.0;
43	loizides	1.2
44	ceballos	1.6	if(fIsData == kFALSE){
45			// get the bosons
46			MCParticleCol *GenBosons = GetObjThisEvt<MCParticleCol>(fMCBosonsName);
47
48			LoadBranch(fMCEvInfoName);
49
50			// only accept the exact process id
51			if (fProcessID == fMCEventInfo->ProcessId()) {
52
53			Double_t ptH = -1.0;
54			for (UInt_t i=0; i<GenBosons->GetEntries(); ++i) {
55			if(GenBosons->At(i)->PdgId() == MCParticle::kH) {
56			ptH = GenBosons->At(i)->Pt();
57			break;
58			}
59	ceballos	1.1	}
60	loizides	1.2
61	ceballos	1.6	if(ptH >= 0) {
62			// calculate bin size
63			Double_t binsize = fPt_histo->GetXaxis()->GetXmax()/fPt_histo->GetNbinsX();
64			// get bin
65			Int_t bin = Int_t((ptH/binsize)) + 1;
66			// overflow protection: use last entry
67			if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX();
68			theWeight = fPt_histo->GetBinContent(bin);
69
70			if (0) {
71			cout << "Bin Size: " << binsize << ", Higgs Pt: " << ptH
72			<< ", Bin: "<< bin << ", KFactor: "<< fPt_histo->GetBinContent(bin) << endl;
73			}
74
75			if (GetFillHist()) {
76			hDHKFactor[0]->Fill(0.5);
77			hDHKFactor[1]->Fill(TMath::Min(ptH,499.999));
78			hDHKFactor[2]->Fill(TMath::Min(ptH,499.999),theWeight);
79	ceballos	1.7	hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
80	ceballos	1.6	}
81	ceballos	1.1	}
82	ceballos	1.7	} else {
83			theWeight = fMCEventInfo->Weight();
84			hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
85			theWeight = theWeight/fabs(theWeight);
86			hDHKFactor[0]->Fill(0.5,theWeight);
87	ceballos	1.1	}
88	ceballos	1.6	// process id distribution
89			if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));
90	ceballos	1.1	}
91
92			TParameter<Double_t> *NNLOWeight = new TParameter<Double_t>("NNLOWeight", theWeight);
93			AddObjThisEvt(NNLOWeight);
94			}
95
96			//--------------------------------------------------------------------------------------------------
97			void HKFactorProducer::SlaveBegin()
98			{
99			// Book branch and histograms if wanted.
100
101	ceballos	1.6	if(fIsData == kFALSE){
102			ReqBranch(fMCEvInfoName, fMCEventInfo);
103			}
104	loizides	1.2
105			if (!fPt_histo) {
106			Info("SlaveBegin", "Using %s as input data file", fInputFileName.Data());
107			fPt_histo = new HWWKfactorList("KFactorList", fInputFileName);
108			}
109	ceballos	1.1
110	loizides	1.2	if (GetFillHist()) {
111	ceballos	1.1	char sb[1024];
112			sprintf(sb,"hDHKFactor_%d", 0); hDHKFactor[0] = new TH1D(sb,sb,1,0,1);
113			sprintf(sb,"hDHKFactor_%d", 1); hDHKFactor[1] = new TH1D(sb,sb,500,0.,500.);
114			sprintf(sb,"hDHKFactor_%d", 2); hDHKFactor[2] = new TH1D(sb,sb,500,0.,500.);
115	ceballos	1.7	sprintf(sb,"hDHKFactor_%d", 3); hDHKFactor[3] = new TH1D(sb,sb,400,-4.0,4.0);
116	ceballos	1.4	sprintf(sb,"hDHKFactor_%d", 4); hDHKFactor[4] = new TH1D(sb,sb,1000,-0.5,999.5);
117			for(Int_t i=0; i<5; i++) AddOutput(hDHKFactor[i]);
118	ceballos	1.1	}
119			}