CMS 3D CMS Logo

CandOneToOneDeltaRMatcher.cc
Go to the documentation of this file.
1 /* \class CandOneToOneDeltaRMatcher
2  *
3  * Producer for simple match map
4  * to match two collections of candidate
5  * with one-to-One matching
6  * minimizing Sum(DeltaR)
7  *
8  */
9 
12 
14 
15 #include <vector>
16 #include <iostream>
17 
19 public:
21  ~CandOneToOneDeltaRMatcher() override;
22 
23 private:
24  void produce(edm::Event&, const edm::EventSetup&) override;
25  double lenght(const std::vector<int>&);
26  std::vector<int> AlgoBruteForce(int, int);
27  std::vector<int> AlgoSwitchMethod(int, int);
28 
31  std::vector<std::vector<float> > AllDist;
33 };
34 
36 
44 
49 
50 #include <Math/VectorUtil.h>
51 #include <TMath.h>
52 
53 using namespace edm;
54 using namespace std;
55 using namespace reco;
56 using namespace ROOT::Math::VectorUtil;
57 using namespace stdcomb;
58 
60  : sourceToken_(consumes<CandidateView>(cfg.getParameter<InputTag>("src"))),
61  matchedToken_(consumes<CandidateView>(cfg.getParameter<InputTag>("matched"))),
62  algoMethod_(cfg.getParameter<string>("algoMethod")) {
63  produces<CandViewMatchMap>("src2mtc");
64  produces<CandViewMatchMap>("mtc2src");
65 }
66 
68 
74 
75  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Source Collection =======";
76  for (CandidateView::const_iterator c = source->begin(); c != source->end(); ++c) {
77  edm::LogVerbatim("CandOneToOneDeltaRMatcher")
78  << " pt source " << c->pt() << " " << c->eta() << " " << c->phi() << endl;
79  }
80  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Matched Collection =======";
81  for (CandidateView::const_iterator c = matched->begin(); c != matched->end(); ++c) {
82  edm::LogVerbatim("CandOneToOneDeltaRMatcher")
83  << " pt source " << c->pt() << " " << c->eta() << " " << c->phi() << endl;
84  }
85 
86  const int nSrc = source->size();
87  const int nMtc = matched->size();
88 
89  const int nMin = min(source->size(), matched->size());
90  const int nMax = max(source->size(), matched->size());
91  if (nMin < 1)
92  return;
93 
94  if (nSrc <= nMtc) {
95  for (CandidateView::const_iterator iSr = source->begin(); iSr != source->end(); iSr++) {
96  vector<float> tempAllDist;
97  for (CandidateView::const_iterator iMt = matched->begin(); iMt != matched->end(); iMt++) {
98  tempAllDist.push_back(DeltaR(iSr->p4(), iMt->p4()));
99  }
100  AllDist.push_back(tempAllDist);
101  tempAllDist.clear();
102  }
103  } else {
104  for (CandidateView::const_iterator iMt = matched->begin(); iMt != matched->end(); iMt++) {
105  vector<float> tempAllDist;
106  for (CandidateView::const_iterator iSr = source->begin(); iSr != source->end(); iSr++) {
107  tempAllDist.push_back(DeltaR(iSr->p4(), iMt->p4()));
108  }
109  AllDist.push_back(tempAllDist);
110  tempAllDist.clear();
111  }
112  }
113 
114  /*
115  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== The DeltaR Matrix =======";
116  for(int m0=0; m0<nMin; m0++) {
117  // for(int m1=0; m1<nMax; m1++) {
118  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << setprecision(2) << fixed << (m1 AllDist[m0][m1] ;
119  //}
120  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "\n";
121  }
122  */
123 
124  // Loop size if Brute Force
125  int nLoopToDo = (int)(TMath::Factorial(nMax) / TMath::Factorial(nMax - nMin));
126  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "nLoop:" << nLoopToDo << endl;
127  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "Choosen Algo is:" << algoMethod_;
128  vector<int> bestCB;
129 
130  // Algo is Brute Force
131  if (algoMethod_ == "BruteForce") {
132  bestCB = AlgoBruteForce(nMin, nMax);
133 
134  // Algo is Switch Method
135  } else if (algoMethod_ == "SwitchMode") {
136  bestCB = AlgoSwitchMethod(nMin, nMax);
137 
138  // Algo is Brute Force if nLoop < 10000
139  } else if (algoMethod_ == "MixMode") {
140  if (nLoopToDo < 10000) {
141  bestCB = AlgoBruteForce(nMin, nMax);
142  } else {
143  bestCB = AlgoSwitchMethod(nMin, nMax);
144  }
145 
146  } else {
147  throw cms::Exception("OneToOne Constructor") << "wrong matching method in ParameterSet";
148  }
149 
150  for (int i1 = 0; i1 < nMin; i1++)
151  edm::LogVerbatim("CandOneToOneDeltaRMatcher")
152  << "min: " << i1 << " " << bestCB[i1] << " " << AllDist[i1][bestCB[i1]];
153 
154  /*
155  auto matchMapSrMt = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( source ),
156  CandidateRefProd( matched ) ) );
157  auto matchMapMtSr = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( matched ),
158  CandidateRefProd( source ) ) );
159 */
160 
161  auto matchMapSrMt = std::make_unique<CandViewMatchMap>();
162  auto matchMapMtSr = std::make_unique<CandViewMatchMap>();
163 
164  for (int c = 0; c != nMin; c++) {
165  if (source->size() <= matched->size()) {
166  matchMapSrMt->insert(source->refAt(c), matched->refAt(bestCB[c]));
167  matchMapMtSr->insert(matched->refAt(bestCB[c]), source->refAt(c));
168  } else {
169  matchMapSrMt->insert(source->refAt(bestCB[c]), matched->refAt(c));
170  matchMapMtSr->insert(matched->refAt(c), source->refAt(bestCB[c]));
171  }
172  }
173 
174  /*
175  for( int c = 0; c != nMin; c ++ ) {
176  if( source->size() <= matched->size() ) {
177  matchMapSrMt->insert( CandidateRef( source, c ), CandidateRef( matched, bestCB[c] ) );
178  matchMapMtSr->insert( CandidateRef( matched, bestCB[c] ), CandidateRef( source, c ) );
179  } else {
180  matchMapSrMt->insert( CandidateRef( source, bestCB[c] ), CandidateRef( matched, c ) );
181  matchMapMtSr->insert( CandidateRef( matched, c ), CandidateRef( source, bestCB[c] ) );
182  }
183  }
184 */
185  evt.put(std::move(matchMapSrMt), "src2mtc");
186  evt.put(std::move(matchMapMtSr), "mtc2src");
187 
188  AllDist.clear();
189 }
190 
191 double CandOneToOneDeltaRMatcher::lenght(const vector<int>& best) {
192  double myLenght = 0;
193  int row = 0;
194  for (vector<int>::const_iterator it = best.begin(); it != best.end(); it++) {
195  myLenght += AllDist[row][*it];
196  row++;
197  }
198  return myLenght;
199 }
200 
201 // this is the Brute Force Algorithm
202 // All the possible combination are checked
203 // The best one is always found
204 // Be carefull when you have high values for nMin and nMax --> the combinatorial could explode!
205 // Sum(DeltaR) is minimized -->
206 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
207 // 0.1 - 0.2 - 0.3 - 3.0
208 // Which one do you prefer? --> BruteForce select always the first
209 
211  vector<int> ca;
212  vector<int> cb;
213  vector<int> bestCB;
214  float totalDeltaR = 0;
215  float BestTotalDeltaR = 1000;
216 
217  ca.reserve(nMax);
218  for (int i1 = 0; i1 < nMax; i1++)
219  ca.push_back(i1);
220  cb.reserve(nMin);
221  for (int i1 = 0; i1 < nMin; i1++)
222  cb.push_back(i1);
223 
224  do {
225  //do your processing on the new combination here
226  for (int cnt = 0; cnt < TMath::Factorial(nMin); cnt++) {
227  totalDeltaR = lenght(cb);
228  if (totalDeltaR < BestTotalDeltaR) {
229  BestTotalDeltaR = totalDeltaR;
230  bestCB = cb;
231  }
232  next_permutation(cb.begin(), cb.end());
233  }
234  } while (next_combination(ca.begin(), ca.end(), cb.begin(), cb.end()));
235 
236  return bestCB;
237 }
238 
239 // This method (Developed originally by Daniele Benedetti) check for the best combination
240 // choosing the minimum DeltaR for each line in AllDist matrix
241 // If no repeated row is found: ie (line,col)=(1,3) and (2,3) --> same as BruteForce
242 // If repetition --> set the higher DeltaR between the 2 repetition to 1000 and re-check best combination
243 // Iterate until no repetition
244 // No guaranted minimum for Sum(DeltaR)
245 // If you have:
246 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
247 // 0.1 - 0.2 - 0.3 - 3.0
248 // SwitchMethod normally select the second solution
249 
251  vector<int> bestCB;
252  for (int i1 = 0; i1 < nMin; i1++) {
253  int minInd = 0;
254  for (int i2 = 1; i2 < nMax; i2++)
255  if (AllDist[i1][i2] < AllDist[i1][minInd])
256  minInd = i2;
257  bestCB.push_back(minInd);
258  }
259 
260  bool inside = true;
261  while (inside) {
262  inside = false;
263  for (int i1 = 0; i1 < nMin; i1++) {
264  for (int i2 = i1 + 1; i2 < nMin; i2++) {
265  if (bestCB[i1] == bestCB[i2]) {
266  inside = true;
267  if (AllDist[i1][(bestCB[i1])] <= AllDist[i2][(bestCB[i2])]) {
268  AllDist[i2][(bestCB[i2])] = 1000;
269  int minInd = 0;
270  for (int i3 = 1; i3 < nMax; i3++)
271  if (AllDist[i2][i3] < AllDist[i2][minInd])
272  minInd = i3;
273  bestCB[i2] = minInd;
274  } else {
275  AllDist[i1][(bestCB[i1])] = 1000;
276  int minInd = 0;
277  for (int i3 = 1; i3 < nMax; i3++)
278  if (AllDist[i1][i3] < AllDist[i1][minInd])
279  minInd = i3;
280  bestCB[i1] = minInd;
281  }
282  } // End if
283  }
284  }
285  } // End while
286 
287  return bestCB;
288 }
289 
292 
muonTagProbeFilters_cff.matched
matched
Definition: muonTagProbeFilters_cff.py:62
CandOneToOneDeltaRMatcher::sourceToken_
edm::EDGetTokenT< reco::CandidateView > sourceToken_
Definition: CandOneToOneDeltaRMatcher.cc:29
testProducerWithPsetDescEmpty_cfi.i2
i2
Definition: testProducerWithPsetDescEmpty_cfi.py:46
Handle.h
CandOneToOneDeltaRMatcher::~CandOneToOneDeltaRMatcher
~CandOneToOneDeltaRMatcher() override
Definition: CandOneToOneDeltaRMatcher.cc:67
testProducerWithPsetDescEmpty_cfi.i3
i3
Definition: testProducerWithPsetDescEmpty_cfi.py:47
MessageLogger.h
EDProducer.h
ESHandle.h
CandOneToOneDeltaRMatcher::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: CandOneToOneDeltaRMatcher.cc:69
metDiagnosticParameterSet_cfi.nMin
nMin
Definition: metDiagnosticParameterSet_cfi.py:10
CandOneToOneDeltaRMatcher::AlgoSwitchMethod
std::vector< int > AlgoSwitchMethod(int, int)
Definition: CandOneToOneDeltaRMatcher.cc:250
min
T min(T a, T b)
Definition: MathUtil.h:58
edm::EDGetTokenT
Definition: EDGetToken.h:33
edm
HLT enums.
Definition: AlignableModifier.h:19
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
CandOneToOneDeltaRMatcher::CandOneToOneDeltaRMatcher
CandOneToOneDeltaRMatcher(const edm::ParameterSet &)
Definition: CandOneToOneDeltaRMatcher.cc:59
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
CandOneToOneDeltaRMatcher::AllDist
std::vector< std::vector< float > > AllDist
Definition: CandOneToOneDeltaRMatcher.cc:31
edm::Handle
Definition: AssociativeIterator.h:50
CandMatchMap.h
CandidateFwd.h
EDMException.h
MakerMacros.h
stdcomb
Definition: combination.h:14
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
CandOneToOneDeltaRMatcher::matchedToken_
edm::EDGetTokenT< reco::CandidateView > matchedToken_
Definition: CandOneToOneDeltaRMatcher.cc:30
source
static const std::string source
Definition: EdmProvDump.cc:47
LeafCandidate.h
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:531
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::View
Definition: CaloClusterFwd.h:14
edm::ParameterSet
Definition: ParameterSet.h:47
Event.h
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
ModuleDef.h
createfilelist.int
int
Definition: createfilelist.py:10
CandOneToOneDeltaRMatcher::algoMethod_
std::string algoMethod_
Definition: CandOneToOneDeltaRMatcher.cc:32
edm::Event::put
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
CandOneToOneDeltaRMatcher
Definition: CandOneToOneDeltaRMatcher.cc:18
edm::EventSetup
Definition: EventSetup.h:57
electronAnalyzer_cfi.DeltaR
DeltaR
Definition: electronAnalyzer_cfi.py:33
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
InputTag.h
looper.cfg
cfg
Definition: looper.py:297
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
stdcomb::next_combination
bool next_combination(BidIt n_begin, BidIt n_end, BidIt r_begin, BidIt r_end)
Definition: combination.h:19
CandOneToOneDeltaRMatcher::lenght
double lenght(const std::vector< int > &)
Definition: CandOneToOneDeltaRMatcher.cc:191
CandOneToOneDeltaRMatcher::AlgoBruteForce
std::vector< int > AlgoBruteForce(int, int)
Definition: CandOneToOneDeltaRMatcher.cc:210
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
Exception
Definition: hltDiff.cc:246
ParameterSetfwd.h
EventSetup.h
edm::EDProducer
Definition: EDProducer.h:35
edm::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
Candidate.h
ParameterSet.h
edm::Event
Definition: Event.h:73
edm::Log
Definition: MessageLogger.h:70
edm::InputTag
Definition: InputTag.h:15
metDiagnosticParameterSet_cfi.nMax
nMax
Definition: metDiagnosticParameterSet_cfi.py:11
combination.h