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

Diff Legend

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