ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/L1RpcTriggerAnalysis/src/GoldenPattern.cc

Comparing UserCode/L1RpcTriggerAnalysis/src/GoldenPattern.cc (file contents):
Revision 1.5 by akalinow, Tue Jun 25 10:49:35 2013 UTC vs.
Revision 1.6 by akalinow, Tue Jun 25 11:06:02 2013 UTC

#	Line 1 | Line 1
1	+	#include <functional>
2	+	#include <numeric>
3	+
4		#include "UserCode/L1RpcTriggerAnalysis/interface/GoldenPattern.h"
5
6		#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#	Line 10 | Line 13
13		#include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
14		#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"
15
16	<
14	<	void GoldenPattern::Result::runNoCheck() const
15	<	{
16	>	void GoldenPattern::Result::runNoCheck() const {
17		double fract = 1;
18	<	for (auto i=posRpcResult.begin(); i!=posRpcResult.end();i++) fract *= norm(POSRPC,i->second);
19	<	for (auto i=posCscResult.begin(); i!=posCscResult.end();i++) fract *= norm(POSCSC,i->second);
20	<	for (auto i=posDtResult.begin();  i!=posDtResult.end();i++)  fract *= norm(POSDT, i->second);
21	<	for (auto i=benCscResult.begin(); i!=benCscResult.end();i++) fract *= norm(BENCSC,i->second);
22	<	for (auto i=benDtResult.begin();  i!=benDtResult.end();i++)  fract *= norm(BENDT, i->second);
23	<	unsigned int nTot = nMatchedPosRpc+nMatchedPosCsc+nMatchedPosDt+nMatchedBenCsc+nMatchedBenDt;
18	>
19	>	for(auto mType=myResults.cbegin();mType!=myResults.cend();++mType){
20	>	for (auto it=mType->second.cbegin(); it!=mType->second.cend();++it) fract *= norm(mType->first,it->second);
21	>	}
22	>	unsigned int nTot = 0;
23	>	for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;
24	>
25	>	theValue*=(myResults[GoldenPattern::POSRPC].size()==nMatchedPoints[GoldenPattern::POSRPC]);
26	>	theValue*=(myResults[GoldenPattern::POSDT].size()==nMatchedPoints[GoldenPattern::POSDT]);
27	>	theValue*=(myResults[GoldenPattern::POSCSC].size()==nMatchedPoints[GoldenPattern::POSCSC]);
28		theValue = ( nTot > 4) ? pow(fract, 1./((double) nTot)) : 0.;
24	–	if (posRpcResult.size() != nMatchedPosRpc) theValue = 0.;
25	–	if (posCscResult.size() != nMatchedPosCsc) theValue = 0.;
26	–	if (posDtResult.size()  != nMatchedPosDt)  theValue = 0.;
29		}
30
31		double GoldenPattern::Result::norm(GoldenPattern::PosBenCase where, double whereInDist) const {
32		double normValue = 2.*(0.5-fabs(whereInDist-0.5));
33		static const double epsilon = 1.e-9;
34	<	if (normValue > epsilon) {
35	<	switch (where) {
34	<	case POSRPC : nMatchedPosRpc++; break;
35	<	case POSCSC : nMatchedPosCsc++; break;
36	<	case BENCSC : nMatchedBenCsc++; break;
37	<	case POSDT  : nMatchedPosDt++;  break;
38	<	case BENDT  : nMatchedBenDt++;  break;
39	<	};
40	<	} else {
41	<	normValue =1.;
42	<	}
34	>	if (normValue > epsilon) ++nMatchedPoints[where];
35	>	else normValue =1.;
36		return normValue;
37		}
38
#	Line 68 | Line 61 | double GoldenPattern::Result::value() co
61
62		unsigned int GoldenPattern::Result::nMatchedTot() const {
63		run();
64	<	return nMatchedPosRpc+nMatchedPosCsc+nMatchedBenCsc+nMatchedPosDt+nMatchedBenDt;
64	>	unsigned int nTot = 0;
65	>	for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;
66	>	return nTot;
67		}
68
69
70		std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
71		{
72		o.run();
73	+
74	+	/*
75		out <<"Result: "
76		<< " value: "<<o.theValue
77	<	<<" nPos+Ben: ("<<o.nMatchedPosCsc<<"/"<<o.posCscResult.size()<<"+"<<o.nMatchedBenCsc<<"/"<<o.benCscResult.size()
78	<	<<", "<<o.nMatchedPosDt <<"/"<<o.posDtResult.size()<<"+"<<o.nMatchedBenDt<<"/"<<o.benDtResult.size()
79	<	<<", "<<o.nMatchedPosRpc<<"/"<<o.posRpcResult.size()<<", tot:"<<o.nMatchedTot()<<")";
77	>	<<" nPos+Ben: (";
78	>	for(auto mType=o.myResults.cbegin();
79	>	mType!=myResults.cend();++mType){
80	>	out<<o.nMatchedPoints[iType]<<"/"<<o.myResults[mType].size()<<", ";
81	>	}
82	>	out <<", tot:"<<o.nMatchedTot()<<")";
83	>	*/
84	>	/*
85	>	<<o.nMatchedPosCsc<<"/"<<o.posCscResult.size()<<"+"<<o.nMatchedBenCsc<<"/"<<o.benCscResult.size()
86	>	<<", "<<o.nMatchedPosDt <<"/"<<o.posDtResult.size()<<"+"<<o.nMatchedBenDt<<"/"<<o.benDtResult.size()
87	>	<<", "<<o.nMatchedPosRpc<<"/"<<o.posRpcResult.size()
88	>	*/
89		return out;
90		}
91
92	<
93	<
88	<	void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq)
89	<	{
90	<	switch ( aCase ) {
91	<	case POSRPC : posRpc[rawId][posOrBen]     += freq; break;
92	<	case POSCSC : posCsc[rawId][posOrBen]     += freq; break;
93	<	case BENCSC : bendingCsc[rawId][posOrBen] += freq; break;
94	<	case POSDT  : posDt[rawId][posOrBen]      += freq; break;
95	<	case BENDT  : bendingDt[rawId][posOrBen]  += freq; break;
96	<	};
92	>	void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq){
93	>	PattCore[aCase][rawId][posOrBen] += freq;
94		}
95
96
97	<	void GoldenPattern::add(const Pattern & p)
98	<	{
97	>	void GoldenPattern::add(const Pattern & p) {
98	>
99		const Pattern::DataType & detdigi = p ;
100	<	for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
100	>	for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
101		uint32_t rawId = is->first;
102		DetId detId(rawId);
103	<	if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
103	>	if (detId.det() != DetId::Muon)
104	>	std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
105		switch (detId.subdetId()) {
106		case MuonSubdetId::RPC: {
107		RPCDigiSpec digi(rawId, is->second);
108	<	posRpc[rawId][digi.halfStrip()]++;
108	>	PattCore[GoldenPattern::POSRPC][rawId][digi.halfStrip()]++;
109		break;
110		}
111		case MuonSubdetId::DT: {
112		DTphDigiSpec digi(rawId, is->second);
113	<	if (digi.bxNum() != 0 || digi.bxCnt() != 0 || digi.ts2() != 0) break;
114	<	posDt[rawId][digi.phi()]++;
115	<	bendingDt[rawId][digi.phiB()]++;
113	>	if (digi.bxNum() != 0 || digi.bxCnt() != 0 || digi.ts2() != 0) break;
114	>	PattCore[GoldenPattern::POSDT][rawId][digi.phi()]++;
115	>	PattCore[GoldenPattern::BENDT][rawId][digi.phiB()]++;
116		break;
117		}
118		case MuonSubdetId::CSC: {
119		CSCDigiSpec digi(rawId, is->second);
120	<	posCsc[rawId][digi.strip()]++;
121	<	bendingCsc[rawId][digi.pattern()]++;
120	>	PattCore[GoldenPattern::POSCSC][rawId][digi.strip()]++;
121	>	PattCore[GoldenPattern::BENCSC][rawId][digi.pattern()]++;
122		break;
123		}
124		default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
#	Line 128 | Line 126 | void GoldenPattern::add(const Pattern &
126		}
127		}
128
129	<	GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const
129	>	GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const
130		{
131		Result result;
134	–
132		const Pattern::DataType & detdigi = p;
133	<	for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
133	>	uint32_t  nTot = 0;
134	>	///Check spatial compatibility of GP with Pattern.
135	>	///Require at least 5 measurements in common detUnits
136	>	for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
137	>	uint32_t rawId = is->first;
138	>	for(auto mType = PattCore.cbegin();mType!=PattCore.cend();++mType){
139	>	nTot+=mType->second.count(rawId);
140	>	}
141	>	}
142	>	if(nTot<5) return result;
143	>
144	>
145	>	SystFreq::const_iterator cit;
146	>	DetFreq::const_iterator idm;
147	>	PosBenCase mType;
148	>	for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
149		uint32_t rawId = is->first;
150		DetId detId(rawId);
151	<	if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
151	>	if (detId.det() != DetId::Muon){
152	>	std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
153	>	return result;
154	>	}
155		if (detId.subdetId() == MuonSubdetId::RPC) {
156		RPCDigiSpec digi(rawId, is->second);
157	<	DetFreq::const_iterator idm = posRpc.find(rawId);
158	<	if (idm != posRpc.end() ) {
157	>	mType = GoldenPattern::POSRPC;
158	>	cit = PattCore.find(mType);
159	>	if(cit==PattCore.cend()){
160	>	std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
161	>	return result;
162	>	}
163	>	idm = cit->second.find(rawId);
164	>	if (idm != cit->second.cend() ) {
165		double f = whereInDistribution(digi.halfStrip(), idm->second);
166	<	result.posRpcResult.push_back( std::make_pair(rawId, f) );
166	>	result.myResults[mType].push_back(std::make_pair(rawId, f));
167		RPCDetId rpc(rawId);
168		if(rpc.station()==1) result.hasStation1 = true;
169		if(rpc.station()==2) result.hasStation2 = true;
170		}
171	<	} else if (detId.subdetId() == MuonSubdetId::DT) {
171	>	}
172	>	else if (detId.subdetId() == MuonSubdetId::DT) {
173		DTphDigiSpec digi(rawId, is->second);
174	<	DetFreq::const_iterator idm = posDt.find(rawId);
175	<	if (idm != posDt.end() ) {
174	>	mType = GoldenPattern::POSDT;
175	>	cit = PattCore.find(mType);
176	>	if(cit==PattCore.cend()){
177	>	std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
178	>	return result;
179	>	}
180	>	idm = cit->second.find(rawId);
181	>	if (idm != cit->second.cend() ) {
182		double f = whereInDistribution(digi.phi(), idm->second);
183	<	result.posDtResult.push_back( std::make_pair(rawId, f) );
183	>	result.myResults[mType].push_back(std::make_pair(rawId, f));
184		DTChamberId dt(rawId);
185		if(dt.station()==1) result.hasStation1 = true;
186		if(dt.station()==2) result.hasStation2 = true;
187		}
188	<	idm = bendingDt.find(rawId);
189	<	if (idm != bendingDt.end() ) {
190	<	double f = whereInDistribution(digi.phiB(), idm->second);
191	<	result.benDtResult.push_back( std::make_pair(rawId, f) );
188	>	mType = GoldenPattern::BENDT;
189	>	cit = PattCore.find(mType);
190	>	if(cit==PattCore.cend()){
191	>	std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
192	>	return result;
193	>	}
194	>	idm = cit->second.find(rawId);
195	>	if (idm != cit->second.cend() ) {
196	>	double f = whereInDistribution(digi.phiB(), idm->second);
197	>	result.myResults[mType].push_back(std::make_pair(rawId, f));
198		}
199	<	} else if (detId.subdetId() == MuonSubdetId::CSC) {
199	>	}
200	>	else if (detId.subdetId() == MuonSubdetId::CSC) {
201		CSCDigiSpec digi(rawId, is->second);
202	<	DetFreq::const_iterator idm = posCsc.find(rawId);
203	<	if (idm != posCsc.end() ) {
204	<	double f = whereInDistribution(digi.strip(), idm->second);
205	<	result.posCscResult.push_back( std::make_pair(rawId, f) );
202	>	mType = GoldenPattern::POSCSC;
203	>	cit = PattCore.find(mType);
204	>	if(cit==PattCore.cend()){
205	>	std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
206	>	return result;
207	>	}
208	>	auto idm = cit->second.find(rawId);
209	>	if (idm != cit->second.cend() ) {
210	>	double f = whereInDistribution(digi.strip(), idm->second);
211	>	result.myResults[mType].push_back(std::make_pair(rawId, f));
212		CSCDetId csc(rawId);
213		if (csc.station()==1) result.hasStation1 = true;
214		if (csc.station()==2) result.hasStation1 = true;
215		}
216	<	idm = bendingCsc.find(rawId);
217	<	if (idm != bendingCsc.end() ) {
216	>	mType = GoldenPattern::BENCSC;
217	>	cit = PattCore.find(mType);
218	>	if(cit==PattCore.cend()){
219	>	std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
220	>	return result;
221	>	}
222	>	idm = cit->second.find(rawId);
223	>	if (idm != cit->second.cend() ) {
224		double f = whereInDistribution(digi.pattern(), idm->second);
225	<	result.benCscResult.push_back( std::make_pair(rawId, f) );
225	>	result.myResults[mType].push_back(std::make_pair(rawId, f));
226		}
227		}
228		}
229	+
230		return result;
231		}
232
233		double GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
234		{
235	+
236		double sum_before = 0;
237		double sum_after = 0;
238		double sum_obj = 0;
239	<	for (MFreq::const_iterator im = m.begin(); im != m.end(); im++) {
239	>	for (MFreq::const_iterator im = m.begin(); im != m.end(); ++im) {
240		if (im->first  < obj) sum_before+= im->second;
241		if (im->first == obj) sum_obj = im->second;
242		if (im->first  > obj) sum_after += im->second;
243		}
244		double sum = std::max(1.,sum_before+sum_after+sum_obj );
245	+	//return sum_obj/sum; //AK
246		return (sum_before+sum_obj/2.)/sum;
247		}
248
249	<	void GoldenPattern::purge()
200	<	{
249	>	void GoldenPattern::purge(){
250
251	<	for (DetFreq::iterator idf=posRpc.begin(); idf != posRpc.end(); idf++) {
252	<	MFreq & mfreq = idf->second;
253	<	MFreq::iterator imf =  mfreq.begin();
254	<	while (imf != mfreq.end() ) {
255	<	bool remove = false;
256	<	int pos = imf->first;
257	<	unsigned int bef2 = (mfreq.find(pos-2) != mfreq.end()) ?  mfreq[pos-2] : 0;
258	<	unsigned int bef1 = (mfreq.find(pos-1) != mfreq.end()) ?  mfreq[pos-1] : 0;
259	<	unsigned int aft1 = (mfreq.find(pos+1) != mfreq.end()) ?  mfreq[pos+1] : 0;
260	<	unsigned int aft2 = (mfreq.find(pos+2) != mfreq.end()) ?  mfreq[pos+2] : 0;
261	<	unsigned int aft3 = (mfreq.find(pos+3) != mfreq.end()) ?  mfreq[pos+3] : 0;
262	<	if (mfreq[pos]==1 && bef1==0 && aft1==0) remove = true;
263	<	if (mfreq[pos]==1 && aft1==1 && aft2==0 && aft3==0 && bef1==0 && bef2==0)  remove = true;
264	<	if (mfreq[pos]==2 && aft1==0 && aft2==0 && bef1==0 && bef2==0)  remove = true;
265	<	//if(remove) std::cout<<"Cleaning pattern: "<<*this<<std::endl;
266	<	if (remove) { mfreq.erase(imf++); } else { ++imf; }
251	>	bool remove = false;
252	>	int pos;
253	>	unsigned int bef2, bef1, aft1, aft2, aft3;
254	>	for (auto isf=PattCore.begin();isf!=PattCore.end();++isf){
255	>	for (auto idf = isf->second.begin(); idf !=isf->second.end();) {
256	>	for (auto imf = idf->second.begin(); imf != idf->second.end();) {
257	>	remove = false;
258	>	pos = imf->first;
259	>	bef2 = (idf->second.find(pos-2) != idf->second.end()) ?  idf->second[pos-2] : 0;
260	>	bef1 = (idf->second.find(pos-1) != idf->second.end()) ?  idf->second[pos-1] : 0;
261	>	aft1 = (idf->second.find(pos+1) != idf->second.end()) ?  idf->second[pos+1] : 0;
262	>	aft2 = (idf->second.find(pos+2) != idf->second.end()) ?  idf->second[pos+2] : 0;
263	>	aft3 = (idf->second.find(pos+3) != idf->second.end()) ?  idf->second[pos+3] : 0;
264	>	if (idf->second[pos]==1 && bef1==0 && aft1==0) remove = true;
265	>	if (idf->second[pos]==1 && aft1==1 && aft2==0 && aft3==0 && bef1==0 && bef2==0)  remove = true;
266	>	if (idf->second[pos]==2 && aft1==0 && aft2==0 && bef1==0 && bef2==0)  remove = true;
267	>	//if(remove) std::cout<<"Cleaning pattern: "<<*this<<std::endl;
268	>	if(remove) {idf->second.erase(imf++); } else { ++imf; }
269	>	}
270	>	if (idf->second.size()==0) isf->second.erase(idf++);  else  ++idf;
271		}
272	+	break;
273		}
274	<	DetFreq::iterator idf=posRpc.begin();
221	<	while (idf != posRpc.end()) if (idf->second.size()==0) posRpc.erase(idf++);  else  ++idf;
222	<
223	<	}
274	>	}
275
276		std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {
226	–	out <<"GoldenPattern"<< o.theKey <<std::endl;
227	–	// RPC
228	–	for (GoldenPattern::DetFreq::const_iterator im = o.posRpc.begin(); im != o.posRpc.end(); im++) {
229	–	out <<"POSRPC Det: "<< im->first<<" Value: ";
230	–	for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
231	–	out << std::endl;
232	–	}
233	–	// Csc pos
234	–	for (GoldenPattern::DetFreq::const_iterator im = o.posCsc.begin(); im != o.posCsc.end(); im++) {
235	–	out <<"POSCSC Det: "<< im->first<<" Value: ";
236	–	for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
237	–	out << std::endl;
238	–	}
239	–	// Csc bending
240	–	for (GoldenPattern::DetFreq::const_iterator im = o.bendingCsc.begin(); im != o.bendingCsc.end(); im++) {
241	–	out <<"BENCSC Det: "<< im->first<<" Value: ";
242	–	for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
243	–	out << std::endl;
244	–	}
245	–	// DT pos
246	–	for (GoldenPattern::DetFreq::const_iterator im = o.posDt.begin(); im != o.posDt.end(); im++) {
247	–	out <<"POSDT Det: "<< im->first<<" Value: ";
248	–	for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
249	–	out << std::endl;
250	–	}
251	–	// DT bending
252	–	for (GoldenPattern::DetFreq::const_iterator im = o.bendingDt.begin(); im != o.bendingDt.end(); im++) {
253	–	out <<"BENDT Det: "<< im->first<<" Value: ";
254	–	for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
255	–	out << std::endl;
256	–	}
277
278	<	return out;
279	<	}
278	>	out <<"GoldenPattern"<< o.theKey <<std::endl;
279	>
280	>	std::vector<std::string> typeInfos = {"POSRPC","POSCSC","BENCSC","POSDT","BENDT"};
281
282	+	for (auto isf=o.PattCore.cbegin();isf!=o.PattCore.cend();++isf){
283	+	for (auto idf = isf->second.cbegin(); idf!=isf->second.cend();++idf) {
284	+	out <<typeInfos[isf->first]<<" Det: "<< idf->first<<" Value: ";
285	+	for (auto imf = idf->second.cbegin(); imf != idf->second.cend();++imf)
286	+	{ out << imf->first<<":"<<imf->second<<", "; }
287	+	out << std::endl;
288	+	}
289	+	}
290	+	return out;
291	+	}
292	+

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines