L1RpcTriggerAnalysis/src/GoldenPattern.cc

#include "UserCode/L1RpcTriggerAnalysis/interface/GoldenPattern.h"

#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "DataFormats/MuonDetId/interface/RPCDetId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h"

#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
#include "UserCode/L1RpcTriggerAnalysis/interface/DTphDigiSpec.h"
#include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"


void GoldenPattern::Result::runNoCheck() const 
{
  double fract = 1;
  for (auto i=posRpcResult.begin(); i!=posRpcResult.end();i++) fract *= norm(POSRPC,i->second); 
  for (auto i=posCscResult.begin(); i!=posCscResult.end();i++) fract *= norm(POSCSC,i->second); 
  for (auto i=posDtResult.begin();  i!=posDtResult.end();i++)  fract *= norm(POSDT, i->second);
  for (auto i=benCscResult.begin(); i!=benCscResult.end();i++) fract *= norm(BENCSC,i->second);
  for (auto i=benDtResult.begin();  i!=benDtResult.end();i++)  fract *= norm(BENDT, i->second);
  unsigned int nTot = nMatchedPosRpc+nMatchedPosCsc+nMatchedPosDt+nMatchedBenCsc+nMatchedBenDt;
  theValue = ( nTot > 2) ? pow(fract, 1./((double) nTot)) : 0.;
  if (posRpcResult.size() != nMatchedPosRpc) theValue = 0.;
  if (posCscResult.size() != nMatchedPosCsc) theValue = 0.;
  if (posDtResult.size()  != nMatchedPosDt)  theValue = 0.;
}

double GoldenPattern::Result::norm(GoldenPattern::PosBenCase where, double whereInDist) const {
  double normValue = 2.*(0.5-fabs(whereInDist-0.5));   
  static const double epsilon = 1.e-9;
  if (normValue > epsilon) {
    switch (where) {
      case POSRPC : nMatchedPosRpc++; break;
      case POSCSC : nMatchedPosCsc++; break;
      case BENCSC : nMatchedBenCsc++; break;
      case POSDT  : nMatchedPosDt++;  break;
      case BENDT  : nMatchedBenDt++;  break;
    };
  } else {
    normValue =1.;
  }
  return normValue; 
}

bool GoldenPattern::Result::operator < (const GoldenPattern::Result &o) const {
  return (value() < o.value() ); 
}

GoldenPattern::Result::operator bool() const {
  return (value() > 0.1);
}

double GoldenPattern::Result::value() const { 
  run(); 
  return theValue; 
}

unsigned int GoldenPattern::Result::nMatchedTot() const {
  run();
  return nMatchedPosRpc+nMatchedPosCsc+nMatchedBenCsc+nMatchedPosDt+nMatchedBenDt;
}

std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
{
  o.run();
  out <<"Result: "
      << " value: "<<o.theValue
      <<" nPos+Ben: ("<<o.nMatchedPosCsc<<"/"<<o.posCscResult.size()<<"+"<<o.nMatchedBenCsc<<"/"<<o.benCscResult.size()
                <<", "<<o.nMatchedPosDt <<"/"<<o.posDtResult.size()<<"+"<<o.nMatchedBenDt<<"/"<<o.benDtResult.size()
                <<", "<<o.nMatchedPosRpc<<"/"<<o.posRpcResult.size()<<", tot:"<<o.nMatchedTot()<<")";
  return out;
}


void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq)
{
  switch ( aCase ) {
    case POSRPC : posRpc[rawId][posOrBen]     += freq; break;
    case POSCSC : posCsc[rawId][posOrBen]     += freq; break;
    case BENCSC : bendingCsc[rawId][posOrBen] += freq; break;
    case POSDT  : posDt[rawId][posOrBen]      += freq; break;
    case BENDT  : bendingDt[rawId][posOrBen]  += freq; break;
  };
}


