L1RpcTriggerAnalysis/src/GoldenPattern.cc

#include <functional>
#include <numeric>

#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 mType=myResults.cbegin();mType!=myResults.cend();++mType){    
    for (auto it=mType->second.cbegin(); it!=mType->second.cend();++it) fract *= norm(mType->first,it->second); 
  }
  unsigned int nTot = 0;
  for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;    

  theValue*=(myResults[GoldenPattern::POSRPC].size()==nMatchedPoints[GoldenPattern::POSRPC]);
  theValue*=(myResults[GoldenPattern::POSDT].size()==nMatchedPoints[GoldenPattern::POSDT]);
  theValue*=(myResults[GoldenPattern::POSCSC].size()==nMatchedPoints[GoldenPattern::POSCSC]);
  theValue = ( nTot > 4) ? pow(fract, 1./((double) nTot)) : 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) ++nMatchedPoints[where];
  else normValue =1.;
  return normValue; 
}

bool GoldenPattern::Result::operator < (const GoldenPattern::Result &o) const {
  if ( *this && o) {
    if (nMatchedTot() < o.nMatchedTot()) return true;
    else if (nMatchedTot() == o.nMatchedTot() && value() < o.value()) return true;
    else return false; 
  } 
  else if (o) {return true; }
  else if (*this) { return false; } 
  else return false;

//  return (value() < o.value() ); 
//  return false;
}

GoldenPattern::Result::operator bool() const {
  return (value() > 0.1); // && hasStation1 && hasStation2);
}

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

unsigned int GoldenPattern::Result::nMatchedTot() const {
  run();
  unsigned int nTot = 0;
  for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;    
  return nTot;   
}


std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
{
  o.run();

  /*
  out <<"Result: "
      << " value: "<<o.theValue
      <<" nPos+Ben: (";
  for(auto mType=o.myResults.cbegin();
      mType!=myResults.cend();++mType){
    out<<o.nMatchedPoints[iType]<<"/"<<o.myResults[mType].size()<<", ";
  }
  out <<", tot:"<<o.nMatchedTot()<<")";
  */
  /*
  <<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()
  */
  return out;
}

void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq){
  PattCore[aCase][rawId][posOrBen] += freq; 
}


