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

Diff Legend

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