TreeFiller/src/FillerMCEventInfo.cc

// $Id: FillerMCEventInfo.cc,v 1.16 2010/03/18 20:21:00 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerMCEventInfo.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataTree/interface/MCEventInfo.h"
#include "MitProd/ObjectService/interface/ObjectService.h"
#include "DataFormats/HepMCCandidate/interface/FlavorHistory.h"
#include "DataFormats/HepMCCandidate/interface/FlavorHistoryEvent.h"
#include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h"

using namespace std;
using namespace edm;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
FillerMCEventInfo::FillerMCEventInfo(const ParameterSet &cfg, const char *name,  bool active) : 
  BaseFiller(cfg,"MCEventInfo",active),
  evtName_(Conf().getUntrackedParameter<string>("evtName",Names::gkMCEvtInfoBrn)),
  genHepMCEvName_(Conf().getUntrackedParameter<string>("genHepMCEventEdmName","generator")),
  genEvtInfoName_(Conf().getUntrackedParameter<string>("genEvtInfoEdmName","generator")),
  flavorHistoryActive_(Conf().getUntrackedParameter<bool>("flavorHistoryActive",false)),
  flavorHistName_(Conf().getUntrackedParameter<string>("flavorHistEdmName","flavorHistoryFilter")),
  eventInfo_(new MCEventInfo())
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
FillerMCEventInfo::~FillerMCEventInfo()
{
  // Destructor.

  delete eventInfo_;
}

//--------------------------------------------------------------------------------------------------
void FillerMCEventInfo::BookDataBlock(TreeWriter &tws)
{
  // Book our branch.

  tws.AddBranch(evtName_,&eventInfo_);
  OS()->add<mithep::MCEventInfo>(eventInfo_,evtName_);
}

