CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PiZeroAnalyzer.cc
Go to the documentation of this file.
1 #include <iostream>
2 //
3 
5 
6 
7 //#define TWOPI 6.283185308
8 //
9 
22 using namespace std;
23 
24 
26 {
27 
28  fName_ = pset.getUntrackedParameter<std::string>("Name");
29  verbosity_ = pset.getUntrackedParameter<int>("Verbosity");
30 
31  prescaleFactor_ = pset.getUntrackedParameter<int>("prescaleFactor",1);
32 
33 
34 
35  barrelEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.getParameter<edm::InputTag>("barrelEcalHits"));
36  endcapEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.getParameter<edm::InputTag>("endcapEcalHits"));
37 
38 
39  standAlone_ = pset.getParameter<bool>("standAlone");
40 
41 
42 
43 
44 
45  // parameters for Pizero finding
46  seleXtalMinEnergy_ = pset.getParameter<double> ("seleXtalMinEnergy");
47  clusSeedThr_ = pset.getParameter<double> ("clusSeedThr");
48  clusEtaSize_ = pset.getParameter<int> ("clusEtaSize");
49  clusPhiSize_ = pset.getParameter<int> ("clusPhiSize");
50  ParameterLogWeighted_ = pset.getParameter<bool> ("ParameterLogWeighted");
51  ParameterX0_ = pset.getParameter<double> ("ParameterX0");
52  ParameterT0_barl_ = pset.getParameter<double> ("ParameterT0_barl");
53  ParameterW0_ = pset.getParameter<double> ("ParameterW0");
54 
55  selePtGammaOne_ = pset.getParameter<double> ("selePtGammaOne");
56  selePtGammaTwo_ = pset.getParameter<double> ("selePtGammaTwo");
57  seleS4S9GammaOne_ = pset.getParameter<double> ("seleS4S9GammaOne");
58  seleS4S9GammaTwo_ = pset.getParameter<double> ("seleS4S9GammaTwo");
59  selePtPi0_ = pset.getParameter<double> ("selePtPi0");
60  selePi0Iso_ = pset.getParameter<double> ("selePi0Iso");
61  selePi0BeltDR_ = pset.getParameter<double> ("selePi0BeltDR");
62  selePi0BeltDeta_ = pset.getParameter<double> ("selePi0BeltDeta");
63  seleMinvMaxPi0_ = pset.getParameter<double> ("seleMinvMaxPi0");
64  seleMinvMinPi0_ = pset.getParameter<double> ("seleMinvMinPi0");
65 
66  parameters_ = pset;
67 
68 
69 }
70 
71 
72 
74 
75 
76 
77 
78 }
79 
80 
82 {
83 
84 
85  nEvt_=0;
86  nEntry_=0;
87 
88  dbe_ = 0;
90 
91 
92 
93  if (dbe_) {
94  if (verbosity_ > 0 ) {
95  dbe_->setVerbose(1);
96  } else {
97  dbe_->setVerbose(0);
98  }
99  }
100  if (dbe_) {
101  if (verbosity_ > 0 ) dbe_->showDirStructure();
102  }
103 
104 
105 
106 
107  //booking all histograms
108 
109  if (dbe_) {
110 
111  currentFolder_.str("");
112  currentFolder_ << "Egamma/PiZeroAnalyzer/";
113  dbe_->setCurrentFolder(currentFolder_.str());
114 
115 
116 
117 
118  hMinvPi0EB_ = dbe_->book1D("Pi0InvmassEB","Pi0 Invariant Mass in EB",100,0.,0.5);
119  hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ",1);
120 
121  hPt1Pi0EB_ = dbe_->book1D("Pt1Pi0EB","Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
122  hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ",1);
123 
124  hPt2Pi0EB_ = dbe_->book1D("Pt2Pi0EB","Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
125  hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ",1);
126 
127  hPtPi0EB_ = dbe_->book1D("PtPi0EB","Pi0 Pt in EB",100,0.,20.);
128  hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ",1);
129 
130  hIsoPi0EB_ = dbe_->book1D("IsoPi0EB","Pi0 Iso in EB",50,0.,1.);
131  hIsoPi0EB_->setAxisTitle("Pi0 Iso",1);
132 
133 
134  }
135 
136 }
137 
138 
139 
140 
141 
142 
144 {
145 
146  using namespace edm;
147 
148  if (nEvt_% prescaleFactor_ ) return;
149  nEvt_++;
150  LogInfo("PiZeroAnalyzer") << "PiZeroAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
151 
152 
153  // Get EcalRecHits
154  bool validEcalRecHits=true;
155  Handle<EcalRecHitCollection> barrelHitHandle;
156  EcalRecHitCollection barrelRecHits;
157  e.getByToken(barrelEcalHits_token_, barrelHitHandle);
158  if (!barrelHitHandle.isValid()) {
159  edm::LogError("PhotonProducer") << "Error! Can't get the product: barrelEcalHits_token_";
160  validEcalRecHits=false;
161  }
162 
163  Handle<EcalRecHitCollection> endcapHitHandle;
164  e.getByToken(endcapEcalHits_token_, endcapHitHandle);
165  EcalRecHitCollection endcapRecHits;
166  if (!endcapHitHandle.isValid()) {
167  edm::LogError("PhotonProducer") << "Error! Can't get the product: endcapEcalHits_token_";
168  validEcalRecHits=false;
169  }
170 
171  if (validEcalRecHits) makePizero(esup, barrelHitHandle, endcapHitHandle);
172 
173 
174 
175 }
176 
178 
179  const EcalRecHitCollection *hitCollection_p = rhEB.product();
180 
181  edm::ESHandle<CaloGeometry> geoHandle;
182  es.get<CaloGeometryRecord>().get(geoHandle);
183 
184  edm::ESHandle<CaloTopology> theCaloTopology;
185  es.get<CaloTopologyRecord>().get(theCaloTopology);
186 
187 
188  const CaloSubdetectorGeometry *geometry_p;
189  const CaloSubdetectorTopology *topology_p;
190  const CaloSubdetectorGeometry *geometryES_p;
191  geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
192  geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
193 
194  // Parameters for the position calculation:
195  edm::ParameterSet posCalcParameters =
196  parameters_.getParameter<edm::ParameterSet>("posCalcParameters");
197  PositionCalc posCalculator_ = PositionCalc(posCalcParameters);
198  //
199  std::map<DetId, EcalRecHit> recHitsEB_map;
200  //
201  std::vector<EcalRecHit> seeds;
202 
203  seeds.clear();
204  //
205  vector<EBDetId> usedXtals;
206  usedXtals.clear();
207  //
209  //
210  static const int MAXCLUS = 2000;
211  int nClus=0;
212  vector<float> eClus;
213  vector<float> etClus;
214  vector<float> etaClus;
215  vector<float> phiClus;
216  vector<EBDetId> max_hit;
217  vector< vector<EcalRecHit> > RecHitsCluster;
218  vector<float> s4s9Clus;
219 
220  // find cluster seeds in EB
221  for(itb=rhEB->begin(); itb!=rhEB->end(); ++itb){
222  EBDetId id(itb->id());
223  double energy = itb->energy();
224  if (energy > seleXtalMinEnergy_) {
225  std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
226  recHitsEB_map.insert(map_entry);
227  }
228  if (energy > clusSeedThr_) seeds.push_back(*itb);
229  } // Eb rechits
230 
231  sort(seeds.begin(), seeds.end(), ecalRecHitLess());
232  for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
233  EBDetId seed_id = itseed->id();
234  std::vector<EBDetId>::const_iterator usedIds;
235 
236  bool seedAlreadyUsed=false;
237  for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
238  if(*usedIds==seed_id){
239  seedAlreadyUsed=true;
240  //cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<endl;
241  break;
242  }
243  }
244  if(seedAlreadyUsed)continue;
245  topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel);
246  std::vector<DetId> clus_v = topology_p->getWindow(seed_id,clusEtaSize_,clusPhiSize_);
247  //std::vector<DetId> clus_used;
248  std::vector<std::pair<DetId, float> > clus_used;
249 
250  vector<EcalRecHit> RecHitsInWindow;
251 
252  double simple_energy = 0;
253 
254  for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
255  // EBDetId EBdet = *det;
256  // cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<endl;
257  bool HitAlreadyUsed=false;
258  for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
259  if(*usedIds==*det){
260  HitAlreadyUsed=true;
261  break;
262  }
263  }
264  if(HitAlreadyUsed)continue;
265  if (recHitsEB_map.find(*det) != recHitsEB_map.end()){
266  // cout<<" Used det "<< EBdet<<endl;
267  std::map<DetId, EcalRecHit>::iterator aHit;
268  aHit = recHitsEB_map.find(*det);
269  usedXtals.push_back(*det);
270  RecHitsInWindow.push_back(aHit->second);
271  clus_used.push_back( std::pair<DetId, float>(*det, 1.) );
272  simple_energy = simple_energy + aHit->second.energy();
273  }
274  }
275 
276  math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_p,geometry_p,geometryES_p);
277  float theta_s = 2. * atan(exp(-clus_pos.eta()));
278  float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
279  float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
280  // float p0z_s = simple_energy * cos(theta_s);
281  float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
282 
283  //cout << " Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
284 
285  eClus.push_back(simple_energy);
286  etClus.push_back(et_s);
287  etaClus.push_back(clus_pos.eta());
288  phiClus.push_back(clus_pos.phi());
289  max_hit.push_back(seed_id);
290  RecHitsCluster.push_back(RecHitsInWindow);
291  //Compute S4/S9 variable
292  //We are not sure to have 9 RecHits so need to check eta and phi:
293  float s4s9_[4];
294  for(int i=0;i<4;i++)s4s9_[i]= itseed->energy();
295  for(unsigned int j=0; j<RecHitsInWindow.size();j++){
296  //cout << " Simple cluster rh, ieta, iphi : "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" "<<((EBDetId)RecHitsInWindow[j].id()).iphi()<<endl;
297  if((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()-1 && seed_id.ieta()!=1 ) || ( seed_id.ieta()==1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()-2))){
298  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
299  s4s9_[0]+=RecHitsInWindow[j].energy();
300  }else{
301  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()){
302  s4s9_[0]+=RecHitsInWindow[j].energy();
303  s4s9_[1]+=RecHitsInWindow[j].energy();
304  }else{
305  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
306  s4s9_[1]+=RecHitsInWindow[j].energy();
307  }
308  }
309  }
310  }else{
311  if(((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()){
312  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
313  s4s9_[0]+=RecHitsInWindow[j].energy();
314  s4s9_[3]+=RecHitsInWindow[j].energy();
315  }else{
316  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
317  s4s9_[1]+=RecHitsInWindow[j].energy();
318  s4s9_[2]+=RecHitsInWindow[j].energy();
319  }
320  }
321  }else{
322  if((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()+1 && seed_id.ieta()!=-1 ) || ( seed_id.ieta()==-1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()+2))){
323  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
324  s4s9_[3]+=RecHitsInWindow[j].energy();
325  }else{
326  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()){
327  s4s9_[2]+=RecHitsInWindow[j].energy();
328  s4s9_[3]+=RecHitsInWindow[j].energy();
329  }else{
330  if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
331  s4s9_[2]+=RecHitsInWindow[j].energy();
332  }
333  }
334  }
335  }else{
336  cout<<" (EBDetId)RecHitsInWindow[j].id()).ieta() "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" seed_id.ieta() "<<seed_id.ieta()<<endl;
337  cout<<" Problem with S4 calculation "<<endl;return;
338  }
339  }
340  }
341  }
342  s4s9Clus.push_back(*max_element( s4s9_,s4s9_+4)/simple_energy);
343  // cout<<" s4s9Clus[0] "<<s4s9_[0]/simple_energy<<" s4s9Clus[1] "<<s4s9_[1]/simple_energy<<" s4s9Clus[2] "<<s4s9_[2]/simple_energy<<" s4s9Clus[3] "<<s4s9_[3]/simple_energy<<endl;
344  // cout<<" Max "<<*max_element( s4s9_,s4s9_+4)/simple_energy<<endl;
345  nClus++;
346  if (nClus == MAXCLUS) return;
347  } // End loop over seed clusters
348 
349  // cout<< " Pi0 clusters: "<<nClus<<endl;
350 
351  // Selection, based on Simple clustering
352  //pi0 candidates
353  static const int MAXPI0S = 200;
354  int npi0_s=0;
355 
356  vector<EBDetId> scXtals;
357  scXtals.clear();
358 
359  if (nClus <= 1) return;
360  for(Int_t i=0 ; i<nClus ; i++){
361  for(Int_t j=i+1 ; j<nClus ; j++){
362  // cout<<" i "<<i<<" etClus[i] "<<etClus[i]<<" j "<<j<<" etClus[j] "<<etClus[j]<<endl;
363  if( etClus[i]>selePtGammaOne_ && etClus[j]>selePtGammaTwo_ && s4s9Clus[i]>seleS4S9GammaOne_ && s4s9Clus[j]>seleS4S9GammaTwo_){
364  float theta_0 = 2. * atan(exp(-etaClus[i]));
365  float theta_1 = 2. * atan(exp(-etaClus[j]));
366 
367  float p0x = eClus[i] * sin(theta_0) * cos(phiClus[i]);
368  float p1x = eClus[j] * sin(theta_1) * cos(phiClus[j]);
369  float p0y = eClus[i] * sin(theta_0) * sin(phiClus[i]);
370  float p1y = eClus[j] * sin(theta_1) * sin(phiClus[j]);
371  float p0z = eClus[i] * cos(theta_0);
372  float p1z = eClus[j] * cos(theta_1);
373 
374  float pt_pi0 = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
375  // cout<<" pt_pi0 "<<pt_pi0<<endl;
376  if (pt_pi0 < selePtPi0_)continue;
377  float m_inv = sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
378  if ( (m_inv<seleMinvMaxPi0_) && (m_inv>seleMinvMinPi0_) ){
379 
380  //New Loop on cluster to measure isolation:
381  vector<int> IsoClus;
382  IsoClus.clear();
383  float Iso = 0;
384  TVector3 pi0vect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
385  for(Int_t k=0 ; k<nClus ; k++){
386  if(k==i || k==j)continue;
387  TVector3 Clusvect = TVector3(eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * cos(phiClus[k]), eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * sin(phiClus[k]) , eClus[k] * cos(2. * atan(exp(-etaClus[k]))));
388  float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
389  float drclpi0 = Clusvect.DeltaR(pi0vect);
390 
391  if((drclpi0<selePi0BeltDR_) && (dretaclpi0<selePi0BeltDeta_) ){
392 
393  Iso = Iso + etClus[k];
394  IsoClus.push_back(k);
395  }
396  }
397 
398 
399  if(Iso/pt_pi0<selePi0Iso_){
400 
401  hMinvPi0EB_->Fill(m_inv);
402  hPt1Pi0EB_->Fill(etClus[i]);
403  hPt2Pi0EB_->Fill(etClus[j]);
404  hPtPi0EB_->Fill(pt_pi0);
405  hIsoPi0EB_->Fill(Iso/pt_pi0);
406 
407 
408  npi0_s++;
409  }
410 
411  if(npi0_s == MAXPI0S) return;
412  }
413  }
414  }
415  }
416 
417 }
418 
419 
420 
422 {
423 
424 
425 
426  bool outputMEsInRootFile = parameters_.getParameter<bool>("OutputMEsInRootFile");
427  std::string outputFileName = parameters_.getParameter<std::string>("OutputFileName");
428  if(outputMEsInRootFile){
429  dbe_->save(outputFileName);
430  }
431 
432  edm::LogInfo("PiZeroAnalyzer") << "Analyzed " << nEvt_ << "\n";
433  return ;
434 }
435 
436 
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:954
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
PiZeroAnalyzer(const edm::ParameterSet &)
void makePizero(const edm::EventSetup &es, const edm::Handle< EcalRecHitCollection > eb, const edm::Handle< EcalRecHitCollection > ee)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< EcalRecHit >::const_iterator const_iterator
return((rh^lh)&mask)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
T sqrt(T t)
Definition: SSEVec.h:48
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
int j
Definition: DBlmapReader.cc:9
void setVerbose(unsigned level)
Definition: DQMStore.cc:631
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
DQMStore * dbe_
virtual void endJob()
int k[5][pyjets_maxn]
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
virtual void analyze(const edm::Event &, const edm::EventSetup &)
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: Handle.h:81
edm::EventID id() const
Definition: EventBase.h:56
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
Definition: DQMStore.cc:2540
virtual ~PiZeroAnalyzer()
tuple cout
Definition: gather_cfg.py:121
void showDirStructure(void) const
Definition: DQMStore.cc:3332
virtual void beginJob()
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:667