TreeFiller/src/FillerDecayParts.cc

// $Id: FillerDecayParts.cc,v 1.8 2008/09/27 05:50:47 loizides Exp $

#include "MitAna/DataTree/interface/DecayParticle.h"
#include "MitAna/DataTree/interface/DaughterData.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "MitEdm/DataFormats/interface/Collections.h"
#include "MitEdm/DataFormats/interface/DecayPart.h"
#include "MitEdm/DataFormats/interface/BasePart.h"
#include "MitEdm/DataFormats/interface/BasePartFwd.h"
#include "MitProd/TreeFiller/interface/FillerDecayParts.h"
#include "MitAna/DataTree/interface/Names.h"

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

//--------------------------------------------------------------------------------------------------
FillerDecayParts::FillerDecayParts(const ParameterSet &cfg, const char *name, bool active) :
  BaseFiller(cfg,name,active),
  edmName_(Conf().getUntrackedParameter<string>("edmName","")),
  mitName_(Conf().getUntrackedParameter<string>("mitName","")),
  vertexMapName_(Conf().getUntrackedParameter<string>("vertexMap","")),
  basePartMapNames_(Conf().exists("basePartMaps") ? 
                    Conf().getUntrackedParameter<vector<string> >("basePartMaps") : 
                    vector<string>()),
  vertexMap_(0),
  decays_(new mithep::DecayParticleArr(250)),
  stableData_(new mithep::StableDataArr(500)),
  decayData_(new mithep::DecayDataArr(500))
{
  // Constructor.
}

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

  delete decays_;
  delete stableData_;
  delete decayData_;
}

//--------------------------------------------------------------------------------------------------
void FillerDecayParts::BookDataBlock(TreeWriter &tws)
{
  // Add tracks branch to tree and get our map.
  std::string stableDataName = mitName_ + "_StableDatas";
  std::string decayDataName  = mitName_ + "_DecayDatas";
  
  tws.AddBranch(mitName_.c_str(),&decays_);
  tws.AddBranch(stableDataName.c_str(), &stableData_);
  tws.AddBranch(decayDataName.c_str(), &decayData_);

  if (!vertexMapName_.empty())
    vertexMap_ = OS()->get<VertexMap>(vertexMapName_.c_str());
  
  for (std::vector<std::string>::const_iterator bmapName = basePartMapNames_.begin();
        bmapName!=basePartMapNames_.end(); ++bmapName) {
    if (!bmapName->empty()) 
      basePartMaps_.push_back(OS()->get<BasePartMap>(bmapName->c_str()));
  }
}