void GoldenPattern::add(const Pattern & p) 
{
  const Pattern::DataType & detdigi = p ;
  for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
    uint32_t rawId = is->first;
    DetId detId(rawId);
    if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
    switch (detId.subdetId()) {
      case MuonSubdetId::RPC: {
        RPCDigiSpec digi(rawId, is->second);
        posRpc[rawId][digi.halfStrip()]++;
        break;
      }
      case MuonSubdetId::DT: {
        DTphDigiSpec digi(rawId, is->second);
        if (digi.bxNum() != 0 || digi.bxCnt() != 0 || digi.ts2() != 0) break;
        posDt[rawId][digi.phi()]++;
        bendingDt[rawId][digi.phiB()]++;
        break;
      }
      case MuonSubdetId::CSC: {
        CSCDigiSpec digi(rawId, is->second);
        posCsc[rawId][digi.strip()]++;
        bendingCsc[rawId][digi.pattern()]++;
        break;
      }
      default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
    };
  } 
}

GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const 
{
  Result result;
  const Pattern::DataType & detdigi = p;
  for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
    uint32_t rawId = is->first;
    DetId detId(rawId);
    if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
    if (detId.subdetId() == MuonSubdetId::RPC) {
      RPCDigiSpec digi(rawId, is->second);
      DetFreq::const_iterator idm = posRpc.find(rawId);
      if (idm != posRpc.end() ) {
          double f = whereInDistribution(digi.halfStrip(), idm->second);
          result.posRpcResult.push_back( std::make_pair(rawId, f) ); 
      }
    } else if (detId.subdetId() == MuonSubdetId::DT) {
      DTphDigiSpec digi(rawId, is->second);
      DetFreq::const_iterator idm = posDt.find(rawId);
      if (idm != posDt.end() ) {
        double f = whereInDistribution(digi.phi(), idm->second);
        result.posDtResult.push_back( std::make_pair(rawId, f) );
      }
      idm = bendingDt.find(rawId);
      if (idm != bendingDt.end() ) {
        double f = whereInDistribution(digi.phiB(), idm->second);
        result.benDtResult.push_back( std::make_pair(rawId, f) );
      }
    } else if (detId.subdetId() == MuonSubdetId::CSC) {
      CSCDigiSpec digi(rawId, is->second);
      DetFreq::const_iterator idm = posCsc.find(rawId);
      if (idm != posCsc.end() ) {
        double f = whereInDistribution(digi.strip(), idm->second);
        result.posCscResult.push_back( std::make_pair(rawId, f) );
      }
      idm = bendingCsc.find(rawId);
      if (idm != bendingCsc.end() ) {
        double f = whereInDistribution(digi.pattern(), idm->second);
        result.benCscResult.push_back( std::make_pair(rawId, f) );
      }
    }
  }
  return result;
}

double GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
{
  double sum_before = 0;
  double sum_after = 0;
  double sum_obj = 0;
  for (MFreq::const_iterator im = m.begin(); im != m.end(); im++) {
    if (im->first  < obj) sum_before+= im->second;
    if (im->first == obj) sum_obj = im->second;  
    if (im->first  > obj) sum_after += im->second;
  }
  double sum = std::max(1.,sum_before+sum_after+sum_obj );
  return (sum_before+sum_obj/2.)/sum; 
}

