VHbbDataFormats/interface/MultiThresholdEfficiency.h

#ifndef MULTITHRESHOLDEFFICIENCY_H
#define MULTITHRESHOLDEFFICIENCY_H
#include <math.h>
#include <iostream>
#include <vector>
using namespace std; 
class MultiThresholdEfficiency 
{
 public:
   MultiThresholdEfficiency(unsigned int nOfThresholds) : thresholds(nOfThresholds) {}
   
   //A vector with, for each object the efficiency of passing threshold N while notpassing N-1 
   template <class Filter> float weight(vector<vector<float> > input);

 private:
  unsigned int thresholds;

};


// NB here sorting is in opposite direction than for btag 0 is tigther , 1 is looser
class Trigger1High2Loose
{
 public:
  static bool filter(std::vector<int> t)
  {
    return t[0] >= 1 && t[1] >= 2;
}
};

class Trigger2SingleThr
{
 public:
  static bool filter(std::vector<int> t)
  {
    return t[0] >= 2;
}
};


template <class Filter> float MultiThresholdEfficiency::weight(vector<vector<float> > objects)
{
 unsigned int nobjects=objects.size();
 std::vector<unsigned int> comb(objects.size());
 for(unsigned int i=0;i < objects.size(); i++) comb[i]=0;
 unsigned int idx=0;
 unsigned int max=thresholds+1; // 
 float p=0,tot=0;
 if(objects.size()==0) return 0.;
 while(comb[objects.size()-1] < max)
 {
//std::cout << std::endl << "New comb" << std::endl;
// for(int i=0;i < objects.size(); i++) {std::cout << comb[i] << " ";}
// std::cout << std::endl;
   std::vector<int> pass;
   for(unsigned int j=0;j<thresholds;j++) pass.push_back(0);

   float thisCombinationProbability=1.;
//   std::cout << "OBJ Probs: ";
   for(size_t j=0;j<nobjects;j++) // loop on objects
    {
       float cumulative=1.;
       for(unsigned int n=0;n<comb[j];n++) cumulative*= (1.-objects[j][n]); // 10 20 70  10/100,20/90   90/100*70/90 
       thisCombinationProbability*=cumulative;
       if(comb[j]< thresholds) {
         thisCombinationProbability*=(objects[j])[comb[j]];     
//         std::cout << cumulative*(objects[j])[comb[j]] << " "; 
       }
/*       else
       { 
              std::cout << cumulative << " " ;
       }*/

 
     for(size_t k=0;k< thresholds; k++ ) // loop on threshold
      {
       bool passed = ( k >= comb[j] ) ;  //k=0, is the tightest, passed only if comb[j] = 0, k=1 pass if comb 0 or 1
       if(passed) pass[k]++;
      }
  }
//   std::cout << endl;
   if(Filter::filter(pass)) p+=thisCombinationProbability;
   tot+=thisCombinationProbability;
//   std::cout << thisCombinationProbability << " " << p << " " << tot<< std::endl;
  while (comb[idx] == max -1  && idx+1 < objects.size()) idx++; // find first object for which we did not already test all configs 
 // next combination
  comb[idx]++;  // test a new config for that jet
  for(unsigned int i=0;i<idx;i++) { comb[i]=0; } // reset the tested configs for all previous objects
  idx=0;
 }
  if( fabs(tot-1.) > 0.01 )
   {
     std::cout << "ERROR, total must be one " << tot << std::endl;
   }
  return p;
}

#endif