//--------------------------------------------------------------------------------------------------
void FillerDecayParts::FillDataBlock(const edm::Event      &evt, 
                                     const edm::EventSetup &setup)
{
  // Fill our EDM DecayPart collection into the MIT DecayParticle collection.
  decays_->Reset();
  stableData_->Reset();
  decayData_->Reset();

  Handle<mitedm::DecayPartCol> hParts;
  GetProduct(edmName_, hParts, evt);  
  const mitedm::DecayPartCol *iParts = hParts.product();
  
  // loop through all decayParts and fill the information
  for (UInt_t i=0; i<iParts->size(); ++i) {
    const mitedm::DecayPart &p =iParts->at(i);                    // for convenience
    mithep::DecayParticle *d = decays_->Allocate();
    new (d) mithep::DecayParticle(p.pid());
    
    d->SetMassError(p.fittedMassError());
    d->SetLxy(p.lxyToPv());
    d->SetLxyError(p.lxyToPvError());
    d->SetDxy(p.dxyToPv());
    d->SetDxyError(p.dxyToPvError());
    d->SetLz(p.lzToPv());
    d->SetLzError(p.lzToPvError());
    d->SetMom(p.fourMomentum());
    d->SetChi2(p.chi2());
    d->SetNdof(p.ndof());

    if (p.primaryVertex().isNonnull())
      d->SetPriVertex(vertexMap_->GetMit(p.primaryVertex()));
    
    if (basePartMaps_.size()) {
       //loop through and add stable daughters
      for (Int_t j=0; j<p.nStableChild(); ++j) {
        const mitedm::StableData &stable = p.getStableData(j);
        mitedm::BasePartPtr thePtr = stable.originalPtr();
        mithep::Particle *daughter = getMitParticle(thePtr);
        
        mithep::StableData *outStable = stableData_->Allocate();
        new (outStable) mithep::StableData(daughter,stable.p3().x(),stable.p3().y(),stable.p3().z());
        
        d->AddDaughterData(outStable);
      }
      //loop through and add decay daughters
      for (Int_t j=0; j<p.nDecayChild(); ++j) {
        const mitedm::DecayData &decay = p.getDecayData(j);
        mitedm::BasePartPtr thePtr = decay.originalPtr();
        mithep::Particle *daughter = getMitParticle(thePtr);
        
        mithep::DecayData *outDecay = decayData_->Allocate();
        new (outDecay) mithep::DecayData(daughter,decay.p4());
        outDecay->SetMassError(decay.massErr());
        outDecay->SetLxy(decay.lxy());
        outDecay->SetLxyError(decay.lxyErr());
        outDecay->SetLz(decay.lz());
        outDecay->SetLzError(decay.lzErr());
        outDecay->SetDxy(decay.dxy());
        outDecay->SetDxyError(decay.dxyErr());
        
        d->AddDaughterData(outDecay);
      }
    }
/*    printf("mitedm::DecayPart              position=%5f,%5f,%5f\n",p.position().x(),
                                                          p.position().y(),
                                                          p.position().z());
    printf("mithep::DecayParticle          position=%5f,%5f,%5f\n",d->Position().X(),
                                                          d->Position().Y(),
                                                          d->Position().Z());
    printf("mithep::DecayParticle relative position=%5f,%5f,%5f\n",d->RelativePosition().X(),
                                                          d->RelativePosition().Y(),
                                                          d->RelativePosition().Z()); */                                                                                                    
    //cout << "MITHEP...\n";d->print();
  }
  decays_->Trim();
  stableData_->Trim();
  decayData_->Trim();
}