void GoldenPattern::add(const Pattern & p) {

  const Pattern::DataType & detdigi = p ;
  for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
    uint32_t rawId = is->first;
    DetId detId(rawId);
    if (detId.det() != DetId::Muon) 
      std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
    switch (detId.subdetId()) {
      case MuonSubdetId::RPC: {
        RPCDigiSpec digi(rawId, is->second);
        PattCore[GoldenPattern::POSRPC][rawId][digi.halfStrip()]++;
        break;
      }
      case MuonSubdetId::DT: {
        DTphDigiSpec digi(rawId, is->second);
        if (digi.bxNum() != 0 || digi.bxCnt() != 0 || digi.ts2() != 0) break;   
        PattCore[GoldenPattern::POSDT][rawId][digi.phi()]++;
        PattCore[GoldenPattern::BENDT][rawId][digi.phiB()]++;
        break;
      }
      case MuonSubdetId::CSC: {
        CSCDigiSpec digi(rawId, is->second);
        PattCore[GoldenPattern::POSCSC][rawId][digi.strip()]++;
        PattCore[GoldenPattern::BENCSC][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;
  uint32_t  nTot = 0;
  ///Check spatial compatibility of GP with Pattern.
  ///Require at least 5 measurements in common detUnits
  for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
    uint32_t rawId = is->first;    
    for(auto mType = PattCore.cbegin();mType!=PattCore.cend();++mType){
      nTot+=mType->second.count(rawId);
    }
  }
  if(nTot<5) return result;


  SystFreq::const_iterator cit;
  DetFreq::const_iterator idm;
  PosBenCase mType;
  for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
    uint32_t rawId = is->first;
    DetId detId(rawId);
    if (detId.det() != DetId::Muon){
      std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
      return result;
    }
    if (detId.subdetId() == MuonSubdetId::RPC) {
      RPCDigiSpec digi(rawId, is->second);
      mType = GoldenPattern::POSRPC;
      cit = PattCore.find(mType);
      if(cit==PattCore.cend()){
        std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
        return result;
      }
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        double f = whereInDistribution(digi.halfStrip(), idm->second);
        result.myResults[mType].push_back(std::make_pair(rawId, f)); 
        RPCDetId rpc(rawId);
        if(rpc.station()==1) result.hasStation1 = true;
        if(rpc.station()==2) result.hasStation2 = true;
      }
    }
    else if (detId.subdetId() == MuonSubdetId::DT) {
      DTphDigiSpec digi(rawId, is->second);
      mType = GoldenPattern::POSDT;
      cit = PattCore.find(mType);
      if(cit==PattCore.cend()){
        std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
        return result;
      }
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        double f = whereInDistribution(digi.phi(), idm->second);
        result.myResults[mType].push_back(std::make_pair(rawId, f)); 
        DTChamberId dt(rawId);
        if(dt.station()==1) result.hasStation1 = true;
        if(dt.station()==2) result.hasStation2 = true;
      }
      mType = GoldenPattern::BENDT;
      cit = PattCore.find(mType);
      if(cit==PattCore.cend()){
        std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
        return result;
      }
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        double f = whereInDistribution(digi.phiB(), idm->second);
        result.myResults[mType].push_back(std::make_pair(rawId, f));
      }
    }
    else if (detId.subdetId() == MuonSubdetId::CSC) {
      CSCDigiSpec digi(rawId, is->second);
      mType = GoldenPattern::POSCSC;
      cit = PattCore.find(mType);
      if(cit==PattCore.cend()){
        std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
        return result;
      }
      auto idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        double f = whereInDistribution(digi.strip(), idm->second);
        result.myResults[mType].push_back(std::make_pair(rawId, f)); 
        CSCDetId csc(rawId);
        if (csc.station()==1) result.hasStation1 = true;
        if (csc.station()==2) result.hasStation1 = true;
      }
      mType = GoldenPattern::BENCSC;
      cit = PattCore.find(mType);
      if(cit==PattCore.cend()){
        std::cout << "PROBLEM: unknown measurement type: "<<mType<<std::endl;
        return result;
      }
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        double f = whereInDistribution(digi.pattern(), idm->second);
        result.myResults[mType].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_obj/sum; //AK
  return (sum_before+sum_obj/2.)/sum; 
}

void GoldenPattern::purge(){

  bool remove = false;
  int pos;
  unsigned int bef2, bef1, aft1, aft2, aft3;
  for (auto isf=PattCore.begin();isf!=PattCore.end();++isf){
    for (auto idf = isf->second.begin(); idf !=isf->second.end();) {
      for (auto imf = idf->second.begin(); imf != idf->second.end();) {
        remove = false;
        pos = imf->first;  
        bef2 = (idf->second.find(pos-2) != idf->second.end()) ?  idf->second[pos-2] : 0;  
        bef1 = (idf->second.find(pos-1) != idf->second.end()) ?  idf->second[pos-1] : 0;  
        aft1 = (idf->second.find(pos+1) != idf->second.end()) ?  idf->second[pos+1] : 0;  
        aft2 = (idf->second.find(pos+2) != idf->second.end()) ?  idf->second[pos+2] : 0;  
        aft3 = (idf->second.find(pos+3) != idf->second.end()) ?  idf->second[pos+3] : 0; 
        if (idf->second[pos]==1 && bef1==0 && aft1==0) remove = true;
        if (idf->second[pos]==1 && aft1==1 && aft2==0 && aft3==0 && bef1==0 && bef2==0)  remove = true;
        if (idf->second[pos]==2 && aft1==0 && aft2==0 && bef1==0 && bef2==0)  remove = true;
        //if(remove) std::cout<<"Cleaning pattern: "<<*this<<std::endl;      
        if(remove) {idf->second.erase(imf++); } else { ++imf; } 
      }
      if (idf->second.size()==0) isf->second.erase(idf++);  else  ++idf;
    }
    break;
  }
}