//--------------------------------------------------------------------------------------------------
void FillerMCEventInfo::FillDataBlock(const edm::Event &event, 
                                      const edm::EventSetup &setup)
{
  // Fill our data structures.

  if (event.isRealData()) {
    PrintErrorAndExit("Expected monte-carlo record, but did not get it. Aborting.");
  }
  
  Handle<GenEventInfoProduct> hEvtInfo;

  if (genEvtInfoName_.empty() || !GetProductSafe(genEvtInfoName_, hEvtInfo, event)) {

    // fall back to hepmc if requested
    if (!genHepMCEvName_.empty()) {
      Handle<edm::HepMCProduct> hHepMCProduct;
      GetProduct(genHepMCEvName_, hHepMCProduct, event);
      const HepMC::GenEvent *genEvt = hHepMCProduct->GetEvent();
      eventInfo_->SetScale(genEvt->event_scale());
      eventInfo_->SetProcessId(genEvt->signal_process_id());
      HepMC::WeightContainer wc = genEvt->weights();
      Double_t weight = 0;
      for (unsigned int i = 0; i< wc.size(); ++i) 
        weight *= wc[i];
      eventInfo_->SetWeight(weight);
      const HepMC::PdfInfo *genPdfInfo = genEvt->pdf_info();
      if (genPdfInfo) {
        eventInfo_->SetId1(genPdfInfo->id1());
        eventInfo_->SetId2(genPdfInfo->id2());
        eventInfo_->SetPdf1(genPdfInfo->pdf1());
        eventInfo_->SetPdf2(genPdfInfo->pdf2());
        eventInfo_->SetScalePdf(genPdfInfo->scalePDF());
        eventInfo_->SetX1(genPdfInfo->x1());
        eventInfo_->SetX2(genPdfInfo->x2());
      }
    }

  } else {

    // use event info product
    eventInfo_->SetScale(hEvtInfo->qScale());
    eventInfo_->SetProcessId(hEvtInfo->signalProcessID());
    eventInfo_->SetWeight(hEvtInfo->weight());
    if (hEvtInfo->hasPDF()) {
      const gen::PdfInfo *pdf = hEvtInfo->pdf();
      eventInfo_->SetId1(pdf->id.first);
      eventInfo_->SetId2(pdf->id.second);
      eventInfo_->SetPdf1(pdf->xPDF.first);
      eventInfo_->SetPdf2(pdf->xPDF.second);
      eventInfo_->SetScalePdf(pdf->scalePDF);
      eventInfo_->SetX1(pdf->x.first);
      eventInfo_->SetX2(pdf->x.second);
    }
  }

  // fill flavor history path if requested
  if (flavorHistoryActive_) {
    Handle<unsigned int> flavorHistoryPath;
    GetProduct(flavorHistName_, flavorHistoryPath, event);
    eventInfo_->SetFlavorHistoryPath(*flavorHistoryPath);
  }
}
Revision:	1.17
Committed:	Fri Mar 30 01:08:40 2012 UTC (13 years, 1 month ago) by paus
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_028, Mit_027a, Mit_027, Mit_026, HEAD
Changes since 1.16:	+2 -2 lines
Log Message:	Initial 5 version.
#	User	Rev	Content
1	paus	1.17	// $Id: FillerMCEventInfo.cc,v 1.16 2010/03/18 20:21:00 bendavid Exp $
2	loizides	1.1
3			#include "MitProd/TreeFiller/interface/FillerMCEventInfo.h"
4	bendavid	1.7	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
5	loizides	1.12	#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
6	loizides	1.1	#include "MitAna/DataTree/interface/Names.h"
7			#include "MitAna/DataTree/interface/MCEventInfo.h"
8	loizides	1.6	#include "MitProd/ObjectService/interface/ObjectService.h"
9	phedex	1.8	#include "DataFormats/HepMCCandidate/interface/FlavorHistory.h"
10			#include "DataFormats/HepMCCandidate/interface/FlavorHistoryEvent.h"
11			#include "PhysicsTools/HepMCCandAlgos/interface/FlavorHistoryProducer.h"
12	loizides	1.1
13			using namespace std;
14			using namespace edm;
15			using namespace mithep;
16
17			//--------------------------------------------------------------------------------------------------
18	loizides	1.3	FillerMCEventInfo::FillerMCEventInfo(const ParameterSet &cfg, const char *name, bool active) :
19	loizides	1.1	BaseFiller(cfg,"MCEventInfo",active),
20	loizides	1.12	evtName_(Conf().getUntrackedParameter<string>("evtName",Names::gkMCEvtInfoBrn)),
21			genHepMCEvName_(Conf().getUntrackedParameter<string>("genHepMCEventEdmName","generator")),
22			genEvtInfoName_(Conf().getUntrackedParameter<string>("genEvtInfoEdmName","generator")),
23	sixie	1.9	flavorHistoryActive_(Conf().getUntrackedParameter<bool>("flavorHistoryActive",false)),
24	loizides	1.12	flavorHistName_(Conf().getUntrackedParameter<string>("flavorHistEdmName","flavorHistoryFilter")),
25	loizides	1.1	eventInfo_(new MCEventInfo())
26			{
27			// Constructor.
28			}
29
30			//--------------------------------------------------------------------------------------------------
31			FillerMCEventInfo::~FillerMCEventInfo()
32			{
33			// Destructor.
34
35			delete eventInfo_;
36			}
37
38			//--------------------------------------------------------------------------------------------------
39	bendavid	1.16	void FillerMCEventInfo::BookDataBlock(TreeWriter &tws)
40	loizides	1.1	{
41	loizides	1.14	// Book our branch.
42	loizides	1.1
43	loizides	1.2	tws.AddBranch(evtName_,&eventInfo_);
44			OS()->add<mithep::MCEventInfo>(eventInfo_,evtName_);
45	loizides	1.1	}
46
47			//--------------------------------------------------------------------------------------------------
48			void FillerMCEventInfo::FillDataBlock(const edm::Event &event,
49			const edm::EventSetup &setup)
50			{
51			// Fill our data structures.
52
53	loizides	1.3	if (event.isRealData()) {
54			PrintErrorAndExit("Expected monte-carlo record, but did not get it. Aborting.");
55			}
56	loizides	1.12
57			Handle<GenEventInfoProduct> hEvtInfo;
58	loizides	1.3
59	loizides	1.12	if (genEvtInfoName_.empty() \|\| !GetProductSafe(genEvtInfoName_, hEvtInfo, event)) {
60	loizides	1.4
61	loizides	1.12	// fall back to hepmc if requested
62			if (!genHepMCEvName_.empty()) {
63			Handle<edm::HepMCProduct> hHepMCProduct;
64			GetProduct(genHepMCEvName_, hHepMCProduct, event);
65			const HepMC::GenEvent *genEvt = hHepMCProduct->GetEvent();
66			eventInfo_->SetScale(genEvt->event_scale());
67			eventInfo_->SetProcessId(genEvt->signal_process_id());
68			HepMC::WeightContainer wc = genEvt->weights();
69			Double_t weight = 0;
70	paus	1.17	for (unsigned int i = 0; i< wc.size(); ++i)
71	loizides	1.12	weight *= wc[i];
72			eventInfo_->SetWeight(weight);
73			const HepMC::PdfInfo *genPdfInfo = genEvt->pdf_info();
74	loizides	1.15	if (genPdfInfo) {
75			eventInfo_->SetId1(genPdfInfo->id1());
76			eventInfo_->SetId2(genPdfInfo->id2());
77			eventInfo_->SetPdf1(genPdfInfo->pdf1());
78			eventInfo_->SetPdf2(genPdfInfo->pdf2());
79			eventInfo_->SetScalePdf(genPdfInfo->scalePDF());
80			eventInfo_->SetX1(genPdfInfo->x1());
81			eventInfo_->SetX2(genPdfInfo->x2());
82			}
83	loizides	1.5	}
84
85	loizides	1.12	} else {
86	loizides	1.5
87	loizides	1.12	// use event info product
88			eventInfo_->SetScale(hEvtInfo->qScale());
89			eventInfo_->SetProcessId(hEvtInfo->signalProcessID());
90			eventInfo_->SetWeight(hEvtInfo->weight());
91			if (hEvtInfo->hasPDF()) {
92			const gen::PdfInfo *pdf = hEvtInfo->pdf();
93			eventInfo_->SetId1(pdf->id.first);
94			eventInfo_->SetId2(pdf->id.second);
95			eventInfo_->SetPdf1(pdf->xPDF.first);
96			eventInfo_->SetPdf2(pdf->xPDF.second);
97			eventInfo_->SetScalePdf(pdf->scalePDF);
98			eventInfo_->SetX1(pdf->x.first);
99			eventInfo_->SetX2(pdf->x.second);
100	loizides	1.5	}
101	loizides	1.4	}
102	phedex	1.8
103	loizides	1.12	// fill flavor history path if requested
104	sixie	1.9	if (flavorHistoryActive_) {
105	loizides	1.12	Handle<unsigned int> flavorHistoryPath;
106			GetProduct(flavorHistName_, flavorHistoryPath, event);
107	sixie	1.9	eventInfo_->SetFlavorHistoryPath(*flavorHistoryPath);
108			}
109	loizides	1.1	}