//--------------------------------------------------------------------------------------------------
mithep::Particle *FillerDecayParts::getMitParticle(mitedm::BasePartPtr ptr) const
{
  mithep::Particle *mitPart = 0;
  for (std::vector<const mithep::BasePartMap*>::const_iterator bmap = basePartMaps_.begin();
        bmap!=basePartMaps_.end(); ++bmap) {
    const mithep::BasePartMap *theMap = *bmap;
    if (theMap->HasMit(ptr))
      mitPart = theMap->GetMit(ptr);
  }
  
  if (!mitPart)
    throw edm::Exception(edm::errors::Configuration, "FillerDecayParts::FillDataBlock()\n")
    << "Error! MITHEP Object " 
    << "not found in AssociationMaps (" << typeid(*this).name() << ")." << std::endl;
    
  return mitPart;
}
Revision:	1.9
Committed:	Tue Sep 30 13:03:16 2008 UTC (16 years, 7 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_004
Changes since 1.8:	+88 -46 lines
Log Message:	updated for new DecayPart and DecayParticle formats and added primary vertex links
#	User	Rev	Content
1	bendavid	1.9	// $Id: FillerDecayParts.cc,v 1.8 2008/09/27 05:50:47 loizides Exp $
2	paus	1.1
3	loizides	1.4	#include "MitAna/DataTree/interface/DecayParticle.h"
4	bendavid	1.6	#include "MitAna/DataTree/interface/DaughterData.h"
5	paus	1.1	#include "FWCore/MessageLogger/interface/MessageLogger.h"
6			#include "DataFormats/Common/interface/Handle.h"
7			#include "DataFormats/TrackReco/interface/Track.h"
8			#include "DataFormats/TrackReco/interface/TrackFwd.h"
9			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
10			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
11			#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
12	loizides	1.8	#include "MitEdm/DataFormats/interface/Collections.h"
13	bendavid	1.2	#include "MitEdm/DataFormats/interface/DecayPart.h"
14			#include "MitEdm/DataFormats/interface/BasePart.h"
15			#include "MitEdm/DataFormats/interface/BasePartFwd.h"
16	paus	1.1	#include "MitProd/TreeFiller/interface/FillerDecayParts.h"
17			#include "MitAna/DataTree/interface/Names.h"
18
19			using namespace std;
20			using namespace edm;
21			using namespace mithep;
22
23			//--------------------------------------------------------------------------------------------------
24	loizides	1.4	FillerDecayParts::FillerDecayParts(const ParameterSet &cfg, const char *name, bool active) :
25	paus	1.1	BaseFiller(cfg,name,active),
26	loizides	1.4	edmName_(Conf().getUntrackedParameter<string>("edmName","")),
27			mitName_(Conf().getUntrackedParameter<string>("mitName","")),
28	bendavid	1.9	vertexMapName_(Conf().getUntrackedParameter<string>("vertexMap","")),
29	loizides	1.8	basePartMapNames_(Conf().exists("basePartMaps") ?
30			Conf().getUntrackedParameter<vector<string> >("basePartMaps") :
31			vector<string>()),
32	bendavid	1.9	vertexMap_(0),
33			decays_(new mithep::DecayParticleArr(250)),
34			stableData_(new mithep::StableDataArr(500)),
35			decayData_(new mithep::DecayDataArr(500))
36	paus	1.1	{
37	loizides	1.4	// Constructor.
38	paus	1.1	}
39
40			//--------------------------------------------------------------------------------------------------
41			FillerDecayParts::~FillerDecayParts()
42			{
43	loizides	1.4	// Destructor.
44
45	paus	1.1	delete decays_;
46	bendavid	1.9	delete stableData_;
47			delete decayData_;
48	paus	1.1	}
49
50			//--------------------------------------------------------------------------------------------------
51			void FillerDecayParts::BookDataBlock(TreeWriter &tws)
52			{
53	loizides	1.4	// Add tracks branch to tree and get our map.
54	bendavid	1.9	std::string stableDataName = mitName_ + "_StableDatas";
55			std::string decayDataName = mitName_ + "_DecayDatas";
56
57			tws.AddBranch(mitName_.c_str(),&decays_);
58			tws.AddBranch(stableDataName.c_str(), &stableData_);
59			tws.AddBranch(decayDataName.c_str(), &decayData_);
60	loizides	1.4
61	bendavid	1.9	if (!vertexMapName_.empty())
62			vertexMap_ = OS()->get<VertexMap>(vertexMapName_.c_str());
63
64	bendavid	1.6	for (std::vector<std::string>::const_iterator bmapName = basePartMapNames_.begin();
65			bmapName!=basePartMapNames_.end(); ++bmapName) {
66			if (!bmapName->empty())
67			basePartMaps_.push_back(OS()->get<BasePartMap>(bmapName->c_str()));
68			}
69	paus	1.1	}
70
71			//--------------------------------------------------------------------------------------------------
72			void FillerDecayParts::FillDataBlock(const edm::Event &evt,
73			const edm::EventSetup &setup)
74			{
75	loizides	1.4	// Fill our EDM DecayPart collection into the MIT DecayParticle collection.
76	paus	1.1	decays_->Reset();
77	bendavid	1.9	stableData_->Reset();
78			decayData_->Reset();
79	loizides	1.4
80	bendavid	1.2	Handle<mitedm::DecayPartCol> hParts;
81	paus	1.1	GetProduct(edmName_, hParts, evt);
82	bendavid	1.2	const mitedm::DecayPartCol *iParts = hParts.product();
83	paus	1.1
84			// loop through all decayParts and fill the information
85	bendavid	1.2	for (UInt_t i=0; i<iParts->size(); ++i) {
86			const mitedm::DecayPart &p =iParts->at(i); // for convenience
87			mithep::DecayParticle *d = decays_->Allocate();
88	bendavid	1.9	new (d) mithep::DecayParticle(p.pid());
89	bendavid	1.2
90	bendavid	1.9	d->SetMassError(p.fittedMassError());
91			d->SetLxy(p.lxyToPv());
92			d->SetLxyError(p.lxyToPvError());
93			d->SetDxy(p.dxyToPv());
94			d->SetDxyError(p.dxyToPvError());
95			d->SetLz(p.lzToPv());
96			d->SetLzError(p.lzToPvError());
97	bendavid	1.2	d->SetMom(p.fourMomentum());
98	mrudolph	1.5	d->SetChi2(p.chi2());
99			d->SetNdof(p.ndof());
100
101	bendavid	1.9	if (p.primaryVertex().isNonnull())
102			d->SetPriVertex(vertexMap_->GetMit(p.primaryVertex()));
103
104	bendavid	1.6	if (basePartMaps_.size()) {
105	bendavid	1.9	//loop through and add stable daughters
106			for (Int_t j=0; j<p.nStableChild(); ++j) {
107			const mitedm::StableData &stable = p.getStableData(j);
108			mitedm::BasePartPtr thePtr = stable.originalPtr();
109			mithep::Particle *daughter = getMitParticle(thePtr);
110
111			mithep::StableData *outStable = stableData_->Allocate();
112			new (outStable) mithep::StableData(daughter,stable.p3().x(),stable.p3().y(),stable.p3().z());
113	bendavid	1.6
114	bendavid	1.9	d->AddDaughterData(outStable);
115			}
116			//loop through and add decay daughters
117			for (Int_t j=0; j<p.nDecayChild(); ++j) {
118			const mitedm::DecayData &decay = p.getDecayData(j);
119			mitedm::BasePartPtr thePtr = decay.originalPtr();
120			mithep::Particle *daughter = getMitParticle(thePtr);
121
122			mithep::DecayData *outDecay = decayData_->Allocate();
123			new (outDecay) mithep::DecayData(daughter,decay.p4());
124			outDecay->SetMassError(decay.massErr());
125			outDecay->SetLxy(decay.lxy());
126			outDecay->SetLxyError(decay.lxyErr());
127			outDecay->SetLz(decay.lz());
128			outDecay->SetLzError(decay.lzErr());
129			outDecay->SetDxy(decay.dxy());
130			outDecay->SetDxyError(decay.dxyErr());
131
132			d->AddDaughterData(outDecay);
133	bendavid	1.2	}
134			}
135	bendavid	1.9	/* printf("mitedm::DecayPart position=%5f,%5f,%5f\n",p.position().x(),
136			p.position().y(),
137			p.position().z());
138			printf("mithep::DecayParticle position=%5f,%5f,%5f\n",d->Position().X(),
139			d->Position().Y(),
140			d->Position().Z());
141			printf("mithep::DecayParticle relative position=%5f,%5f,%5f\n",d->RelativePosition().X(),
142			d->RelativePosition().Y(),
143			d->RelativePosition().Z()); */
144			//cout << "MITHEP...\n";d->print();
145	paus	1.1	}
146			decays_->Trim();
147	bendavid	1.9	stableData_->Trim();
148			decayData_->Trim();
149			}
150
151			//--------------------------------------------------------------------------------------------------
152			mithep::Particle *FillerDecayParts::getMitParticle(mitedm::BasePartPtr ptr) const
153			{
154			mithep::Particle *mitPart = 0;
155			for (std::vector<const mithep::BasePartMap*>::const_iterator bmap = basePartMaps_.begin();
156			bmap!=basePartMaps_.end(); ++bmap) {
157			const mithep::BasePartMap theMap = bmap;
158			if (theMap->HasMit(ptr))
159			mitPart = theMap->GetMit(ptr);
160			}
161
162			if (!mitPart)
163			throw edm::Exception(edm::errors::Configuration, "FillerDecayParts::FillDataBlock()\n")
164			<< "Error! MITHEP Object "
165			<< "not found in AssociationMaps (" << typeid(*this).name() << ")." << std::endl;
166
167			return mitPart;
168	paus	1.1	}