Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.13 2011/07/26 13:43:06 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>
#include <TTree.h>
#include <TFile.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),
  fMakePDFNtuple(kFALSE),
  fDoHiggsPtReweighting(kFALSE),
  fPt_histo(0),
  fMCEventInfo(0),
  fOutputFile(0),
  fOutputName("ntuple.root")
{
  // 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->GetXaxis()->GetXmin())/fPt_histo->GetNbinsX();
        Int_t bin = 0;
        // underflow protection: use underflow entry
        if(ptH >= fPt_histo->GetXaxis()->GetXmin()){
          bin = Int_t((ptH-fPt_histo->GetXaxis()->GetXmin())/binsize) + 1;
        }
        // overflow protection: use overflow entry
        if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX()+1;
        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 if (fProcessID == 999){ // for MCatNLO we care about positive or negative weights
      theWeight = fMCEventInfo->Weight();
      if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
      theWeight = theWeight/fabs(theWeight);
      if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
    }
    else if (fProcessID == 998){ // for other samples we care about the actual weights
      theWeight = fMCEventInfo->Weight();
      if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
      if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
    }
    else {
      if (GetFillHist()) hDHKFactor[0]->Fill(0.5,1.0);
    }
    // process id distribution
    if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));

    if (fMakePDFNtuple == kTRUE){
      fTreeVariables[0] = fMCEventInfo->Weight();
      fTreeVariables[1] = fMCEventInfo->Scale();
      fTreeVariables[2] = fMCEventInfo->Id1();
      fTreeVariables[3] = fMCEventInfo->X1();
      fTreeVariables[4] = fMCEventInfo->Pdf1();
      fTreeVariables[5] = fMCEventInfo->Id2();
      fTreeVariables[6] = fMCEventInfo->X2();
      fTreeVariables[7] = fMCEventInfo->Pdf2();
      fTree->Fill();
    }
  }

  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 (fDoHiggsPtReweighting) {
    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]);
  }

  //***********************************************************************************************
  //Create Ntuple Tree  
  //***********************************************************************************************
  if (fMakePDFNtuple == kTRUE){
    printf("... init PDF ntuple ...\n");
    fOutputFile = new TFile(fOutputName, "RECREATE");
    fTree = new TTree("PDFTree", "PDFTree");
    const char* TreeFormat;
    TreeFormat = "weight/F:lq:lid1:lx1:lpdf1:lid2:lx2:lpdf2";
    fTree->Branch("H", &fTreeVariables,TreeFormat);    
  }
}