void GoldenPattern::purge()
{
  for (DetFreq::iterator idf=posRpc.begin(); idf != posRpc.end(); idf++) {
    MFreq & mfreq = idf->second;
    MFreq::iterator imf =  mfreq.begin();
    while (imf != mfreq.end() ) {
      bool remove = false;
      int pos = imf->first;  
      unsigned int bef2 = (mfreq.find(pos-2) != mfreq.end()) ?  mfreq[pos-2] : 0;  
      unsigned int bef1 = (mfreq.find(pos-1) != mfreq.end()) ?  mfreq[pos-1] : 0;  
      unsigned int aft1 = (mfreq.find(pos+1) != mfreq.end()) ?  mfreq[pos+1] : 0;  
      unsigned int aft2 = (mfreq.find(pos+2) != mfreq.end()) ?  mfreq[pos+2] : 0;  
      unsigned int aft3 = (mfreq.find(pos+3) != mfreq.end()) ?  mfreq[pos+3] : 0;  
      if (mfreq[pos]==1 && bef1==0 && aft1==0) remove = true;
      if (mfreq[pos]==1 && aft1==1 && aft2==0 && aft3==0 && bef1==0 && bef2==0)  remove = true;
      if (mfreq[pos]==2 && aft1==0 && aft2==0 && bef1==0 && bef2==0)  remove = true;
      if (remove) { mfreq.erase(imf++); } else { ++imf; }
    }
  }
  DetFreq::iterator idf=posRpc.begin();  
  while (idf != posRpc.end()) if (idf->second.size()==0) posRpc.erase(idf++);  else  ++idf; 
 
}

std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {
   out <<"GoldenPattern"<< o.theKey <<std::endl;
   // RPC
   for (GoldenPattern::DetFreq::const_iterator im = o.posRpc.begin(); im != o.posRpc.end(); im++) {
     out <<"RPC Det:"<< im->first<<" Pos: ";
     for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
     out << std::endl;
   }
   // DT pos
   for (GoldenPattern::DetFreq::const_iterator im = o.posDt.begin(); im != o.posDt.end(); im++) {
     out <<"DT  Det:"<< im->first<<" Pos: ";
     for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
     out << std::endl;
   }
   // DT bending 
   for (GoldenPattern::DetFreq::const_iterator im = o.bendingDt.begin(); im != o.bendingDt.end(); im++) {
     out <<"DT  Det:"<< im->first<<" Ben: ";
     for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
     out << std::endl;
   }

   // Csc pos
   for (GoldenPattern::DetFreq::const_iterator im = o.posCsc.begin(); im != o.posCsc.end(); im++) {
     out <<"Csc  Det:"<< im->first<<" Pos: ";
     for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
     out << std::endl;
   }
   // Csc bending 
   for (GoldenPattern::DetFreq::const_iterator im = o.bendingCsc.begin(); im != o.bendingCsc.end(); im++) {
     out <<"Csc  Det:"<< im->first<<" Ben: ";
     for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it !=  im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
     out << std::endl;
   }

   return out;
  }

