ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerSuperClusters.cc
Revision: 1.17
Committed: Fri Dec 28 17:27:21 2012 UTC (12 years, 4 months ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029b, Mit_030_pre1, Mit_029a
Changes since 1.16: +9 -7 lines
Log Message:
Added Embedded and PFAOD functionality

File Contents

# User Rev Content
1 pharris 1.17 // $Id: FillerSuperClusters.cc,v 1.16 2012/05/05 16:49:59 paus Exp $
2 sixie 1.1
3     #include "MitProd/TreeFiller/interface/FillerSuperClusters.h"
4     #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
5     #include "DataFormats/EgammaReco/interface/SuperCluster.h"
6     #include "MitAna/DataTree/interface/BasicCluster.h"
7 paus 1.16 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
8     #include "DataFormats/EgammaReco/interface/PreshowerClusterFwd.h"
9 loizides 1.5 #include "MitAna/DataTree/interface/SuperClusterCol.h"
10 sixie 1.1 #include "MitAna/DataTree/interface/Names.h"
11 loizides 1.5 #include "MitProd/ObjectService/interface/ObjectService.h"
12 sixie 1.10 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
13 paus 1.16 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
14     #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
15     #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
16     #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
17     #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
18     #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
19     #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
20     #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
21     #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
22     #include "Geometry/Records/interface/CaloGeometryRecord.h"
23 sixie 1.1
24     using namespace std;
25     using namespace edm;
26     using namespace mithep;
27    
28     //--------------------------------------------------------------------------------------------------
29     FillerSuperClusters::FillerSuperClusters(const ParameterSet &cfg, const char *name, bool active) :
30 paus 1.16 BaseFiller (cfg,name,active),
31     edmName_ (Conf().getUntrackedParameter<string>("edmName","hybridSuperClusters")),
32     mitName_ (Conf().getUntrackedParameter<string>("mitName","SuperClusters")),
33     basicClusterMapName_ (Conf().getUntrackedParameter<string>("basicClusterMapName",
34     "BasicClusterMap")),
35     psClusterMapName_ (Conf().getUntrackedParameter<string>("psClusterMapName",
36     "")),
37 bendavid 1.11 caloTowerDetIdMapName_(Conf().getUntrackedParameter<string>("caloTowerDetIdMapName",
38 paus 1.16 "CaloTowerDetIdMap")),
39     superClusterMapName_ (Conf().getUntrackedParameter<string>("superClusterMapName",
40     "SuperClusterMap")),
41 bendavid 1.2 superClusterIdMapName_(Conf().getUntrackedParameter<string>("superClusterIdMapName",
42 loizides 1.7 "SuperClusterIdMap")),
43 paus 1.16 caloTowerName_ (Conf().getUntrackedParameter<string>("caloTowerName","towerMaker")),
44     basicClusterMap_ (0),
45     psClusterMap_ (0),
46     caloTowerDetIdMap_ (0),
47     superClusters_ (new mithep::SuperClusterArr(25)),
48     superClusterMap_ (new mithep::SuperClusterMap),
49     superClusterIdMap_ (new mithep::SuperClusterIdMap)
50 sixie 1.1 {
51     // Constructor.
52     }
53    
54     //--------------------------------------------------------------------------------------------------
55     FillerSuperClusters::~FillerSuperClusters()
56     {
57     // Destructor.
58    
59     delete superClusters_;
60     delete superClusterMap_;
61 bendavid 1.2 delete superClusterIdMap_;
62 sixie 1.1 }
63    
64     //--------------------------------------------------------------------------------------------------
65 bendavid 1.9 void FillerSuperClusters::BookDataBlock(TreeWriter &tws)
66 sixie 1.1 {
67 loizides 1.4 // Add super cluster branch to tree and get pointers to maps.
68    
69     tws.AddBranch(mitName_,&superClusters_);
70     OS()->add<SuperClusterArr>(superClusters_,mitName_);
71    
72     if (!basicClusterMapName_.empty()) {
73     basicClusterMap_ = OS()->get<BasicClusterMap>(basicClusterMapName_);
74     if (basicClusterMap_)
75     AddBranchDep(mitName_,basicClusterMap_->GetBrName());
76     }
77 paus 1.16
78     if (!psClusterMapName_.empty()) {
79     psClusterMap_ = OS()->get<PsClusterMap>(psClusterMapName_);
80     if (psClusterMap_)
81     AddBranchDep(mitName_,psClusterMap_->GetBrName());
82     }
83    
84 bendavid 1.11 if (!caloTowerDetIdMapName_.empty()) {
85     caloTowerDetIdMap_ = OS()->get<CaloTowerDetIdMap>(caloTowerDetIdMapName_);
86     if (caloTowerDetIdMap_)
87     AddBranchDep(mitName_,caloTowerDetIdMap_->GetBrName());
88     }
89    
90 loizides 1.4 if (!superClusterMapName_.empty()) {
91     superClusterMap_->SetBrName(mitName_);
92     OS()->add<SuperClusterMap>(superClusterMap_,superClusterMapName_);
93     }
94     if (!superClusterIdMapName_.empty()) {
95     superClusterIdMap_->SetBrName(mitName_);
96     OS()->add<SuperClusterIdMap>(superClusterIdMap_,superClusterIdMapName_);
97     }
98 sixie 1.1 }
99    
100     //--------------------------------------------------------------------------------------------------
101     void FillerSuperClusters::FillDataBlock(const edm::Event &event,
102 loizides 1.4 const edm::EventSetup &setup)
103 sixie 1.1 {
104 loizides 1.4 // Fill the collection.
105    
106 bendavid 1.3 superClusters_->Delete();
107 sixie 1.1 superClusterMap_->Reset();
108 bendavid 1.2 superClusterIdMap_->Reset();
109 sixie 1.1
110     Handle<reco::SuperClusterCollection> hSuperClusterProduct;
111     GetProduct(edmName_, hSuperClusterProduct, event);
112 pharris 1.17
113 sixie 1.10 Handle<CaloTowerCollection> hCaloTowerProduct;
114 pharris 1.17 if(caloTowerName_.size() > 0) GetProduct(caloTowerName_, hCaloTowerProduct, event);
115 sixie 1.10
116 sixie 1.1 superClusterMap_->SetEdmProductId(hSuperClusterProduct.id().id());
117     const reco::SuperClusterCollection inSuperClusters = *(hSuperClusterProduct.product());
118    
119 paus 1.16 //lazy tools for supercluster ecal timing
120     EcalClusterLazyTools lazyTools(event, setup, edm::InputTag("reducedEcalRecHitsEB"),
121     edm::InputTag("reducedEcalRecHitsEE"));
122    
123     //mustache id
124     reco::Mustache mustache;
125    
126     //es shape variables
127     edm::Handle<EcalRecHitCollection> hESRecHits;
128     event.getByLabel("reducedEcalRecHitsES" , hESRecHits);
129    
130     edm::ESHandle<CaloGeometry> pGeometry;
131     setup.get<CaloGeometryRecord>().get(pGeometry);
132    
133     edm::ESHandle<CaloTopology> pTopology;
134     setup.get<CaloTopologyRecord>().get(pTopology);
135    
136     const CaloSubdetectorGeometry *geometryES = pGeometry->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
137     EcalPreshowerTopology topology_p(pGeometry);
138    
139     //map of preshower rechits for shape calculations
140     std::map<DetId, EcalRecHit> esmap;
141     if (hESRecHits.isValid()) {
142     EcalRecHitCollection::const_iterator it;
143     for (it = hESRecHits->begin(); it != hESRecHits->end(); ++it) {
144     // remove bad ES rechits
145     if (it->recoFlag()==1 || it->recoFlag()==14 || (it->recoFlag()<=10 && it->recoFlag()>=5)) continue;
146     //Make the map of DetID, EcalRecHit pairs
147     esmap.insert(std::make_pair(it->id(), *it));
148     }
149     }
150    
151 sixie 1.1 // loop through all super clusters
152     for (reco::SuperClusterCollection::const_iterator inSC = inSuperClusters.begin();
153     inSC != inSuperClusters.end(); ++inSC) {
154    
155     mithep::SuperCluster *outSC = superClusters_->Allocate();
156     new (outSC) mithep::SuperCluster();
157    
158     outSC->SetXYZ(inSC->x(),inSC->y(),inSC->z());
159     outSC->SetEnergy(inSC->energy());
160     outSC->SetRawEnergy(inSC->rawEnergy());
161     outSC->SetPreshowerEnergy(inSC->preshowerEnergy());
162     outSC->SetPhiWidth(inSC->phiWidth());
163     outSC->SetEtaWidth(inSC->etaWidth());
164    
165 sixie 1.10 //Compute Hadronic Energy behind the supercluster (within DR < 0.15)
166 pharris 1.17 if(caloTowerName_.size() > 0) {
167     EgammaTowerIsolation towerIsoDepth1(0.15,0.,0.,1,hCaloTowerProduct.product()) ;
168     EgammaTowerIsolation towerIsoDepth2(0.15,0.,0.,2,hCaloTowerProduct.product()) ;
169     outSC->SetHcalDepth1Energy(towerIsoDepth1.getTowerESum(&(*inSC)));
170     outSC->SetHcalDepth2Energy(towerIsoDepth2.getTowerESum(&(*inSC)));
171     }
172 sixie 1.10
173 paus 1.16 //ecal timing information
174     outSC->SetTime(lazyTools.SuperClusterTime(*inSC, event));
175     outSC->SetSeedTime(lazyTools.SuperClusterSeedTime(*inSC));
176    
177     //preshower shape variables
178     if (inSC->seed()->hitsAndFractions().at(0).first.subdetId()==EcalEndcap) {
179     std::vector<float> phoESShape = getESShape(getESHits(inSC->x(), inSC->y(), inSC->z(), esmap, *pGeometry.product(), &topology_p, 0));
180    
181     if (phoESShape.size()>=2) {
182     outSC->SetPsEffWidthSigmaXX(phoESShape[0]);
183     outSC->SetPsEffWidthSigmaYY(phoESShape[1]);
184     }
185    
186     }
187    
188    
189 loizides 1.4 // set the seed
190 sixie 1.1 if (basicClusterMap_ && inSC->seed().isNonnull())
191     outSC->SetSeed(basicClusterMap_->GetMit(inSC->seed()));
192    
193 paus 1.16 // add basic clusters that belong to this super cluster and tag with mustache id
194 bendavid 1.13 if (basicClusterMap_) {
195 paus 1.16 std::vector<reco::CaloCluster> bcs;
196     std::vector<unsigned int> insidebcs;
197     std::vector<unsigned int> outsidebcs;
198 bendavid 1.13 for(reco::CaloCluster_iterator bc = inSC->clustersBegin(); bc != inSC->clustersEnd(); ++bc) {
199 paus 1.16 if (bc->isNonnull()) {
200 bendavid 1.13 outSC->AddCluster(basicClusterMap_->GetMit(*bc));
201 paus 1.16 bcs.push_back(**bc);
202     }
203     }
204    
205     mustache.MustacheClust(bcs,insidebcs,outsidebcs);
206     for (unsigned int iclus=0; iclus<insidebcs.size(); ++iclus) {
207     const_cast<mithep::BasicCluster*>(outSC->Cluster(insidebcs.at(iclus)))->SetInsideMustache(kTRUE);
208     }
209     for (unsigned int iclus=0; iclus<outsidebcs.size(); ++iclus) {
210     const_cast<mithep::BasicCluster*>(outSC->Cluster(outsidebcs.at(iclus)))->SetInsideMustache(kFALSE);
211     }
212     }
213    
214     // add preshower clusters that belong to this super cluster and tag with mustache id
215     if (psClusterMap_) {
216     std::vector<reco::CaloCluster> bcs;
217     std::vector<unsigned int> insidebcs;
218     std::vector<unsigned int> outsidebcs;
219     for(reco::CaloCluster_iterator bc = inSC->preshowerClustersBegin(); bc != inSC->preshowerClustersEnd(); ++bc) {
220     if (bc->isNonnull()) {
221     outSC->AddPsClust(psClusterMap_->GetMit(*bc));
222     bcs.push_back(**bc);
223     }
224     }
225    
226     mustache.MustacheClust(bcs,insidebcs,outsidebcs);
227     for (unsigned int iclus=0; iclus<insidebcs.size(); ++iclus) {
228     const_cast<mithep::PsCluster*>(outSC->PsClust(insidebcs.at(iclus)))->SetInsideMustache(kTRUE);
229     }
230     for (unsigned int iclus=0; iclus<outsidebcs.size(); ++iclus) {
231     const_cast<mithep::PsCluster*>(outSC->PsClust(outsidebcs.at(iclus)))->SetInsideMustache(kFALSE);
232 bendavid 1.13 }
233 sixie 1.1 }
234    
235 paus 1.16 //compute preshower energy sums by plane
236     double psenergyplane1 = 0.;
237     double psenergyplane2 = 0.;
238     for(reco::CaloCluster_iterator bc = inSC->preshowerClustersBegin(); bc != inSC->preshowerClustersEnd(); ++bc) {
239     const reco::PreshowerCluster *pscluster = dynamic_cast<const reco::PreshowerCluster*>(bc->get());
240     if (!pscluster) continue;
241    
242     if (pscluster->plane()==1) psenergyplane1 += pscluster->energy();
243     else if (pscluster->plane()==2) psenergyplane2 += pscluster->energy();
244     }
245     outSC->SetPreshowerEnergyPlane1(psenergyplane1);
246     outSC->SetPreshowerEnergyPlane2(psenergyplane2);
247    
248 bendavid 1.11 //add super cluster det ids to the id map and also fill supercluster-calotower associations
249 bendavid 1.6 const std::vector< std::pair<DetId, float> > &pairs = inSC->hitsAndFractions();
250     for (std::vector< std::pair<DetId, float> >::const_iterator ipair = pairs.begin();
251     ipair < pairs.end(); ++ipair) {
252 bendavid 1.2
253 bendavid 1.6 const DetId &ihit = ipair->first;
254 bendavid 1.11
255     if (caloTowerDetIdMap_) {
256     if (caloTowerDetIdMap_->HasMit(ihit)) {
257 bendavid 1.12 const mithep::CaloTower *matchedTower = caloTowerDetIdMap_->GetMit(ihit);
258     if (!outSC->HasTower(matchedTower)) {
259     outSC->AddTower(matchedTower);
260     }
261 bendavid 1.11 }
262     }
263    
264 bendavid 1.6 superClusterIdMap_->Add(ihit,outSC);
265 bendavid 1.2 }
266 bendavid 1.11
267     // add super cluster to the map
268     reco::SuperClusterRef theRef(hSuperClusterProduct, inSC-inSuperClusters.begin());
269     superClusterMap_->Add(theRef, outSC);
270    
271 sixie 1.1 }
272     superClusters_->Trim();
273     }
274 paus 1.16
275     //horrible code below copied from globe for preshower cluster shape calculations
276     std::vector<float> FillerSuperClusters::getESHits(double X, double Y, double Z, std::map<DetId, EcalRecHit> rechits_map, const CaloGeometry& geometry, CaloSubdetectorTopology *topology_p, int row) {
277     std::vector<float> esHits;
278    
279     const GlobalPoint point(X,Y,Z);
280    
281     const CaloSubdetectorGeometry *geometry_p ;
282     geometry_p = geometry.getSubdetectorGeometry (DetId::Ecal,EcalPreshower) ;
283    
284     DetId esId1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
285     DetId esId2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
286     ESDetId esDetId1 = (esId1 == DetId(0)) ? ESDetId(0) : ESDetId(esId1);
287     ESDetId esDetId2 = (esId2 == DetId(0)) ? ESDetId(0) : ESDetId(esId2);
288    
289     std::map<DetId, EcalRecHit>::iterator it;
290     ESDetId next;
291     ESDetId strip1;
292     ESDetId strip2;
293    
294     strip1 = esDetId1;
295     strip2 = esDetId2;
296    
297     EcalPreshowerNavigator theESNav1(strip1, topology_p);
298     theESNav1.setHome(strip1);
299    
300     EcalPreshowerNavigator theESNav2(strip2, topology_p);
301     theESNav2.setHome(strip2);
302    
303     if (row == 1) {
304     if (strip1 != ESDetId(0)) strip1 = theESNav1.north();
305     if (strip2 != ESDetId(0)) strip2 = theESNav2.east();
306     } else if (row == -1) {
307     if (strip1 != ESDetId(0)) strip1 = theESNav1.south();
308     if (strip2 != ESDetId(0)) strip2 = theESNav2.west();
309     }
310    
311     // Plane 1
312     if (strip1 == ESDetId(0)) {
313     for (int i=0; i<31; ++i) esHits.push_back(0);
314     } else {
315    
316     it = rechits_map.find(strip1);
317     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
318     else esHits.push_back(0);
319     //cout<<"center : "<<strip1<<" "<<it->second.energy()<<endl;
320    
321     // east road
322     for (int i=0; i<15; ++i) {
323     next = theESNav1.east();
324     if (next != ESDetId(0)) {
325     it = rechits_map.find(next);
326     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
327     else esHits.push_back(0);
328     //cout<<"east "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
329     } else {
330     for (int j=i; j<15; j++) esHits.push_back(0);
331     break;
332     //cout<<"east "<<i<<" : "<<next<<" "<<0<<endl;
333     }
334     }
335    
336     // west road
337     theESNav1.setHome(strip1);
338     theESNav1.home();
339     for (int i=0; i<15; ++i) {
340     next = theESNav1.west();
341     if (next != ESDetId(0)) {
342     it = rechits_map.find(next);
343     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
344     else esHits.push_back(0);
345     //cout<<"west "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
346     } else {
347     for (int j=i; j<15; j++) esHits.push_back(0);
348     break;
349     //cout<<"west "<<i<<" : "<<next<<" "<<0<<endl;
350     }
351     }
352     }
353    
354     if (strip2 == ESDetId(0)) {
355     for (int i=0; i<31; ++i) esHits.push_back(0);
356     } else {
357    
358     it = rechits_map.find(strip2);
359     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
360     else esHits.push_back(0);
361     //cout<<"center : "<<strip2<<" "<<it->second.energy()<<endl;
362    
363     // north road
364     for (int i=0; i<15; ++i) {
365     next = theESNav2.north();
366     if (next != ESDetId(0)) {
367     it = rechits_map.find(next);
368     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
369     else esHits.push_back(0);
370     //cout<<"north "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
371     } else {
372     for (int j=i; j<15; j++) esHits.push_back(0);
373     break;
374     //cout<<"north "<<i<<" : "<<next<<" "<<0<<endl;
375     }
376     }
377    
378     // south road
379     theESNav2.setHome(strip2);
380     theESNav2.home();
381     for (int i=0; i<15; ++i) {
382     next = theESNav2.south();
383     if (next != ESDetId(0)) {
384     it = rechits_map.find(next);
385     if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
386     else esHits.push_back(0);
387     //cout<<"south "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
388     } else {
389     for (int j=i; j<15; j++) esHits.push_back(0);
390     break;
391     //cout<<"south "<<i<<" : "<<next<<" "<<0<<endl;
392     }
393     }
394     }
395    
396     return esHits;
397     }
398    
399     std::vector<float> FillerSuperClusters::getESShape(std::vector<float> ESHits0)
400     {
401     std::vector<float> esShape;
402    
403     const int nBIN = 21;
404     float esRH_F[nBIN];
405     float esRH_R[nBIN];
406     for (int idx=0; idx<nBIN; idx++) {
407     esRH_F[idx] = 0.;
408     esRH_R[idx] = 0.;
409     }
410    
411     for(int ibin=0; ibin<((nBIN+1)/2); ibin++) {
412     if (ibin==0) {
413     esRH_F[(nBIN-1)/2] = ESHits0[ibin];
414     esRH_R[(nBIN-1)/2] = ESHits0[ibin+31];
415     } else {
416     esRH_F[(nBIN-1)/2+ibin] = ESHits0[ibin];
417     esRH_F[(nBIN-1)/2-ibin] = ESHits0[ibin+15];
418     esRH_R[(nBIN-1)/2+ibin] = ESHits0[ibin+31];
419     esRH_R[(nBIN-1)/2-ibin] = ESHits0[ibin+31+15];
420     }
421     }
422    
423     // ---- Effective Energy Deposit Width ---- //
424     double EffWidthSigmaXX = 0.;
425     double EffWidthSigmaYY = 0.;
426     double totalEnergyXX = 0.;
427     double totalEnergyYY = 0.;
428     double EffStatsXX = 0.;
429     double EffStatsYY = 0.;
430     for (int id_X=0; id_X<21; id_X++) {
431     totalEnergyXX += esRH_F[id_X];
432     EffStatsXX += esRH_F[id_X]*(id_X-10)*(id_X-10);
433     totalEnergyYY += esRH_R[id_X];
434     EffStatsYY += esRH_R[id_X]*(id_X-10)*(id_X-10);
435     }
436     EffWidthSigmaXX = (totalEnergyXX>0.) ? sqrt(fabs(EffStatsXX / totalEnergyXX)) : 0.;
437     EffWidthSigmaYY = (totalEnergyYY>0.) ? sqrt(fabs(EffStatsYY / totalEnergyYY)) : 0.;
438    
439     esShape.push_back(EffWidthSigmaXX);
440     esShape.push_back(EffWidthSigmaYY);
441    
442     return esShape;
443     }
444