TreeFiller/src/FillerCaloTowers.cc

// $Id: FillerCaloTowers.cc,v 1.17 2010/06/25 15:18:38 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerCaloTowers.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "MitAna/DataTree/interface/CaloTowerCol.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitProd/ObjectService/interface/ObjectService.h"

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

//--------------------------------------------------------------------------------------------------
FillerCaloTowers::FillerCaloTowers(const ParameterSet &cfg, const char *name, bool active) : 
  BaseFiller(cfg,name,active),
  edmName_(Conf().getUntrackedParameter<string>("edmName","towerMaker")),
  mitName_(Conf().getUntrackedParameter<string>("mitName","CaloTowers")),
  caloTowerMapName_(Conf().getUntrackedParameter<string>("caloTowerMapName", "CaloTowerMap")),
  caloTowerDetIdMapName_(Conf().getUntrackedParameter<string>("caloTowerDetIdMapName", "CaloTowerDetIdMap")),
  caloTowers_(new mithep::CaloTowerArr(1000)),
  caloTowerMap_(new mithep::CaloTowerMap),
  caloTowerDetIdMap_(new mithep::CaloTowerDetIdMap)
{
  // Constructor.
}

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

  delete caloTowers_;
  delete caloTowerMap_;
  delete caloTowerDetIdMap_;
}

//--------------------------------------------------------------------------------------------------
void FillerCaloTowers::BookDataBlock(TreeWriter &tws)
{
  // Add CaloTower branch to tree.

  Int_t brsize = tws.GetDefaultBrSize();
  if (brsize<32*1024)
    brsize=32*1024;

  tws.AddBranch(mitName_,&caloTowers_,brsize);
  OS()->add<CaloTowerArr>(caloTowers_,mitName_);

  if (!caloTowerMapName_.empty()) {
    caloTowerMap_->SetBrName(mitName_);
    OS()->add<CaloTowerMap>(caloTowerMap_,caloTowerMapName_);
  }
  
  if (!caloTowerDetIdMapName_.empty()) {
    caloTowerDetIdMap_->SetBrName(mitName_);
    OS()->add<CaloTowerDetIdMap>(caloTowerDetIdMap_,caloTowerDetIdMapName_);
  }
}

