CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PiZeroDiscriminatorProducer.cc
Go to the documentation of this file.
1 // system include files
2 #include <vector>
3 
4 // user include files
5 
8 
9 // Reconstruction Classes
12 
14 #include <fstream>
15 #include <sstream>
16 
18 
19 // Class for Cluster Shape Algorithm
22 
25 
26 // to compute on-the-fly cluster shapes
27 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
34 
35 
36 using namespace std;
37 using namespace reco;
38 using namespace edm;
40 
42  // use configuration file to setup input/output collection names
43 
44  pPreshowerShapeClustersXToken_ = consumes<PreshowerClusterShapeCollection>(edm::InputTag(ps.getParameter<std::string>("preshClusterShapeProducer"), ps.getParameter<std::string>("preshClusterShapeCollectionX")));
45  pPreshowerShapeClustersYToken_ = consumes<PreshowerClusterShapeCollection>(edm::InputTag(ps.getParameter<std::string>("preshClusterShapeProducer"), ps.getParameter<std::string>("preshClusterShapeCollectionY")));
46 
47  correctedPhotonToken_ = consumes<PhotonCollection>(edm::InputTag(ps.getParameter<string>("corrPhoProducer"), ps.getParameter<string>("correctedPhotonCollection")));
48 
49  barrelRecHitCollection_ = ps.getParameter<edm::InputTag>("barrelRecHitCollection");
50  barrelRecHitCollectionToken_ = consumes<EcalRecHitCollection>(barrelRecHitCollection_);
51  endcapRecHitCollection_ = ps.getParameter<edm::InputTag>("endcapRecHitCollection");
52  endcapRecHitCollectionToken_ = consumes<EcalRecHitCollection>(endcapRecHitCollection_);
53 
54  EScorr_ = ps.getParameter<int>("EScorr");
55 
56  preshNst_ = ps.getParameter<int>("preshPi0Nstrip");
57 
58  preshStripECut_ = ps.getParameter<double>("preshStripEnergyCut");
59 
60  w0_ = ps.getParameter<double>("w0");
61 
62  PhotonPi0DiscriminatorAssociationMap_ = ps.getParameter<string>("Pi0Association");
63 
64  string debugString = ps.getParameter<string>("debugLevel");
65 
66  if (debugString == "DEBUG") debugL_pi0 = pDEBUG;
67  else if (debugString == "INFO") debugL_pi0 = pINFO;
68  else debugL_pi0 = pERROR;
69 
70  string tmpPath = ps.getUntrackedParameter<string>("pathToWeightFiles","RecoEcal/EgammaClusterProducers/data/");
71 
72  presh_pi0_algo = new EndcapPiZeroDiscriminatorAlgo(preshStripECut_, preshNst_, tmpPath.c_str());
73 
74  produces< PhotonPi0DiscriminatorAssociationMap >(PhotonPi0DiscriminatorAssociationMap_);
75 
76 
77  nEvt_ = 0;
78 
79 }
80 
81 
83  delete presh_pi0_algo;
84 }
85 
86 
88 
89  ostringstream ostr; // use this stream for all messages in produce
90 
91  if ( debugL_pi0 <= pDEBUG )
92  cout << "\n PiZeroDiscriminatorProducer: ....... Event " << evt.id() << " with Number = " << nEvt_+1
93  << " is analyzing ....... " << endl << endl;
94 
95  // Get ES clusters in X plane
96  Handle<reco::PreshowerClusterShapeCollection> pPreshowerShapeClustersX;
97  evt.getByToken(pPreshowerShapeClustersXToken_, pPreshowerShapeClustersX);
98  const reco::PreshowerClusterShapeCollection *clustersX = pPreshowerShapeClustersX.product();
99  if ( debugL_pi0 <= pDEBUG ) {
100  cout << "\n PiZeroDiscriminatorProducer: pPreshowerShapeClustersX->size() = " << clustersX->size() << endl;
101  }
102  // Get ES clusters in Y plane
103  Handle<reco::PreshowerClusterShapeCollection> pPreshowerShapeClustersY;
104  evt.getByToken(pPreshowerShapeClustersYToken_, pPreshowerShapeClustersY);
105  const reco::PreshowerClusterShapeCollection *clustersY = pPreshowerShapeClustersY.product();
106  if ( debugL_pi0 <= pDEBUG ) {
107  cout << "\n PiZeroDiscriminatorProducer: pPreshowerShapeClustersY->size() = " << clustersY->size() << endl;
108  }
109  auto_ptr<PhotonPi0DiscriminatorAssociationMap> Pi0Assocs_p(new PhotonPi0DiscriminatorAssociationMap);
110 
112  evt.getByToken( barrelRecHitCollectionToken_, pEBRecHits );
113  const EcalRecHitCollection *ebRecHits = pEBRecHits.product();
114 
116  evt.getByToken( endcapRecHitCollectionToken_, pEERecHits );
117  const EcalRecHitCollection *eeRecHits = pEERecHits.product();
118 
119  ESHandle<CaloGeometry> pGeometry;
120  es.get<CaloGeometryRecord>().get(pGeometry);
121  const CaloGeometry *geometry = pGeometry.product();
122 
123  ESHandle<CaloTopology> pTopology;
124  es.get<CaloTopologyRecord>().get(pTopology);
125  const CaloTopology *topology = pTopology.product();
126 
127  //make cycle over Photon Collection
128  Handle<PhotonCollection> correctedPhotonHandle;
129  evt.getByToken(correctedPhotonToken_, correctedPhotonHandle);
130  const PhotonCollection corrPhoCollection = *(correctedPhotonHandle.product());
131  if ( debugL_pi0 <= pDEBUG ) {
132  cout << " PiZeroDiscriminatorProducer: Photon Collection size : " << corrPhoCollection.size() << endl;
133  }
134 
135  for( PhotonCollection::const_iterator iPho = corrPhoCollection.begin(); iPho != corrPhoCollection.end(); iPho++) {
136 
137  float Phot_En = iPho->energy();
138  float Phot_Et = Phot_En*sin(2*atan(exp(-iPho->eta())));
139  float Phot_eta = iPho->eta();
140  float Phot_phi = iPho->phi();
141  float Phot_R9 = iPho->r9();
142 
143  if ( debugL_pi0 <= pDEBUG ) {
144  cout << " PiZeroDiscriminatorProducer: Photon index : " << iPho - corrPhoCollection.begin()
145  << " with Energy = " << Phot_En
146  << " Et = " << Phot_Et
147  << " ETA = " << Phot_eta
148  << " PHI = " << Phot_phi
149  << " R9 = " << Phot_R9 << endl;
150  }
151  SuperClusterRef it_super = iPho->superCluster();
152 
153  float SC_En = it_super->energy();
154  float SC_Et = SC_En*sin(2*atan(exp(-it_super->eta())));
155  float SC_eta = it_super->eta();
156  float SC_phi = it_super->phi();
157 
158  if ( debugL_pi0 <= pDEBUG ) {
159  cout << " PiZeroDiscriminatorProducer: superE = " << SC_En
160  << " superEt = " << SC_Et
161  << " superETA = " << SC_eta
162  << " superPHI = " << SC_phi << endl;
163  }
164 
165  // New way to get ClusterShape info
166  // Find the entry in the map corresponding to the seed BasicCluster of the SuperCluster
167  // DetId id = it_super->seed()->hitsAndFractions()[0].first;
168 
169  // get on-the-fly the cluster shapes
170 // EcalClusterLazyTools lazyTool( evt, es, barrelRecHitCollection_, endcapRecHitCollection_ );
171 
172  float nnoutput = -1.;
173  if(fabs(SC_eta) >= 1.65 && fabs(SC_eta) <= 2.5) { // Use Preshower region only
174  const GlobalPoint pointSC(it_super->x(),it_super->y(),it_super->z()); // get the centroid of the SC
175  if ( debugL_pi0 <= pDEBUG ) cout << "SC centroind = " << pointSC << endl;
176  double SC_seed_energy = it_super->seed()->energy();
177 
178  const CaloClusterPtr seed = it_super->seed();
179 
180  EcalClusterTools::eMax( *seed, ebRecHits );
181 
182  double SC_seed_Shape_E1 = EcalClusterTools::eMax( *seed, eeRecHits );
183  double SC_seed_Shape_E3x3 = EcalClusterTools::e3x3( *seed, eeRecHits, topology );
184  double SC_seed_Shape_E5x5 = EcalClusterTools::e5x5( *seed, eeRecHits, topology );
185 
186  if ( debugL_pi0 <= pDEBUG ) {
187  cout << "PiZeroDiscriminatorProducer: ( SeedBC_energy, E1, E3x3, E5x5) = "
188  << SC_seed_energy << " "
189  << SC_seed_Shape_E1 << " "
190  << SC_seed_Shape_E3x3 << " "
191  << SC_seed_Shape_E5x5 << endl;
192  }
193 
194 // Get the Preshower 2-planes energy vectors associated with the given SC
195  vector<float> vout_stripE1;
196  vector<float> vout_stripE2;
197  for(reco::PreshowerClusterShapeCollection::const_iterator esClus = clustersX->begin();
198  esClus !=clustersX->end(); esClus++) {
199  SuperClusterRef sc_ref = esClus->superCluster();
200  float dR = sqrt((SC_eta-sc_ref->eta())*(SC_eta-sc_ref->eta()) +
201  (SC_phi-sc_ref->phi())*(SC_phi-sc_ref->phi()));
202  if(dR < 0.01 ) {
203 
204  vout_stripE1 = esClus->getStripEnergies();
205 
206  }
207  }
208  for(reco::PreshowerClusterShapeCollection::const_iterator esClus = clustersY->begin();
209  esClus !=clustersY->end(); esClus++) {
210  SuperClusterRef sc_ref = esClus->superCluster();
211  float dR = sqrt((SC_eta-sc_ref->eta())*(SC_eta-sc_ref->eta()) +
212  (SC_phi-sc_ref->phi())*(SC_phi-sc_ref->phi()));
213  if(dR < 0.01 ) {
214 
215  vout_stripE2 = esClus->getStripEnergies();
216 
217  }
218  }
219 
220  if(vout_stripE1.size() == 0 || vout_stripE2.size() == 0 ) {
221  if ( debugL_pi0 <= pDEBUG )
222  cout << " PiZeroDiscriminatorProducer: Attention!!!!! Not Valid ES NN input Variables Return NNout = -1" << endl;
223  Pi0Assocs_p->insert(Ref<PhotonCollection>(correctedPhotonHandle,iPho - corrPhoCollection.begin()), nnoutput);
224  continue;
225  }
226 
227  if ( debugL_pi0 <= pDEBUG ) {
228  cout << "PiZeroDiscriminatorProducer : vout_stripE1.size = " << vout_stripE1.size()
229  << " vout_stripE2.size = " << vout_stripE2.size() << endl;
230  cout << "PiZeroDiscriminatorProducer : ES_input_vector = " ;
231  for(int k1=0;k1<11;k1++) {
232  cout << vout_stripE1[k1] << " " ;
233  }
234  for(int k1=0;k1<11;k1++) {
235  cout << vout_stripE2[k1] << " " ;
236  }
237  cout << endl;
238  }
239 
240  bool valid_NNinput = presh_pi0_algo->calculateNNInputVariables(vout_stripE1, vout_stripE2,
241  SC_seed_Shape_E1, SC_seed_Shape_E3x3, SC_seed_Shape_E5x5, EScorr_);
242 
243  if(!valid_NNinput) {
244  if ( debugL_pi0 <= pDEBUG )
245  cout << " PiZeroDiscriminatorProducer: Attention!!!!! Not Valid ES NN input Variables Return NNout = -1" << endl;
246  Pi0Assocs_p->insert(Ref<PhotonCollection>(correctedPhotonHandle,iPho - corrPhoCollection.begin()), nnoutput);
247  continue;
248  }
249 
250  float* nn_input_var = presh_pi0_algo->get_input_vector();
251 
252  if ( debugL_pi0 <= pDEBUG ) {
253  cout << " PiZeroDiscriminatorProducer: NN_ESEndcap_input_vector+Et+Eta+Phi+R9 = " ;
254  for(int k1=0;k1<25;k1++) {
255  cout << nn_input_var[k1] << " " ;
256  }
257  cout << SC_Et << " " << SC_eta << " " << SC_phi << " " << Phot_R9 << endl;
258  }
259 
260  nnoutput = presh_pi0_algo->GetNNOutput(SC_Et);
261 
262  if ( debugL_pi0 <= pDEBUG ) {
263  cout << " PiZeroDiscriminatorProducer: Event : " << evt.id()
264  << " SC id = " << iPho - corrPhoCollection.begin()
265  << " with Pt = " << SC_Et
266  << " eta = " << SC_eta
267  << " phi = " << SC_phi
268  << " contains: " << it_super->clustersSize() << " BCs "
269  << " has NNout = " << nnoutput << endl;
270  }
271 
272  Pi0Assocs_p->insert(Ref<PhotonCollection>(correctedPhotonHandle,iPho - corrPhoCollection.begin()), nnoutput);
273 
274  } else if((fabs(SC_eta) <= 1.4442) || (fabs(SC_eta) < 1.65 && fabs(SC_eta) >= 1.566) || fabs(SC_eta) >= 2.5) {
275 
276  const CaloClusterPtr seed = it_super->seed();
277 
278  double SC_seed_Shape_E1 = EcalClusterTools::eMax( *seed, ebRecHits );
279  double SC_seed_Shape_E3x3 = EcalClusterTools::e3x3( *seed, ebRecHits, topology );
280  double SC_seed_Shape_E5x5 = EcalClusterTools::e5x5( *seed, ebRecHits, topology );
281  double SC_seed_Shape_E2 = EcalClusterTools::e2nd( *seed, ebRecHits );
282 
283  std::vector<float> vCov = EcalClusterTools::covariances( *seed, ebRecHits , topology, geometry, w0_ );
284 
285  double SC_seed_Shape_cEE = vCov[0];
286  double SC_seed_Shape_cEP = vCov[1];
287  double SC_seed_Shape_cPP = vCov[2];
288 
289  double SC_seed_Shape_E2x2 = EcalClusterTools::e2x2( *seed, ebRecHits, topology );
290  double SC_seed_Shape_E3x2 = EcalClusterTools::e3x2( *seed, ebRecHits, topology );
291 
292  double SC_seed_Shape_E3x2r = 0.0;
293  double SC_seed_Shape_ELeft = EcalClusterTools::eLeft( *seed, ebRecHits, topology );
294  double SC_seed_Shape_ERight = EcalClusterTools::eRight( *seed, ebRecHits, topology );
295  double SC_seed_Shape_ETop = EcalClusterTools::eTop( *seed, ebRecHits, topology );
296  double SC_seed_Shape_EBottom = EcalClusterTools::eBottom( *seed, ebRecHits, topology );
297 
298  double DA = SC_seed_Shape_E2x2 - SC_seed_Shape_E2 - SC_seed_Shape_E1;
299 
300  if(SC_seed_Shape_E2==SC_seed_Shape_ETop || SC_seed_Shape_E2==SC_seed_Shape_EBottom) {
301  if( SC_seed_Shape_ELeft > SC_seed_Shape_ERight ) {
302  SC_seed_Shape_E3x2r = (DA - SC_seed_Shape_ELeft)/(0.25+SC_seed_Shape_ELeft);
303  } else {
304  SC_seed_Shape_E3x2r = (DA - SC_seed_Shape_ERight)/(0.25+SC_seed_Shape_ERight);
305  }
306 
307  } else if(SC_seed_Shape_E2==SC_seed_Shape_ELeft || SC_seed_Shape_E2==SC_seed_Shape_ERight) {
308 
309  if( SC_seed_Shape_ETop > SC_seed_Shape_EBottom ) {
310  SC_seed_Shape_E3x2r = (DA - SC_seed_Shape_ETop)/(0.25+SC_seed_Shape_ETop);
311  } else {
312  SC_seed_Shape_E3x2r = (DA - SC_seed_Shape_EBottom)/(0.25+SC_seed_Shape_EBottom);
313  }
314 
315  }
316 
317  double SC_seed_Shape_xcog = EcalClusterTools::eRight( *seed, ebRecHits, topology ) - EcalClusterTools::e2x5Left( *seed, ebRecHits, topology );
318  double SC_seed_Shape_ycog = EcalClusterTools::e2x5Top( *seed, ebRecHits, topology ) - EcalClusterTools::e2x5Bottom( *seed, ebRecHits, topology );
319 
320 
321  if ( debugL_pi0 <= pDEBUG ) {
322  cout << "PiZeroDiscriminatorProduce: lazyTool (E1,E3x3,E5x5,E2,cEE,cEP,cPP,E2x2,E3x2_E3x2r,Xcog,Ycog,E2x5Bottom,E2x5Top,Et,Eta,PhiR9) = ( "
323  << SC_seed_Shape_E1 << " "
324  << SC_seed_Shape_E3x3 << " "
325  << SC_seed_Shape_E5x5 << " "
326  << SC_seed_Shape_E2 << " "
327  << SC_seed_Shape_cEE << " "
328  << SC_seed_Shape_cEP << " "
329  << SC_seed_Shape_cPP << " "
330  << SC_seed_Shape_E2x2 << " "
331  << SC_seed_Shape_E3x2 << " "
332  << SC_seed_Shape_E3x2r << " "
333  << SC_seed_Shape_xcog << " "
334  << SC_seed_Shape_ycog << " "
335  << EcalClusterTools::e2x5Bottom( *seed, ebRecHits, topology ) << " "
336  << EcalClusterTools::e2x5Top( *seed, ebRecHits, topology ) << " "
337  << SC_Et << " "
338  << SC_eta << " "
339  << SC_phi << " "
340  << Phot_R9 << " )" << endl;
341  }
342 
343  float SC_et = it_super->energy()*sin(2*atan(exp(-it_super->eta())));
344 
345  presh_pi0_algo->calculateBarrelNNInputVariables(SC_et, SC_seed_Shape_E1, SC_seed_Shape_E3x3,
346  SC_seed_Shape_E5x5, SC_seed_Shape_E2,
347  SC_seed_Shape_cEE, SC_seed_Shape_cEP,
348  SC_seed_Shape_cPP, SC_seed_Shape_E2x2,
349  SC_seed_Shape_E3x2, SC_seed_Shape_E3x2r,
350  SC_seed_Shape_xcog, SC_seed_Shape_ycog);
351 
352  float* nn_input_var = presh_pi0_algo->get_input_vector();
353 
354  if ( debugL_pi0 <= pDEBUG ) {
355  cout << " PiZeroDiscriminatorProducer : NN_barrel_nonESEndcap_variables+Et+Eta+Phi+R9 = " ;
356  for(int k3=0;k3<12;k3++) {
357  cout << nn_input_var[k3] << " " ;
358  }
359  cout << SC_Et << " " << SC_eta << " " << SC_phi << " " << Phot_R9 << endl;
360 
361  }
362 
363  nnoutput = presh_pi0_algo->GetBarrelNNOutput(SC_et);
364 
365 
366  if ( debugL_pi0 <= pDEBUG ) {
367  cout << "PiZeroDiscriminatorProducer : Event : " << evt.id()
368  << " SC id = " << iPho - corrPhoCollection.begin()
369  << " with Pt = " << SC_Et
370  << " eta = " << SC_eta
371  << " phi = " << SC_phi
372  << " contains: " << it_super->clustersSize() << " BCs "
373  << " has NNout = " << nnoutput
374  << endl;
375  }
376 
377  Pi0Assocs_p->insert(Ref<PhotonCollection>(correctedPhotonHandle,iPho - corrPhoCollection.begin()), nnoutput);
378  } else { Pi0Assocs_p->insert(Ref<PhotonCollection>(correctedPhotonHandle,iPho - corrPhoCollection.begin()), -1.);}
379  } // end of cycle over Photons
380 
381  evt.put(Pi0Assocs_p,PhotonPi0DiscriminatorAssociationMap_);
382  if ( debugL_pi0 <= pDEBUG ) cout << "PiZeroDiscriminatorProducer: PhotonPi0DiscriminatorAssociationMap added to the event" << endl;
383 
384  nEvt_++;
385 
386  LogDebug("PiZeroDiscriminatorDebug") << ostr.str();
387 
388 
389 }
#define LogDebug(id)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
CaloTopology const * topology(0)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:113
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< PreshowerClusterShape > PreshowerClusterShapeCollection
collection of PreshowerClusterShape objects
T const * product() const
Definition: Handle.h:81
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:86
ESHandle< TrackerGeometry > geometry
edm::EventID id() const
Definition: EventBase.h:56
tuple cout
Definition: gather_cfg.py:121
PiZeroDiscriminatorProducer(const edm::ParameterSet &ps)
virtual void produce(edm::Event &evt, const edm::EventSetup &es)