//--------------------------------------------------------------------------------------------------
void HKFactorProducer::SlaveTerminate(){
  if (fMakePDFNtuple == kTRUE){
    fOutputFile->cd();
    fTree->Write();
    fOutputFile->Close();
 }
}
Revision:	1.14
Committed:	Tue Jul 26 15:01:21 2011 UTC (13 years, 9 months ago) by sixie
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024
Changes since 1.13:	+7 -4 lines
Log Message:	add bool option to turn off loading kfactor
#	User	Rev	Content
1	sixie	1.14	// $Id: HKFactorProducer.cc,v 1.13 2011/07/26 13:43:06 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	ceballos	1.9	#include <TTree.h>
11	ceballos	1.11	#include <TFile.h>
12	ceballos	1.1
13			using namespace mithep;
14
15			ClassImp(mithep::HKFactorProducer)
16
17			//--------------------------------------------------------------------------------------------------
18			HKFactorProducer::HKFactorProducer(const char name, const char title) :
19			BaseMod(name,title),
20			fProcessID(102),
21	loizides	1.2	fInputFileName("setme"),
22	ceballos	1.1	fMCBosonsName(ModNames::gkMCBosonsName),
23	loizides	1.2	fMCEvInfoName(Names::gkMCEvtInfoBrn),
24	ceballos	1.6	fIsData(kFALSE),
25	ceballos	1.9	fMakePDFNtuple(kFALSE),
26	sixie	1.14	fDoHiggsPtReweighting(kFALSE),
27	loizides	1.2	fPt_histo(0),
28	ceballos	1.11	fMCEventInfo(0),
29			fOutputFile(0),
30			fOutputName("ntuple.root")
31	ceballos	1.1	{
32	loizides	1.2	// Constructor
33	ceballos	1.1	}
34
35	loizides	1.2	//--------------------------------------------------------------------------------------------------
36	ceballos	1.1	HKFactorProducer::~HKFactorProducer()
37			{
38	loizides	1.2	// Destructor
39	ceballos	1.1	delete fPt_histo;
40			}
41
42			//--------------------------------------------------------------------------------------------------
43			void HKFactorProducer::Process()
44			{
45			// Process entries of the tree.
46
47			Double_t theWeight = 1.0;
48	loizides	1.2
49	ceballos	1.6	if(fIsData == kFALSE){
50			// get the bosons
51			MCParticleCol *GenBosons = GetObjThisEvt<MCParticleCol>(fMCBosonsName);
52
53			LoadBranch(fMCEvInfoName);
54
55			// only accept the exact process id
56			if (fProcessID == fMCEventInfo->ProcessId()) {
57
58			Double_t ptH = -1.0;
59			for (UInt_t i=0; i<GenBosons->GetEntries(); ++i) {
60			if(GenBosons->At(i)->PdgId() == MCParticle::kH) {
61			ptH = GenBosons->At(i)->Pt();
62			break;
63			}
64	ceballos	1.1	}
65	loizides	1.2
66	ceballos	1.6	if(ptH >= 0) {
67			// calculate bin size
68	ceballos	1.10	Double_t binsize = (fPt_histo->GetXaxis()->GetXmax()-fPt_histo->GetXaxis()->GetXmin())/fPt_histo->GetNbinsX();
69			Int_t bin = 0;
70			// underflow protection: use underflow entry
71			if(ptH >= fPt_histo->GetXaxis()->GetXmin()){
72			bin = Int_t((ptH-fPt_histo->GetXaxis()->GetXmin())/binsize) + 1;
73			}
74			// overflow protection: use overflow entry
75			if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX()+1;
76	ceballos	1.6	theWeight = fPt_histo->GetBinContent(bin);
77
78			if (0) {
79			cout << "Bin Size: " << binsize << ", Higgs Pt: " << ptH
80			<< ", Bin: "<< bin << ", KFactor: "<< fPt_histo->GetBinContent(bin) << endl;
81			}
82
83			if (GetFillHist()) {
84			hDHKFactor[0]->Fill(0.5);
85			hDHKFactor[1]->Fill(TMath::Min(ptH,499.999));
86			hDHKFactor[2]->Fill(TMath::Min(ptH,499.999),theWeight);
87	ceballos	1.7	hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
88	ceballos	1.6	}
89	ceballos	1.1	}
90	ceballos	1.13	}
91			else if (fProcessID == 999){ // for MCatNLO we care about positive or negative weights
92	ceballos	1.7	theWeight = fMCEventInfo->Weight();
93	ceballos	1.8	if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
94	ceballos	1.7	theWeight = theWeight/fabs(theWeight);
95	ceballos	1.8	if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
96			}
97	ceballos	1.13	else if (fProcessID == 998){ // for other samples we care about the actual weights
98			theWeight = fMCEventInfo->Weight();
99			if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
100			if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
101			}
102	ceballos	1.8	else {
103			if (GetFillHist()) hDHKFactor[0]->Fill(0.5,1.0);
104	ceballos	1.1	}
105	ceballos	1.6	// process id distribution
106			if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));
107	ceballos	1.9
108			if (fMakePDFNtuple == kTRUE){
109			fTreeVariables[0] = fMCEventInfo->Weight();
110			fTreeVariables[1] = fMCEventInfo->Scale();
111			fTreeVariables[2] = fMCEventInfo->Id1();
112			fTreeVariables[3] = fMCEventInfo->X1();
113			fTreeVariables[4] = fMCEventInfo->Pdf1();
114			fTreeVariables[5] = fMCEventInfo->Id2();
115			fTreeVariables[6] = fMCEventInfo->X2();
116			fTreeVariables[7] = fMCEventInfo->Pdf2();
117			fTree->Fill();
118			}
119	ceballos	1.1	}
120
121			TParameter<Double_t> *NNLOWeight = new TParameter<Double_t>("NNLOWeight", theWeight);
122			AddObjThisEvt(NNLOWeight);
123			}
124
125			//--------------------------------------------------------------------------------------------------
126			void HKFactorProducer::SlaveBegin()
127			{
128			// Book branch and histograms if wanted.
129
130	ceballos	1.6	if(fIsData == kFALSE){
131			ReqBranch(fMCEvInfoName, fMCEventInfo);
132			}
133	loizides	1.2
134	sixie	1.14	if (fDoHiggsPtReweighting) {
135			if (!fPt_histo) {
136			Info("SlaveBegin", "Using %s as input data file", fInputFileName.Data());
137			fPt_histo = new HWWKfactorList("KFactorList", fInputFileName);
138			}
139	loizides	1.2	}
140	ceballos	1.1
141	loizides	1.2	if (GetFillHist()) {
142	ceballos	1.1	char sb[1024];
143			sprintf(sb,"hDHKFactor_%d", 0); hDHKFactor[0] = new TH1D(sb,sb,1,0,1);
144			sprintf(sb,"hDHKFactor_%d", 1); hDHKFactor[1] = new TH1D(sb,sb,500,0.,500.);
145			sprintf(sb,"hDHKFactor_%d", 2); hDHKFactor[2] = new TH1D(sb,sb,500,0.,500.);
146	ceballos	1.7	sprintf(sb,"hDHKFactor_%d", 3); hDHKFactor[3] = new TH1D(sb,sb,400,-4.0,4.0);
147	ceballos	1.4	sprintf(sb,"hDHKFactor_%d", 4); hDHKFactor[4] = new TH1D(sb,sb,1000,-0.5,999.5);
148			for(Int_t i=0; i<5; i++) AddOutput(hDHKFactor[i]);
149	ceballos	1.1	}
150	ceballos	1.9
151			//***********************************************************************************************
152			//Create Ntuple Tree
153			//***********************************************************************************************
154			if (fMakePDFNtuple == kTRUE){
155			printf("... init PDF ntuple ...\n");
156	ceballos	1.11	fOutputFile = new TFile(fOutputName, "RECREATE");
157	ceballos	1.9	fTree = new TTree("PDFTree", "PDFTree");
158			const char* TreeFormat;
159			TreeFormat = "weight/F:lq:lid1:lx1:lpdf1:lid2:lx2:lpdf2";
160			fTree->Branch("H", &fTreeVariables,TreeFormat);
161			}
162	ceballos	1.1	}
163	ceballos	1.11
164			//--------------------------------------------------------------------------------------------------
165			void HKFactorProducer::SlaveTerminate(){
166	ceballos	1.12	if (fMakePDFNtuple == kTRUE){
167			fOutputFile->cd();
168			fTree->Write();
169			fOutputFile->Close();
170			}
171	ceballos	1.11	}