//--------------------------------------------------------------------------------------------------
void FillerCaloTowers::FillDataBlock(const edm::Event      &event, 
                                     const edm::EventSetup &setup)
{
  // Fill the CaloTower info.

  caloTowers_->Delete();
  caloTowerMap_->Reset();
  caloTowerDetIdMap_->Reset();

  Handle<CaloTowerCollection> hCaloTowerProduct;
  GetProduct(edmName_, hCaloTowerProduct, event);
  const CaloTowerCollection inCaloTowers = *(hCaloTowerProduct.product());  

  for (CaloTowerCollection::const_iterator inCaloTower = inCaloTowers.begin(); 
       inCaloTower != inCaloTowers.end(); ++inCaloTower) {
    
    mithep::CaloTower *outCaloTower = caloTowers_->Allocate();
    new (outCaloTower) mithep::CaloTower();
       
    double deltaE = inCaloTower->energy() - inCaloTower->emEnergy() - inCaloTower->hadEnergy();
       
    if ( deltaE > 1e-3)
      throw edm::Exception(edm::errors::Configuration, "FillerCaloTowers::FillDataBlock()\n")
        << "Error! reco::CaloTower default energy does not exclude HO " 
        << "as assumed by mithep::CaloTower, deltaE = " << deltaE << std::endl;

    ThreeVector avgPos;
    ThreeVector emPos(inCaloTower->emPosition().x(),
                      inCaloTower->emPosition().y(),
                      inCaloTower->emPosition().z());
    ThreeVector hadPos(inCaloTower->hadPosition().x(),
                       inCaloTower->hadPosition().y(),
                       inCaloTower->hadPosition().z());
    double energyWeight = 0.0;
    if (inCaloTower->emEnergy()>0.0) {
      avgPos += inCaloTower->emEnergy()*emPos;
      energyWeight += inCaloTower->emEnergy();
    }
    if (inCaloTower->hadEnergy()>0.0) {
      avgPos += inCaloTower->hadEnergy()*hadPos;
      energyWeight += inCaloTower->hadEnergy();
    }
    if (energyWeight > 0.0)
      avgPos = (1.0/energyWeight)*avgPos;
    else
      avgPos = emPos;
      
    ThreeVectorC finalPos(avgPos.Rho(),inCaloTower->eta(),inCaloTower->phi());

    outCaloTower->SetPosition(finalPos.x(),finalPos.y(),finalPos.z());

    outCaloTower->SetEmEnergy(inCaloTower->emEnergy());
    outCaloTower->SetHadEnergy(inCaloTower->hadEnergy());
    outCaloTower->SetOuterEnergy(inCaloTower->outerEnergy());
    outCaloTower->SetMass(inCaloTower->mass());
    
    //fill detid associations (used for supercluster-calotower assocations)
    for (std::vector<DetId>::const_iterator itId = inCaloTower->constituents().begin(); itId!=inCaloTower->constituents().end(); ++itId) {
      caloTowerDetIdMap_->Add(*itId,outCaloTower);
    }
    
    //fill edm-mit association
    caloTowerMap_->Add(inCaloTower->id(),outCaloTower);
    
    if (verbose_>1) {
      printf("reco::CaloTower:   pt=%5f, eta=%5f, phi=%5f, mass=%5f\n",
             inCaloTower->pt(),inCaloTower->eta(),inCaloTower->phi(),inCaloTower->mass());
      printf("mithep::CaloTower: pt=%5f, eta=%5f, phi=%5f, mass=%5f\n",
             outCaloTower->Pt(),outCaloTower->Eta(),outCaloTower->Phi(),outCaloTower->Mass());
      printf("emPos : rho=%5f, eta=%5f, phi=%5f\n",emPos.Rho(),emPos.Eta(),emPos.Phi());
      printf("hadPos: rho=%5f, eta=%5f, phi=%5f\n",hadPos.Rho(),hadPos.Eta(),hadPos.Phi());
      printf("avgPos: rho=%5f, eta=%5f, phi=%5f\n",avgPos.Rho(),avgPos.Eta(),avgPos.Phi());
      printf("finPos: rho=%5f, eta=%5f, phi=%5f\n",finalPos.Rho(),finalPos.Eta(),finalPos.Phi());
    }
  }
  caloTowers_->Trim();
}
Revision:	1.18
Committed:	Wed Aug 18 01:36:11 2010 UTC (14 years, 8 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d
Changes since 1.17:	+2 -1 lines
Log Message:	fixes for broken calotower-supercluster associations
#	User	Rev	Content
1	bendavid	1.18	// $Id: FillerCaloTowers.cc,v 1.17 2010/06/25 15:18:38 bendavid Exp $
2	sixie	1.1
3			#include "MitProd/TreeFiller/interface/FillerCaloTowers.h"
4			#include "DataFormats/CaloTowers/interface/CaloTower.h"
5			#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
6	loizides	1.11	#include "FWCore/Framework/interface/ESHandle.h"
7	sixie	1.1	#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
8			#include "Geometry/Records/interface/IdealGeometryRecord.h"
9	loizides	1.11	#include "MitAna/DataTree/interface/CaloTowerCol.h"
10			#include "MitAna/DataTree/interface/Names.h"
11			#include "MitProd/ObjectService/interface/ObjectService.h"
12	sixie	1.1
13			using namespace std;
14			using namespace edm;
15			using namespace mithep;
16
17			//--------------------------------------------------------------------------------------------------
18			FillerCaloTowers::FillerCaloTowers(const ParameterSet &cfg, const char *name, bool active) :
19			BaseFiller(cfg,name,active),
20			edmName_(Conf().getUntrackedParameter<string>("edmName","towerMaker")),
21			mitName_(Conf().getUntrackedParameter<string>("mitName","CaloTowers")),
22	loizides	1.9	caloTowerMapName_(Conf().getUntrackedParameter<string>("caloTowerMapName", "CaloTowerMap")),
23	bendavid	1.17	caloTowerDetIdMapName_(Conf().getUntrackedParameter<string>("caloTowerDetIdMapName", "CaloTowerDetIdMap")),
24	bendavid	1.3	caloTowers_(new mithep::CaloTowerArr(1000)),
25	bendavid	1.17	caloTowerMap_(new mithep::CaloTowerMap),
26			caloTowerDetIdMap_(new mithep::CaloTowerDetIdMap)
27	sixie	1.1	{
28			// Constructor.
29			}
30
31			//--------------------------------------------------------------------------------------------------
32			FillerCaloTowers::~FillerCaloTowers()
33			{
34			// Destructor.
35
36			delete caloTowers_;
37	bendavid	1.3	delete caloTowerMap_;
38	bendavid	1.17	delete caloTowerDetIdMap_;
39	sixie	1.1	}
40
41			//--------------------------------------------------------------------------------------------------
42	bendavid	1.16	void FillerCaloTowers::BookDataBlock(TreeWriter &tws)
43	sixie	1.1	{
44			// Add CaloTower branch to tree.
45
46	loizides	1.10	Int_t brsize = tws.GetDefaultBrSize();
47			if (brsize<32*1024)
48			brsize=32*1024;
49
50			tws.AddBranch(mitName_,&caloTowers_,brsize);
51	loizides	1.9	OS()->add<CaloTowerArr>(caloTowers_,mitName_);
52	sixie	1.1
53	loizides	1.9	if (!caloTowerMapName_.empty()) {
54			caloTowerMap_->SetBrName(mitName_);
55			OS()->add<CaloTowerMap>(caloTowerMap_,caloTowerMapName_);
56			}
57	bendavid	1.17
58			if (!caloTowerDetIdMapName_.empty()) {
59			caloTowerDetIdMap_->SetBrName(mitName_);
60			OS()->add<CaloTowerDetIdMap>(caloTowerDetIdMap_,caloTowerDetIdMapName_);
61			}
62	sixie	1.1	}
63
64			//--------------------------------------------------------------------------------------------------
65			void FillerCaloTowers::FillDataBlock(const edm::Event &event,
66	loizides	1.9	const edm::EventSetup &setup)
67	sixie	1.1	{
68	loizides	1.4	// Fill the CaloTower info.
69	sixie	1.1
70	bendavid	1.8	caloTowers_->Delete();
71	bendavid	1.3	caloTowerMap_->Reset();
72	bendavid	1.18	caloTowerDetIdMap_->Reset();
73	sixie	1.1
74			Handle<CaloTowerCollection> hCaloTowerProduct;
75			GetProduct(edmName_, hCaloTowerProduct, event);
76			const CaloTowerCollection inCaloTowers = *(hCaloTowerProduct.product());
77
78			for (CaloTowerCollection::const_iterator inCaloTower = inCaloTowers.begin();
79			inCaloTower != inCaloTowers.end(); ++inCaloTower) {
80
81			mithep::CaloTower *outCaloTower = caloTowers_->Allocate();
82			new (outCaloTower) mithep::CaloTower();
83
84	bendavid	1.6	double deltaE = inCaloTower->energy() - inCaloTower->emEnergy() - inCaloTower->hadEnergy();
85	bendavid	1.13
86	bendavid	1.7	if ( deltaE > 1e-3)
87	bendavid	1.6	throw edm::Exception(edm::errors::Configuration, "FillerCaloTowers::FillDataBlock()\n")
88	loizides	1.9	<< "Error! reco::CaloTower default energy does not exclude HO "
89			<< "as assumed by mithep::CaloTower, deltaE = " << deltaE << std::endl;
90	bendavid	1.12
91			ThreeVector avgPos;
92			ThreeVector emPos(inCaloTower->emPosition().x(),
93			inCaloTower->emPosition().y(),
94			inCaloTower->emPosition().z());
95			ThreeVector hadPos(inCaloTower->hadPosition().x(),
96			inCaloTower->hadPosition().y(),
97			inCaloTower->hadPosition().z());
98			double energyWeight = 0.0;
99			if (inCaloTower->emEnergy()>0.0) {
100			avgPos += inCaloTower->emEnergy()*emPos;
101			energyWeight += inCaloTower->emEnergy();
102			}
103			if (inCaloTower->hadEnergy()>0.0) {
104			avgPos += inCaloTower->hadEnergy()*hadPos;
105			energyWeight += inCaloTower->hadEnergy();
106			}
107			if (energyWeight > 0.0)
108			avgPos = (1.0/energyWeight)*avgPos;
109			else
110			avgPos = emPos;
111
112			ThreeVectorC finalPos(avgPos.Rho(),inCaloTower->eta(),inCaloTower->phi());
113
114			outCaloTower->SetPosition(finalPos.x(),finalPos.y(),finalPos.z());
115
116	sixie	1.1	outCaloTower->SetEmEnergy(inCaloTower->emEnergy());
117			outCaloTower->SetHadEnergy(inCaloTower->hadEnergy());
118			outCaloTower->SetOuterEnergy(inCaloTower->outerEnergy());
119	bendavid	1.12	outCaloTower->SetMass(inCaloTower->mass());
120	bendavid	1.17
121			//fill detid associations (used for supercluster-calotower assocations)
122			for (std::vector<DetId>::const_iterator itId = inCaloTower->constituents().begin(); itId!=inCaloTower->constituents().end(); ++itId) {
123			caloTowerDetIdMap_->Add(*itId,outCaloTower);
124			}
125
126			//fill edm-mit association
127	bendavid	1.3	caloTowerMap_->Add(inCaloTower->id(),outCaloTower);
128	bendavid	1.12
129			if (verbose_>1) {
130	loizides	1.14	printf("reco::CaloTower: pt=%5f, eta=%5f, phi=%5f, mass=%5f\n",
131			inCaloTower->pt(),inCaloTower->eta(),inCaloTower->phi(),inCaloTower->mass());
132			printf("mithep::CaloTower: pt=%5f, eta=%5f, phi=%5f, mass=%5f\n",
133			outCaloTower->Pt(),outCaloTower->Eta(),outCaloTower->Phi(),outCaloTower->Mass());
134	bendavid	1.12	printf("emPos : rho=%5f, eta=%5f, phi=%5f\n",emPos.Rho(),emPos.Eta(),emPos.Phi());
135			printf("hadPos: rho=%5f, eta=%5f, phi=%5f\n",hadPos.Rho(),hadPos.Eta(),hadPos.Phi());
136			printf("avgPos: rho=%5f, eta=%5f, phi=%5f\n",avgPos.Rho(),avgPos.Eta(),avgPos.Phi());
137			printf("finPos: rho=%5f, eta=%5f, phi=%5f\n",finalPos.Rho(),finalPos.Eta(),finalPos.Phi());
138			}
139	sixie	1.1	}
140			caloTowers_->Trim();
141			}