std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {

 out <<"GoldenPattern"<< o.theKey <<std::endl;

 std::vector<std::string> typeInfos = {"POSRPC","POSCSC","BENCSC","POSDT","BENDT"};

 for (auto isf=o.PattCore.cbegin();isf!=o.PattCore.cend();++isf){
   for (auto idf = isf->second.cbegin(); idf!=isf->second.cend();++idf) {      
     out <<typeInfos[isf->first]<<" Det: "<< idf->first<<" Value: ";
     for (auto imf = idf->second.cbegin(); imf != idf->second.cend();++imf) 
       { out << imf->first<<":"<<imf->second<<", "; }
     out << std::endl;
   }
 }
 return out;
}
 
Revision:	1.6
Committed:	Tue Jun 25 11:06:02 2013 UTC (11 years, 10 months ago) by akalinow
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+159 -127 lines
Log Message:	*some further interface modifications by AK
#	User	Rev	Content
1	akalinow	1.6	#include <functional>
2			#include <numeric>
3
4	konec	1.1	#include "UserCode/L1RpcTriggerAnalysis/interface/GoldenPattern.h"
5
6			#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
7			#include "DataFormats/MuonDetId/interface/RPCDetId.h"
8			#include "DataFormats/MuonDetId/interface/CSCDetId.h"
9			#include "DataFormats/MuonDetId/interface/DTChamberId.h"
10
11			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
12			#include "UserCode/L1RpcTriggerAnalysis/interface/DTphDigiSpec.h"
13			#include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
14			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"
15
16	akalinow	1.6	void GoldenPattern::Result::runNoCheck() const {
17			double fract = 1;
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	konec	1.1
25	akalinow	1.6	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	konec	1.4	theValue = ( nTot > 4) ? pow(fract, 1./((double) nTot)) : 0.;
29	konec	1.1	}
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	akalinow	1.6	if (normValue > epsilon) ++nMatchedPoints[where];
35			else normValue =1.;
36	konec	1.1	return normValue;
37			}
38
39			bool GoldenPattern::Result::operator < (const GoldenPattern::Result &o) const {
40	konec	1.2	if ( *this && o) {
41			if (nMatchedTot() < o.nMatchedTot()) return true;
42			else if (nMatchedTot() == o.nMatchedTot() && value() < o.value()) return true;
43			else return false;
44			}
45			else if (o) {return true; }
46			else if (*this) { return false; }
47			else return false;
48
49			// return (value() < o.value() );
50			// return false;
51	konec	1.1	}
52
53			GoldenPattern::Result::operator bool() const {
54	konec	1.3	return (value() > 0.1); // && hasStation1 && hasStation2);
55	konec	1.1	}
56
57			double GoldenPattern::Result::value() const {
58			run();
59			return theValue;
60			}
61
62			unsigned int GoldenPattern::Result::nMatchedTot() const {
63			run();
64	akalinow	1.6	unsigned int nTot = 0;
65			for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;
66			return nTot;
67	konec	1.1	}
68
69	akalinow	1.5
70	konec	1.1	std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
71			{
72			o.run();
73	akalinow	1.6
74			/*
75	konec	1.1	out <<"Result: "
76			<< " value: "<<o.theValue
77	akalinow	1.6	<<" 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	konec	1.1	return out;
90			}
91
92	akalinow	1.6	void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq){
93			PattCore[aCase][rawId][posOrBen] += freq;
94			}
95	konec	1.1
96
97	akalinow	1.6	void GoldenPattern::add(const Pattern & p) {
98	konec	1.1
99			const Pattern::DataType & detdigi = p ;
100	akalinow	1.6	for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
101	konec	1.1	uint32_t rawId = is->first;
102			DetId detId(rawId);
103	akalinow	1.6	if (detId.det() != DetId::Muon)
104			std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
105	konec	1.1	switch (detId.subdetId()) {
106			case MuonSubdetId::RPC: {
107			RPCDigiSpec digi(rawId, is->second);
108	akalinow	1.6	PattCore[GoldenPattern::POSRPC][rawId][digi.halfStrip()]++;
109	konec	1.1	break;
110			}
111			case MuonSubdetId::DT: {
112			DTphDigiSpec digi(rawId, is->second);
113	akalinow	1.6	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	konec	1.1	break;
117			}
118			case MuonSubdetId::CSC: {
119			CSCDigiSpec digi(rawId, is->second);
120	akalinow	1.6	PattCore[GoldenPattern::POSCSC][rawId][digi.strip()]++;
121			PattCore[GoldenPattern::BENCSC][rawId][digi.pattern()]++;
122	konec	1.1	break;
123			}
124			default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
125			};
126			}
127			}
128
129	akalinow	1.6	GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const
130	konec	1.1	{
131			Result result;
132	akalinow	1.6	const Pattern::DataType & detdigi = p;
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	akalinow	1.5
145	akalinow	1.6	SystFreq::const_iterator cit;
146			DetFreq::const_iterator idm;
147			PosBenCase mType;
148			for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
149	konec	1.1	uint32_t rawId = is->first;
150			DetId detId(rawId);
151	akalinow	1.6	if (detId.det() != DetId::Muon){
152			std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
153			return result;
154			}
155	konec	1.1	if (detId.subdetId() == MuonSubdetId::RPC) {
156			RPCDigiSpec digi(rawId, is->second);
157	akalinow	1.6	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	konec	1.3	double f = whereInDistribution(digi.halfStrip(), idm->second);
166	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
167	konec	1.3	RPCDetId rpc(rawId);
168			if(rpc.station()==1) result.hasStation1 = true;
169			if(rpc.station()==2) result.hasStation2 = true;
170	konec	1.1	}
171	akalinow	1.6	}
172			else if (detId.subdetId() == MuonSubdetId::DT) {
173	konec	1.1	DTphDigiSpec digi(rawId, is->second);
174	akalinow	1.6	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	konec	1.1	double f = whereInDistribution(digi.phi(), idm->second);
183	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
184	konec	1.3	DTChamberId dt(rawId);
185			if(dt.station()==1) result.hasStation1 = true;
186			if(dt.station()==2) result.hasStation2 = true;
187	konec	1.1	}
188	akalinow	1.6	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	konec	1.1	}
199	akalinow	1.6	}
200			else if (detId.subdetId() == MuonSubdetId::CSC) {
201	konec	1.1	CSCDigiSpec digi(rawId, is->second);
202	akalinow	1.6	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	konec	1.3	CSCDetId csc(rawId);
213			if (csc.station()==1) result.hasStation1 = true;
214			if (csc.station()==2) result.hasStation1 = true;
215	konec	1.1	}
216	akalinow	1.6	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	konec	1.1	double f = whereInDistribution(digi.pattern(), idm->second);
225	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
226	konec	1.1	}
227			}
228			}
229	akalinow	1.6
230	konec	1.1	return result;
231			}
232
233			double GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
234			{
235	akalinow	1.6
236	konec	1.1	double sum_before = 0;
237			double sum_after = 0;
238			double sum_obj = 0;
239	akalinow	1.6	for (MFreq::const_iterator im = m.begin(); im != m.end(); ++im) {
240	konec	1.1	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	akalinow	1.6	//return sum_obj/sum; //AK
246	konec	1.1	return (sum_before+sum_obj/2.)/sum;
247			}
248
249	akalinow	1.6	void GoldenPattern::purge(){
250	akalinow	1.5
251	akalinow	1.6	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	konec	1.1	}
272	akalinow	1.6	break;
273	konec	1.1	}
274	akalinow	1.6	}
275	konec	1.1
276			std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {
277	akalinow	1.6
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	konec	1.1	out << std::endl;
288			}
289	akalinow	1.6	}
290			return out;
291			}
292