CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GenXSecAnalyzer.cc
Go to the documentation of this file.
3 #include <iostream>
4 #include <iomanip>
5 
7  hepidwtup_(-1),
8  theProcesses_size(0),
9  hasHepMCFilterInfo_(false),
10  xsec_(0,0),
11  filterOnlyEffStat_(0,0,0,0,0.,0.,0.,0.),
12  hepMCFilterEffStat_(0,0,0,0,0.,0.,0.,0.)
13 {
14  eventEffStat_.clear();
15  jetMatchEffStat_.clear();
16  xsecBeforeMatching_.clear();
17  xsecAfterMatching_.clear();
18  products_.clear();
19 }
20 
22 {
23 }
24 
25 void
27  eventEffStat_.clear();
28  jetMatchEffStat_.clear();
29  xsecBeforeMatching_.clear();
30  xsecAfterMatching_.clear();
31  products_.clear();
32 }
33 
34 void
36 {
37 }
38 
39 // ------------ method called once each job just after ending the event loop ------------
40 
41 void
43  // add up information of GenFilterInfo from different luminosity blocks
45  if(iLumi.getByLabel("genFilterEfficiencyProducer", genFilter))
46  filterOnlyEffStat_.mergeProduct(*genFilter);
47 
48  edm::Handle<GenFilterInfo> hepMCFilter;
49  hasHepMCFilterInfo_ = iLumi.getByLabel("generator",hepMCFilter);
51  hepMCFilterEffStat_.mergeProduct(*hepMCFilter);
52 
54  bool foundLumiInfo = iLumi.getByLabel("generator",genLumiInfo);
55  if (!foundLumiInfo) return;
56 
57  hepidwtup_ = genLumiInfo->getHEPIDWTUP();
58 
59  std::vector<GenLumiInfoProduct::ProcessInfo> theProcesses = genLumiInfo->getProcessInfos();
60  theProcesses_size = theProcesses.size();
61  // if it's a pure parton-shower generator, check there should be only one element in thisProcessInfos
62  // the error of lheXSec is -1
63  if(hepidwtup_== -1)
64  {
65  if(theProcesses_size!=1){
66  edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "Pure parton shower has thisProcessInfos size!=1";
67  return;
68  }
69  if(theProcesses[0].lheXSec().value()<1e-6){
70  edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "cross section value = " << theProcesses[0].lheXSec().value();
71  return;
72  }
73  }
74 
75  // now determine if this LuminosityBlock has the same lheXSec as existing products
76  bool sameMC = false;
77  for(unsigned int i=0; i < products_.size(); i++){
78 
79  if(products_[i].mergeProduct(*genLumiInfo))
80  sameMC = true;
81  else if(!products_[i].samePhysics(*genLumiInfo))
82  {
83  edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "Merging samples that come from different physics processes";
84  return;
85  }
86 
87  }
88 
89  if(!sameMC)
90  products_.push_back(*genLumiInfo);
91 
92 
93  // initialize jetMatchEffStat with nprocess+1 elements
94  if(jetMatchEffStat_.size()==0){
95 
96  for(unsigned int ip=0; ip < theProcesses_size; ip++)
97  {
98  jetMatchEffStat_.push_back(GenFilterInfo(0,0,0,0,0.,0.,0.,0.));
99  }
100  jetMatchEffStat_.push_back(GenFilterInfo(0,0,0,0,0.,0.,0.,0.));
101  }
102 
103  // initialize event-level statistics with nprocess+1 elements
104  if(eventEffStat_.size()==0){
105 
106  for(unsigned int ip=0; ip < theProcesses_size; ip++)
107  {
108  eventEffStat_.push_back(GenFilterInfo(0,0,0,0,0.,0.,0.,0.));
109  }
110  eventEffStat_.push_back(GenFilterInfo(0,0,0,0,0.,0.,0.,0.));
111  }
112 
113  // doing generic summing for jet matching statistics
114  for(unsigned int ip=0; ip < theProcesses_size; ip++)
115  {
116  GenLumiInfoProduct::FinalStat temp_killed = theProcesses[ip].killed();
117  GenLumiInfoProduct::FinalStat temp_selected = theProcesses[ip].selected();
118  double passw = temp_killed.sum();
119  double passw2 = temp_killed.sum2();
120  double totalw = temp_selected.sum();
121  double totalw2 = temp_selected.sum2();
122  // matching statistics for each process
123  jetMatchEffStat_[ip].mergeProduct(GenFilterInfo(
124  theProcesses[ip].nPassPos(),
125  theProcesses[ip].nPassNeg(),
126  theProcesses[ip].nTotalPos(),
127  theProcesses[ip].nTotalNeg(),
128  passw,
129  passw2,
130  totalw,
131  totalw2)
132  );
133  // matching statistics for all processes
135  theProcesses[ip].nPassPos(),
136  theProcesses[ip].nPassNeg(),
137  theProcesses[ip].nTotalPos(),
138  theProcesses[ip].nTotalNeg(),
139  passw,
140  passw2,
141  totalw,
142  totalw2)
143  );
144 
145 
146 
147  // event-level statistics for each process
148  eventEffStat_[ip].mergeProduct(GenFilterInfo(
149  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
150  0,
151  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
152  0,
153  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
154  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
155  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
156  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
157  );
158  // event-level statistics for all processes
160  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
161  0,
162  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
163  0,
164  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
165  theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
166  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
167  theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
168  );
169 
170 
171  }
172 
173 
174  return;
175 
176 }
177 
178 void
180 {
181  // for pure parton shower generator
182 
183  if(hepidwtup_== -1)
184  {
185  double sigSum = 0.0;
186  double totalN = 0.0;
187  for(unsigned int i=0; i < products_.size(); i++){
188 
189  GenLumiInfoProduct::ProcessInfo proc = products_[i].getProcessInfos()[0];
190  double hepxsec_value = proc.lheXSec().value();
191  sigSum += proc.tried().sum() * hepxsec_value;
192  totalN += proc.tried().sum();
193  }
194  // average over number of events since the cross sections have no errors
195  double sigAve = totalN>1e-6? sigSum/totalN: 0;
196  xsecBeforeMatching_.push_back(GenLumiInfoProduct::XSec(sigAve,-1));
197  xsecAfterMatching_.push_back(GenLumiInfoProduct::XSec(sigAve,-1));
198  }
199  // for ME+parton shower MC
200  else{
201 
202  const unsigned int sizeOfInfos= theProcesses_size+1;
203  // for computing cross sectiona before matching
204  double sum_numerator_before[sizeOfInfos];
205  double sum_denominator_before[sizeOfInfos];
206 
207  // for computing cross sectiona after matching
208  double sum_numerator_after[sizeOfInfos];
209  double sum_denominator_after[sizeOfInfos];
210 
211  // initialize every element with zero
212  for(unsigned int i=0; i < sizeOfInfos; i++)
213  {
214  sum_numerator_before[i]=0;
215  sum_denominator_before[i]=0;
216  sum_numerator_after[i]=0;
217  sum_denominator_after[i]=0;
218  }
219 
220  // loop over different MC samples
221  for(unsigned int i=0; i < products_.size(); i++){
222 
223  // sum of cross sections and errors over different processes
224  double sigSelSum = 0.0;
225  double errSel2Sum = 0.0;
226  double sigSum = 0.0;
227  double err2Sum = 0.0;
228 
229  // loop over different processes for each sample
230  unsigned int vectorSize = products_[i].getProcessInfos().size();
231  for(unsigned int ip=0; ip < vectorSize; ip++){
232  GenLumiInfoProduct::ProcessInfo proc = products_[i].getProcessInfos()[ip];
233  double hepxsec_value = proc.lheXSec().value();
234  double hepxsec_error = proc.lheXSec().error();
235 
236  // skips computation if jet matching efficiency=0
237  if (proc.killed().n()<1)
238  continue;
239 
240  // computing jet matching efficiency for this process
241  double fracAcc = 0;
242  double ntotal = proc.nTotalPos()-proc.nTotalNeg();
243  double npass = proc.nPassPos() -proc.nPassNeg();
244  switch(hepidwtup_){
245  case 3: case -3:
246  fracAcc = ntotal > 1e-6? npass/ntotal: -1;
247  break;
248  default:
249  fracAcc = proc.selected().sum() > 1e-6? proc.killed().sum() / proc.selected().sum():-1;
250  break;
251  }
252 
253  if(fracAcc<1e-6)continue;
254 
255  // cross section after matching for this particular process
256  double sigmaFin = hepxsec_value * fracAcc;
257 
258  // computing error on jet matching efficiency
259  double relErr = 1.0;
260  double efferr2=0;
261  switch(hepidwtup_) {
262  case 3: case -3:
263  {
264  double ntotal_pos = proc.nTotalPos();
265  double effp = ntotal_pos > 1e-6?
266  (double)proc.nPassPos()/ntotal_pos:0;
267  double effp_err2 = ntotal_pos > 1e-6?
268  (1-effp)*effp/ntotal_pos: 0;
269 
270  double ntotal_neg = proc.nTotalNeg();
271  double effn = ntotal_neg > 1e-6?
272  (double)proc.nPassNeg()/ntotal_neg:0;
273  double effn_err2 = ntotal_neg > 1e-6?
274  (1-effn)*effn/ntotal_neg: 0;
275 
276  efferr2 = ntotal > 0 ?
277  (ntotal_pos*ntotal_pos*effp_err2 +
278  ntotal_neg*ntotal_neg*effn_err2)/ntotal/ntotal:0;
279  break;
280  }
281  default:
282  {
283  double denominator = pow(proc.selected().sum(),4);
284  double passw = proc.killed().sum();
285  double passw2 = proc.killed().sum2();
286  double failw = proc.selected().sum() - passw;
287  double failw2 = proc.selected().sum2() - passw2;
288  double numerator = (passw2*failw*failw + failw2*passw*passw);
289 
290  efferr2 = denominator>1e-6?
291  numerator/denominator:0;
292  break;
293  }
294  }
295  double delta2Veto = efferr2/fracAcc/fracAcc;
296 
297  // computing total error on cross section after matching efficiency
298 
299  double sigma2Sum, sigma2Err;
300  sigma2Sum = hepxsec_value * hepxsec_value;
301  sigma2Err = hepxsec_error * hepxsec_error;
302 
303  // sum of cross sections before matching and errors over different samples for each process
304  sum_denominator_before[ip] += (sigma2Err > 0)? 1/sigma2Err: 0;
305  sum_numerator_before[ip] += (sigma2Err > 0)? hepxsec_value/sigma2Err: 0;
306 
307 
308  double delta2Sum = delta2Veto
309  + sigma2Err / sigma2Sum;
310  relErr = (delta2Sum > 0.0 ?
311  std::sqrt(delta2Sum) : 0.0);
312  double deltaFin = sigmaFin * relErr;
313 
314  // sum of cross sections and errors over different processes
315  sigSelSum += hepxsec_value;
316  errSel2Sum += sigma2Err;
317  sigSum += sigmaFin;
318  err2Sum += deltaFin * deltaFin;
319 
320  // sum of cross sections after matching and errors over different samples for each process
321  sum_denominator_after[ip] += (deltaFin > 0)? 1/(deltaFin * deltaFin): 0;
322  sum_numerator_after[ip] += (deltaFin > 0)? sigmaFin/(deltaFin * deltaFin): 0;
323 
324 
325  } // end of loop over different processes
326 
327  sum_denominator_before[sizeOfInfos-1] += (errSel2Sum> 0)? 1/errSel2Sum: 0;
328  sum_numerator_before[sizeOfInfos-1] += (errSel2Sum> 0)? sigSelSum/errSel2Sum: 0;
329 
330  sum_denominator_after[sizeOfInfos-1] += (err2Sum>0)? 1/err2Sum: 0;
331  sum_numerator_after[sizeOfInfos-1] += (err2Sum>0)? sigSum/err2Sum: 0;
332 
333 
334  } // end of loop over different samples
335 
336 
337  for(unsigned int i=0; i<sizeOfInfos; i++)
338  {
339  double final_value = (sum_denominator_before[i]>0) ? (sum_numerator_before[i]/sum_denominator_before[i]):0;
340  double final_error = (sum_denominator_before[i]>0) ? (1/sqrt(sum_denominator_before[i])):-1;
341  xsecBeforeMatching_.push_back(GenLumiInfoProduct::XSec(final_value, final_error));
342 
343  double final_value2 = (sum_denominator_after[i]>0) ? (sum_numerator_after[i]/sum_denominator_after[i]):0;
344  double final_error2 = (sum_denominator_after[i]>0) ? 1/sqrt(sum_denominator_after[i]):-1;
345  xsecAfterMatching_.push_back(GenLumiInfoProduct::XSec(final_value2, final_error2));
346  }
347 
348  }
349  return;
350 }
351 
352 
353 void
355 
356  edm::LogPrint("GenXSecAnalyzer") << "\n"
357  << "------------------------------------" << "\n"
358  << "GenXsecAnalyzer:" << "\n"
359  << "------------------------------------";
360 
361  if(!products_.size()) {
362  edm::LogPrint("GenXSecAnalyzer") << "------------------------------------" << "\n"
363  << "Cross-section summary not available" << "\n"
364  << "------------------------------------";
365  return;
366  }
367 
368 
369  compute();
370 
371 
372  edm::LogPrint("GenXSecAnalyzer")
373  << "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \n"
374  << "Overall cross-section summary:" << "\n"
375  << "--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
376  edm::LogPrint("GenXSecAnalyzer")
377  << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw \txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
378 
379  const int sizeOfInfos= theProcesses_size+1;
380  const int last = sizeOfInfos-1;
381  std::string title[sizeOfInfos];
382 
383  for(int i=0; i < sizeOfInfos; i++){
384 
385  double jetmatch_eff=0;
386  double jetmatch_err=0;
387 
388  if(i==last)
389  {
390  title[i] = "Total";
391 
392  edm::LogPrint("GenXSecAnalyzer")
393  << "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";
394  jetmatch_eff = xsecBeforeMatching_[i].value()>0? xsecAfterMatching_[i].value()/xsecBeforeMatching_[i].value(): 0;
395  jetmatch_err = (xsecBeforeMatching_[i].value()>0 && xsecAfterMatching_[i].value()>0 &&
397  )? jetmatch_eff*sqrt(pow(xsecAfterMatching_[i].error()/xsecAfterMatching_[i].value(),2)-
399  }
400  else
401  {
402  std::ostringstream titless;
403  titless << i;
404  title[i] = titless.str();
405  jetmatch_eff = jetMatchEffStat_[i].filterEfficiency(hepidwtup_);
406  jetmatch_err = jetMatchEffStat_[i].filterEfficiencyError(hepidwtup_);
407  }
408 
409 
410  edm::LogPrint("GenXSecAnalyzer")
411  << title[i] << "\t\t"
412  << std::scientific << std::setprecision(3)
413  << xsecBeforeMatching_[i].value() << " +/- "
414  << xsecBeforeMatching_[i].error() << "\t\t"
415  << eventEffStat_[i].numEventsPassed() << "\t"
416  << jetMatchEffStat_[i].numPassPositiveEvents() << "\t"
417  << jetMatchEffStat_[i].numPassNegativeEvents() << "\t"
418  << eventEffStat_[i].numEventsTotal() << "\t"
419  << jetMatchEffStat_[i].numTotalPositiveEvents() << "\t"
420  << jetMatchEffStat_[i].numTotalNegativeEvents() << "\t"
421  << std::scientific << std::setprecision(3)
422  << xsecAfterMatching_[i].value() << " +/- "
423  << xsecAfterMatching_[i].error() << "\t\t"
424  << std::fixed << std::setprecision(1)
425  << (jetmatch_eff*100) << " +/- " << (jetmatch_err*100) << "\t"
426  << std::fixed << std::setprecision(1)
427  << (eventEffStat_[i].filterEfficiency(+3) * 100) << " +/- " << ( eventEffStat_[i].filterEfficiencyError(+3) * 100);
428 
429 
430  }
431 
432  edm::LogPrint("GenXSecAnalyzer")
433  << "--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
434 
435  edm::LogPrint("GenXSecAnalyzer")
436  << "Before matching: total cross section = "
437  << std::scientific << std::setprecision(3)
438  << xsecBeforeMatching_[last].value() << " +- " << xsecBeforeMatching_[last].error() << " pb";
439 
440  double xsec_match = xsecAfterMatching_[last].value();
441  double error_match = xsecAfterMatching_[last].error();
442 
443  edm::LogPrint("GenXSecAnalyzer")
444  << "After matching: total cross section = "
445  << std::scientific << std::setprecision(3)
446  << xsec_match << " +- " << error_match << " pb";
447 
448 
449  // hepMC filter efficiency
450  double hepMCFilter_eff = 1.0;
451  double hepMCFilter_err = 0.0;
452  if( hasHepMCFilterInfo_){
455 
456  edm::LogPrint("GenXSecAnalyzer")
457  << "HepMC filter efficiency (taking into account weights)= "
458  << "(" << hepMCFilterEffStat_.sumPassWeights() << ")"
459  << " / "
460  << "(" << hepMCFilterEffStat_.sumWeights() << ")"
461  << " = "
462  << std::scientific << std::setprecision(3)
463  << hepMCFilter_eff << " +- " << hepMCFilter_err;
464 
467  double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass/ hepMCFilter_event_total : 0;
468  double hepMCFilter_event_err = hepMCFilter_event_total > 0 ?
469  sqrt((1-hepMCFilter_event_eff)*hepMCFilter_event_eff/hepMCFilter_event_total): -1;
470  edm::LogPrint("GenXSecAnalyzer")
471  << "HepMC filter efficiency (event-level)= "
472  << "(" << hepMCFilter_event_pass << ")"
473  << " / "
474  << "(" << hepMCFilter_event_total << ")"
475  << " = "
476  << std::scientific << std::setprecision(3)
477  << hepMCFilter_event_eff << " +- " << hepMCFilter_event_err;
478  }
479 
480  // gen-particle filter efficiency
481  double filterOnly_eff = filterOnlyEffStat_.filterEfficiency(hepidwtup_);
482  double filterOnly_err = filterOnlyEffStat_.filterEfficiencyError(hepidwtup_);
483 
484  edm::LogPrint("GenXSecAnalyzer")
485  << "Filter efficiency (taking into account weights)= "
486  << "(" << filterOnlyEffStat_.sumPassWeights() << ")"
487  << " / "
488  << "(" << filterOnlyEffStat_.sumWeights() << ")"
489  << " = "
490  << std::scientific << std::setprecision(3)
491  << filterOnly_eff << " +- " << filterOnly_err;
492 
495  double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass/filterOnly_event_total : 0;
496  double filterOnly_event_err = filterOnly_event_total > 0 ?
497  sqrt((1-filterOnly_event_eff)*filterOnly_event_eff/filterOnly_event_total): -1;
498  edm::LogPrint("GenXSecAnalyzer")
499  << "Filter efficiency (event-level)= "
500  << "(" << filterOnly_event_pass << ")"
501  << " / "
502  << "(" << filterOnly_event_total << ")"
503  << " = "
504  << std::scientific << std::setprecision(3)
505  << filterOnly_event_eff << " +- " << filterOnly_event_err;
506 
507 
508  double xsec_final = xsec_match*filterOnly_eff*hepMCFilter_eff;
509  double error_final = xsec_final*sqrt(error_match*error_match/xsec_match/xsec_match+
510  filterOnly_err*filterOnly_err/filterOnly_eff/filterOnly_eff +
511  hepMCFilter_err*hepMCFilter_err/hepMCFilter_eff/hepMCFilter_eff
512  );
513 
514  edm::LogPrint("GenXSecAnalyzer")
515  << "After filter: final cross section = "
516  << std::scientific << std::setprecision(3)
517  << xsec_final << " +- " << error_final << " pb";
518 
519  xsec_ = GenLumiInfoProduct::XSec(xsec_final,error_final);
520 
521 
522 }
523 
std::vector< GenLumiInfoProduct::XSec > xsecBeforeMatching_
std::vector< GenFilterInfo > jetMatchEffStat_
int i
Definition: DBlmapReader.cc:9
std::vector< GenLumiInfoProduct > products_
unsigned int numTotalPositiveEvents() const
Definition: GenFilterInfo.h:28
TrainProcessor *const proc
Definition: MVATrainer.cc:101
double filterEfficiency(int idwtup=+3) const
GenFilterInfo filterOnlyEffStat_
bool getByLabel(std::string const &label, Handle< PROD > &result) const
unsigned int theProcesses_size
GenLumiInfoProduct::XSec xsec_
GenFilterInfo hepMCFilterEffStat_
bool mergeProduct(GenFilterInfo const &other)
std::vector< GenFilterInfo > eventEffStat_
double sumPassWeights() const
Definition: GenFilterInfo.h:34
T sqrt(T t)
Definition: SSEVec.h:46
virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
unsigned int numTotalNegativeEvents() const
Definition: GenFilterInfo.h:31
double filterEfficiencyError(int idwtup=+3) const
virtual void beginJob()
unsigned int numPassPositiveEvents() const
Definition: GenFilterInfo.h:27
virtual void endJob()
GenXSecAnalyzer(const edm::ParameterSet &)
std::vector< GenLumiInfoProduct::XSec > xsecAfterMatching_
virtual void analyze(const edm::Event &, const edm::EventSetup &)
double sumWeights() const
Definition: GenFilterInfo.h:40
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
unsigned int numPassNegativeEvents() const
Definition: GenFilterInfo.h:30