Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.9 2010/10/19 22:25:25 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>

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