Revision:	1.1
Committed:	Fri May 17 13:04:54 2013 UTC (11 years, 11 months ago) by konec
Content type:	text/plain
Branch:	MAIN
Log Message:	extenstion for OTF
#	User	Rev	Content
1	konec	1.1	#include "UserCode/L1RpcTriggerAnalysis/interface/GoldenPattern.h"
2
3			#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
4			#include "DataFormats/MuonDetId/interface/RPCDetId.h"
5			#include "DataFormats/MuonDetId/interface/CSCDetId.h"
6			#include "DataFormats/MuonDetId/interface/DTChamberId.h"
7
8			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
9			#include "UserCode/L1RpcTriggerAnalysis/interface/DTphDigiSpec.h"
10			#include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
11			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"
12
13
14			void GoldenPattern::Result::runNoCheck() const
15			{
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.;
42			}
43			return normValue;
44			}
45
46			bool GoldenPattern::Result::operator < (const GoldenPattern::Result &o) const {
47			return (value() < o.value() );
48			}
49
50			GoldenPattern::Result::operator bool() const {
51			return (value() > 0.1);
52			}
53
54			double GoldenPattern::Result::value() const {
55			run();
56			return theValue;
57			}
58
59			unsigned int GoldenPattern::Result::nMatchedTot() const {
60			run();
61			return nMatchedPosRpc+nMatchedPosCsc+nMatchedBenCsc+nMatchedPosDt+nMatchedBenDt;
62			}
63
64			std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
65			{
66			o.run();
67			out <<"Result: "
68			<< " value: "<<o.theValue
69			<<" nPos+Ben: ("<<o.nMatchedPosCsc<<"/"<<o.posCscResult.size()<<"+"<<o.nMatchedBenCsc<<"/"<<o.benCscResult.size()
70			<<", "<<o.nMatchedPosDt <<"/"<<o.posDtResult.size()<<"+"<<o.nMatchedBenDt<<"/"<<o.benDtResult.size()
71			<<", "<<o.nMatchedPosRpc<<"/"<<o.posRpcResult.size()<<", tot:"<<o.nMatchedTot()<<")";
72			return out;
73			}
74
75
76
77			void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq)
78			{
79			switch ( aCase ) {
80			case POSRPC : posRpc[rawId][posOrBen] += freq; break;
81			case POSCSC : posCsc[rawId][posOrBen] += freq; break;
82			case BENCSC : bendingCsc[rawId][posOrBen] += freq; break;
83			case POSDT : posDt[rawId][posOrBen] += freq; break;
84			case BENDT : bendingDt[rawId][posOrBen] += freq; break;
85			};
86			}
87
88
89			void GoldenPattern::add(const Pattern & p)
90			{
91			const Pattern::DataType & detdigi = p ;
92			for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
93			uint32_t rawId = is->first;
94			DetId detId(rawId);
95			if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
96			switch (detId.subdetId()) {
97			case MuonSubdetId::RPC: {
98			RPCDigiSpec digi(rawId, is->second);
99			posRpc[rawId][digi.halfStrip()]++;
100			break;
101			}
102			case MuonSubdetId::DT: {
103			DTphDigiSpec digi(rawId, is->second);
104			if (digi.bxNum() != 0 \|\| digi.bxCnt() != 0 \|\| digi.ts2() != 0) break;
105			posDt[rawId][digi.phi()]++;
106			bendingDt[rawId][digi.phiB()]++;
107			break;
108			}
109			case MuonSubdetId::CSC: {
110			CSCDigiSpec digi(rawId, is->second);
111			posCsc[rawId][digi.strip()]++;
112			bendingCsc[rawId][digi.pattern()]++;
113			break;
114			}
115			default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
116			};
117			}
118			}
119
120			GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const
121			{
122			Result result;
123			const Pattern::DataType & detdigi = p;
124			for (Pattern::DataType::const_iterator is = detdigi.begin(); is != detdigi.end(); ++is) {
125			uint32_t rawId = is->first;
126			DetId detId(rawId);
127			if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
128			if (detId.subdetId() == MuonSubdetId::RPC) {
129			RPCDigiSpec digi(rawId, is->second);
130			DetFreq::const_iterator idm = posRpc.find(rawId);
131			if (idm != posRpc.end() ) {
132			double f = whereInDistribution(digi.halfStrip(), idm->second);
133			result.posRpcResult.push_back( std::make_pair(rawId, f) );
134			}
135			} else if (detId.subdetId() == MuonSubdetId::DT) {
136			DTphDigiSpec digi(rawId, is->second);
137			DetFreq::const_iterator idm = posDt.find(rawId);
138			if (idm != posDt.end() ) {
139			double f = whereInDistribution(digi.phi(), idm->second);
140			result.posDtResult.push_back( std::make_pair(rawId, f) );
141			}
142			idm = bendingDt.find(rawId);
143			if (idm != bendingDt.end() ) {
144			double f = whereInDistribution(digi.phiB(), idm->second);
145			result.benDtResult.push_back( std::make_pair(rawId, f) );
146			}
147			} else if (detId.subdetId() == MuonSubdetId::CSC) {
148			CSCDigiSpec digi(rawId, is->second);
149			DetFreq::const_iterator idm = posCsc.find(rawId);
150			if (idm != posCsc.end() ) {
151			double f = whereInDistribution(digi.strip(), idm->second);
152			result.posCscResult.push_back( std::make_pair(rawId, f) );
153			}
154			idm = bendingCsc.find(rawId);
155			if (idm != bendingCsc.end() ) {
156			double f = whereInDistribution(digi.pattern(), idm->second);
157			result.benCscResult.push_back( std::make_pair(rawId, f) );
158			}
159			}
160			}
161			return result;
162			}
163
164			double GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
165			{
166			double sum_before = 0;
167			double sum_after = 0;
168			double sum_obj = 0;
169			for (MFreq::const_iterator im = m.begin(); im != m.end(); im++) {
170			if (im->first < obj) sum_before+= im->second;
171			if (im->first == obj) sum_obj = im->second;
172			if (im->first > obj) sum_after += im->second;
173			}
174			double sum = std::max(1.,sum_before+sum_after+sum_obj );
175			return (sum_before+sum_obj/2.)/sum;
176			}
177
178			void GoldenPattern::purge()
179			{
180			for (DetFreq::iterator idf=posRpc.begin(); idf != posRpc.end(); idf++) {
181			MFreq & mfreq = idf->second;
182			MFreq::iterator imf = mfreq.begin();
183			while (imf != mfreq.end() ) {
184			bool remove = false;
185			int pos = imf->first;
186			unsigned int bef2 = (mfreq.find(pos-2) != mfreq.end()) ? mfreq[pos-2] : 0;
187			unsigned int bef1 = (mfreq.find(pos-1) != mfreq.end()) ? mfreq[pos-1] : 0;
188			unsigned int aft1 = (mfreq.find(pos+1) != mfreq.end()) ? mfreq[pos+1] : 0;
189			unsigned int aft2 = (mfreq.find(pos+2) != mfreq.end()) ? mfreq[pos+2] : 0;
190			unsigned int aft3 = (mfreq.find(pos+3) != mfreq.end()) ? mfreq[pos+3] : 0;
191			if (mfreq[pos]==1 && bef1==0 && aft1==0) remove = true;
192			if (mfreq[pos]==1 && aft1==1 && aft2==0 && aft3==0 && bef1==0 && bef2==0) remove = true;
193			if (mfreq[pos]==2 && aft1==0 && aft2==0 && bef1==0 && bef2==0) remove = true;
194			if (remove) { mfreq.erase(imf++); } else { ++imf; }
195			}
196			}
197			DetFreq::iterator idf=posRpc.begin();
198			while (idf != posRpc.end()) if (idf->second.size()==0) posRpc.erase(idf++); else ++idf;
199
200			}
201
202			std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {
203			out <<"GoldenPattern"<< o.theKey <<std::endl;
204			// RPC
205			for (GoldenPattern::DetFreq::const_iterator im = o.posRpc.begin(); im != o.posRpc.end(); im++) {
206			out <<"RPC Det:"<< im->first<<" Pos: ";
207			for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it != im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
208			out << std::endl;
209			}
210			// DT pos
211			for (GoldenPattern::DetFreq::const_iterator im = o.posDt.begin(); im != o.posDt.end(); im++) {
212			out <<"DT Det:"<< im->first<<" Pos: ";
213			for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it != im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
214			out << std::endl;
215			}
216			// DT bending
217			for (GoldenPattern::DetFreq::const_iterator im = o.bendingDt.begin(); im != o.bendingDt.end(); im++) {
218			out <<"DT Det:"<< im->first<<" Ben: ";
219			for (GoldenPattern::MFreq::const_iterator it = im->second.begin(); it != im->second.end(); it++) { out << it->first<<":"<<it->second<<", "; }
220			out << std::endl;
221			}
222
223			// Csc pos
224			for (GoldenPattern::DetFreq::const_iterator im = o.posCsc.begin(); im != o.posCsc.end(); im++) {
225			out <<"Csc 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			// Csc bending
230			for (GoldenPattern::DetFreq::const_iterator im = o.bendingCsc.begin(); im != o.bendingCsc.end(); im++) {
231			out <<"Csc Det:"<< im->first<<" Ben: ";
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
236			return out;
237			}
238