Mods/src/HKFactorProducer.cc

// $Id: HKFactorProducer.cc,v 1.8 2010/09/27 15:49:24 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->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 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.9
Committed:	Tue Oct 19 22:25:25 2010 UTC (14 years, 6 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, Mit_016, Mit_015b, Mit_015a
Changes since 1.8:	+27 -1 lines
Log Message:	new stuff for PDFs
#	User	Rev	Content
1	ceballos	1.9	// $Id: HKFactorProducer.cc,v 1.8 2010/09/27 15:49:24 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			Double_t binsize = fPt_histo->GetXaxis()->GetXmax()/fPt_histo->GetNbinsX();
66			// get bin
67			Int_t bin = Int_t((ptH/binsize)) + 1;
68			// overflow protection: use last entry
69			if(bin > fPt_histo->GetNbinsX()) bin=fPt_histo->GetNbinsX();
70			theWeight = fPt_histo->GetBinContent(bin);
71
72			if (0) {
73			cout << "Bin Size: " << binsize << ", Higgs Pt: " << ptH
74			<< ", Bin: "<< bin << ", KFactor: "<< fPt_histo->GetBinContent(bin) << endl;
75			}
76
77			if (GetFillHist()) {
78			hDHKFactor[0]->Fill(0.5);
79			hDHKFactor[1]->Fill(TMath::Min(ptH,499.999));
80			hDHKFactor[2]->Fill(TMath::Min(ptH,499.999),theWeight);
81	ceballos	1.7	hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
82	ceballos	1.6	}
83	ceballos	1.1	}
84	ceballos	1.8	} else if (fProcessID == 999){
85	ceballos	1.7	theWeight = fMCEventInfo->Weight();
86	ceballos	1.8	if (GetFillHist()) hDHKFactor[3]->Fill(TMath::Max(TMath::Min(theWeight,3.999),-3.999));
87	ceballos	1.7	theWeight = theWeight/fabs(theWeight);
88	ceballos	1.8	if (GetFillHist()) hDHKFactor[0]->Fill(0.5,theWeight);
89			}
90			else {
91			if (GetFillHist()) hDHKFactor[0]->Fill(0.5,1.0);
92	ceballos	1.1	}
93	ceballos	1.6	// process id distribution
94			if (GetFillHist()) hDHKFactor[4]->Fill(TMath::Min((Double_t)fMCEventInfo->ProcessId(),999.499));
95	ceballos	1.9
96			if (fMakePDFNtuple == kTRUE){
97			fTreeVariables[0] = fMCEventInfo->Weight();
98			fTreeVariables[1] = fMCEventInfo->Scale();
99			fTreeVariables[2] = fMCEventInfo->Id1();
100			fTreeVariables[3] = fMCEventInfo->X1();
101			fTreeVariables[4] = fMCEventInfo->Pdf1();
102			fTreeVariables[5] = fMCEventInfo->Id2();
103			fTreeVariables[6] = fMCEventInfo->X2();
104			fTreeVariables[7] = fMCEventInfo->Pdf2();
105			fTree->Fill();
106			}
107	ceballos	1.1	}
108
109			TParameter<Double_t> *NNLOWeight = new TParameter<Double_t>("NNLOWeight", theWeight);
110			AddObjThisEvt(NNLOWeight);
111			}
112
113			//--------------------------------------------------------------------------------------------------
114			void HKFactorProducer::SlaveBegin()
115			{
116			// Book branch and histograms if wanted.
117
118	ceballos	1.6	if(fIsData == kFALSE){
119			ReqBranch(fMCEvInfoName, fMCEventInfo);
120			}
121	loizides	1.2
122			if (!fPt_histo) {
123			Info("SlaveBegin", "Using %s as input data file", fInputFileName.Data());
124			fPt_histo = new HWWKfactorList("KFactorList", fInputFileName);
125			}
126	ceballos	1.1
127	loizides	1.2	if (GetFillHist()) {
128	ceballos	1.1	char sb[1024];
129			sprintf(sb,"hDHKFactor_%d", 0); hDHKFactor[0] = new TH1D(sb,sb,1,0,1);
130			sprintf(sb,"hDHKFactor_%d", 1); hDHKFactor[1] = new TH1D(sb,sb,500,0.,500.);
131			sprintf(sb,"hDHKFactor_%d", 2); hDHKFactor[2] = new TH1D(sb,sb,500,0.,500.);
132	ceballos	1.7	sprintf(sb,"hDHKFactor_%d", 3); hDHKFactor[3] = new TH1D(sb,sb,400,-4.0,4.0);
133	ceballos	1.4	sprintf(sb,"hDHKFactor_%d", 4); hDHKFactor[4] = new TH1D(sb,sb,1000,-0.5,999.5);
134			for(Int_t i=0; i<5; i++) AddOutput(hDHKFactor[i]);
135	ceballos	1.1	}
136	ceballos	1.9
137			//***********************************************************************************************
138			//Create Ntuple Tree
139			//***********************************************************************************************
140			if (fMakePDFNtuple == kTRUE){
141			printf("... init PDF ntuple ...\n");
142			fTree = new TTree("PDFTree", "PDFTree");
143			const char* TreeFormat;
144			TreeFormat = "weight/F:lq:lid1:lx1:lpdf1:lid2:lx2:lpdf2";
145			fTree->Branch("H", &fTreeVariables,TreeFormat);
146			AddOutput(fTree);
147			}
148	ceballos	1.1	}