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