Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.14 2011/07/26 15:01:21 sixie 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),
  fEmbedWeightName(Names::gkEmbedWeightBrn),
  fIsData(kFALSE),
  fMakePDFNtuple(kFALSE),
  fDoHiggsPtReweighting(kFALSE),
  fPt_histo(0),
  fMCEventInfo(0),
  fEmbedWeight(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();
    }
  }
  else if (fProcessID == 997){ // for tau embedding samples
    LoadBranch(fEmbedWeightName);
    theWeight = fEmbedWeight->At(fEmbedWeight->GetEntries()-1)->Weight();
    if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
    if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
  }

  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(fProcessID == 997){
    ReqBranch(fEmbedWeightName, fEmbedWeight);
  }

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