Revision:	1.3
Committed:	Thu Nov 10 00:01:40 2011 UTC (13 years, 5 months ago) by arizzi
Content type:	text/plain
Branch:	MAIN
CVS Tags:	EDMV42_Step2_V8, EDMV42_Step2_V7, EDMV42_Step2_V6, EDMV42_Step2_V5a, EDMV42_Step2_V5, tauCandV42, hbbsubstructDev_11, hbbsubstructDev_10, hbbsubstructDev_9, hbbsubstructDev_8, hbbsubstructDev_7, hbbsubstructDev_6, hbbsubstructDev_5, hbbsubstructDev_4, hbbsubstructDev_3, hbbsubstructDev_2, hbbsubstructDev_1, hbbsubstructDev, V21TauCand_0, EDMV42_Step2_V4a, EDMV42_Step2_V4, EDMV42_Step2_V3, EDMV42_Step2_V2, EDMV42_Step2_V1, EdmV42, EdmV41alpha1, EdmV40alpha1, EdmV40alpha, V21emuCand, EdmV33Jun12v2_consistent, Step2ForV33_v2, Step2ForV33_v1, EdmV33Jun12v2, EdmV33Jun12v1, EdmV33Jun12v0, Step2ForV32_v2, Step2ForV32_v1, Step2ForV32_v0, Step2ForV31_v0, EdmV32May24v0, EdmV31May21v1, EdmV31May17v0, May14thStep2, EdmV30Apr10, EdmV21Apr10v2, EdmV22May9, EdmV21Apr06, EdmV21Apr10, EdmV21Apr04, EdmV21Apr03, EdmV21Apr2, EdmV21Mar30, EdmV20Mar12, AR_Nov10Ntuple, HEAD
Branch point for:	V42TauCandidate, hbbsubstructDevPostHCP, V21TauCand, V21emuCandidate
Changes since 1.2:	+9 -0 lines
Log Message:	- 3D PU reweighting with latest recipe and targetting the full 2011B run period - PU reweight fixed in 2011A against Fall11 MC - run boundary handling to properly handle the two periods of 2011B - factorized MET triggers - events failing the "at least 2 jets pt > 20" are recovered with a pt15 cut (NB: the election of higest dijet is first done with 20GeV threshold, the 15 GeV is consiered only if the 20 fails) - apply JER smearing by default - add leptons generator info - add jet by jet btag SF info - add electron id80NoIso
#	User	Rev	Content
1	arizzi	1.1	#ifndef MULTITHRESHOLDEFFICIENCY_H
2			#define MULTITHRESHOLDEFFICIENCY_H
3			#include <math.h>
4			#include <iostream>
5			#include <vector>
6			using namespace std;
7			class MultiThresholdEfficiency
8			{
9			public:
10	arizzi	1.2	MultiThresholdEfficiency(unsigned int nOfThresholds) : thresholds(nOfThresholds) {}
11	arizzi	1.1
12			//A vector with, for each object the efficiency of passing threshold N while notpassing N-1
13	arizzi	1.2	template <class Filter> float weight(vector<vector<float> > input);
14	arizzi	1.1
15			private:
16	arizzi	1.2	unsigned int thresholds;
17	arizzi	1.1
18			};
19
20
21	arizzi	1.2	// NB here sorting is in opposite direction than for btag 0 is tigther , 1 is looser
22			class Trigger1High2Loose
23	arizzi	1.1	{
24	arizzi	1.2	public:
25			static bool filter(std::vector<int> t)
26			{
27			return t[0] >= 1 && t[1] >= 2;
28	arizzi	1.1	}
29	arizzi	1.2	};
30
31	arizzi	1.3	class Trigger2SingleThr
32			{
33			public:
34			static bool filter(std::vector<int> t)
35			{
36			return t[0] >= 2;
37			}
38			};
39
40	arizzi	1.1
41
42
43	arizzi	1.2	template <class Filter> float MultiThresholdEfficiency::weight(vector<vector<float> > objects)
44	arizzi	1.1	{
45	arizzi	1.2	unsigned int nobjects=objects.size();
46			std::vector<unsigned int> comb(objects.size());
47			for(unsigned int i=0;i < objects.size(); i++) comb[i]=0;
48			unsigned int idx=0;
49			unsigned int max=thresholds+1; //
50	arizzi	1.1	float p=0,tot=0;
51			if(objects.size()==0) return 0.;
52			while(comb[objects.size()-1] < max)
53			{
54			//std::cout << std::endl << "New comb" << std::endl;
55			// for(int i=0;i < objects.size(); i++) {std::cout << comb[i] << " ";}
56			// std::cout << std::endl;
57			std::vector<int> pass;
58	arizzi	1.2	for(unsigned int j=0;j<thresholds;j++) pass.push_back(0);
59	arizzi	1.1
60			float thisCombinationProbability=1.;
61			// std::cout << "OBJ Probs: ";
62			for(size_t j=0;j<nobjects;j++) // loop on objects
63			{
64			float cumulative=1.;
65	arizzi	1.2	for(unsigned int n=0;n<comb[j];n++) cumulative= (1.-objects[j][n]); // 10 20 70 10/100,20/90 90/10070/90
66	arizzi	1.1	thisCombinationProbability*=cumulative;
67			if(comb[j]< thresholds) {
68			thisCombinationProbability*=(objects[j])[comb[j]];
69			// std::cout << cumulative*(objects[j])[comb[j]] << " ";
70			}
71			/* else
72			{
73			std::cout << cumulative << " " ;
74			}*/
75
76
77			for(size_t k=0;k< thresholds; k++ ) // loop on threshold
78			{
79			bool passed = ( k >= comb[j] ) ; //k=0, is the tightest, passed only if comb[j] = 0, k=1 pass if comb 0 or 1
80			if(passed) pass[k]++;
81			}
82			}
83			// std::cout << endl;
84	arizzi	1.2	if(Filter::filter(pass)) p+=thisCombinationProbability;
85	arizzi	1.1	tot+=thisCombinationProbability;
86			// std::cout << thisCombinationProbability << " " << p << " " << tot<< std::endl;
87			while (comb[idx] == max -1 && idx+1 < objects.size()) idx++; // find first object for which we did not already test all configs
88			// next combination
89			comb[idx]++; // test a new config for that jet
90	arizzi	1.2	for(unsigned int i=0;i<idx;i++) { comb[i]=0; } // reset the tested configs for all previous objects
91	arizzi	1.1	idx=0;
92			}
93			if( fabs(tot-1.) > 0.01 )
94			{
95			std::cout << "ERROR, total must be one " << tot << std::endl;
96			}
97			return p;
98			}
99
100			#endif
101