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;
}
};


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.2
Committed:	Wed Sep 14 14:31:20 2011 UTC (13 years, 7 months ago) by arizzi
Content type:	text/plain
Branch:	MAIN
CVS Tags:	AR_step2_Oct25, AR_step2_oct19, EdmV11Oct2011, AR_Step2_Oct13, AR_Oct9Ntuple, AR_Oct7_step2ntuple, AR_Oct5Ntuple, EdmV10Oct2011, EdmV9Sept2011, Sept19th2011_2, Sept19th2011, Sept19th, VHNtupleV9_AR1, VHSept15_AR1, Sept14th2011_AR1
Changes since 1.1:	+21 -19 lines
Log Message:	fixes
#	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.1
32
33
34	arizzi	1.2	template <class Filter> float MultiThresholdEfficiency::weight(vector<vector<float> > objects)
35	arizzi	1.1	{
36	arizzi	1.2	unsigned int nobjects=objects.size();
37			std::vector<unsigned int> comb(objects.size());
38			for(unsigned int i=0;i < objects.size(); i++) comb[i]=0;
39			unsigned int idx=0;
40			unsigned int max=thresholds+1; //
41	arizzi	1.1	float p=0,tot=0;
42			if(objects.size()==0) return 0.;
43			while(comb[objects.size()-1] < max)
44			{
45			//std::cout << std::endl << "New comb" << std::endl;
46			// for(int i=0;i < objects.size(); i++) {std::cout << comb[i] << " ";}
47			// std::cout << std::endl;
48			std::vector<int> pass;
49	arizzi	1.2	for(unsigned int j=0;j<thresholds;j++) pass.push_back(0);
50	arizzi	1.1
51			float thisCombinationProbability=1.;
52			// std::cout << "OBJ Probs: ";
53			for(size_t j=0;j<nobjects;j++) // loop on objects
54			{
55			float cumulative=1.;
56	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
57	arizzi	1.1	thisCombinationProbability*=cumulative;
58			if(comb[j]< thresholds) {
59			thisCombinationProbability*=(objects[j])[comb[j]];
60			// std::cout << cumulative*(objects[j])[comb[j]] << " ";
61			}
62			/* else
63			{
64			std::cout << cumulative << " " ;
65			}*/
66
67
68			for(size_t k=0;k< thresholds; k++ ) // loop on threshold
69			{
70			bool passed = ( k >= comb[j] ) ; //k=0, is the tightest, passed only if comb[j] = 0, k=1 pass if comb 0 or 1
71			if(passed) pass[k]++;
72			}
73			}
74			// std::cout << endl;
75	arizzi	1.2	if(Filter::filter(pass)) p+=thisCombinationProbability;
76	arizzi	1.1	tot+=thisCombinationProbability;
77			// std::cout << thisCombinationProbability << " " << p << " " << tot<< std::endl;
78			while (comb[idx] == max -1 && idx+1 < objects.size()) idx++; // find first object for which we did not already test all configs
79			// next combination
80			comb[idx]++; // test a new config for that jet
81	arizzi	1.2	for(unsigned int i=0;i<idx;i++) { comb[i]=0; } // reset the tested configs for all previous objects
82	arizzi	1.1	idx=0;
83			}
84			if( fabs(tot-1.) > 0.01 )
85			{
86			std::cout << "ERROR, total must be one " << tot << std::endl;
87			}
88			return p;
89			}
90
91			#endif
92