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 {
  float 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./((float) nTot)) : 0.;
}

float GoldenPattern::Result::norm(GoldenPattern::PosBenCase where, float whereInDist) const {
  float normValue = 2.*(0.5-fabs(whereInDist-0.5));   
  //normValue = whereInDist; //AK
  static const float epsilon = 1.e-9;
  if (normValue > epsilon) ++nMatchedPoints[where];
  else normValue = 0.05;
  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);
}

float 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 cit=o.myResults.cbegin();cit!=o.myResults.cend();++cit){
    out<<o.nMatchedPoints[cit->first]<<"/"<<cit->second.size()<<", ";
  }
  out <<", tot:"<<o.nMatchedTot()<<")";
  
  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 << "GoldenPattern::compare 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())  return result; //AK: Ugly, FIX.
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        float 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()) return result;
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        float 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()) return result;
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        float 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()) return result;
      auto idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        float 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()) return result;
      idm = cit->second.find(rawId);
      if (idm != cit->second.cend() ) {
        float f = whereInDistribution(digi.pattern(), idm->second);
        result.myResults[mType].push_back(std::make_pair(rawId, f));
      }
    }
  }

  return result;
}

float GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
{

  float sum_before = 0;
  float sum_after = 0;
  float 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;
  }
  float sum = std::max((float)1.,sum_before+sum_after+sum_obj );
  //return sum_obj/sum; //AK
  return (sum_before+sum_obj/2.)/sum; 
}

bool GoldenPattern::purge(){

  bool remove = false;
  int pos;
  unsigned int bef2, bef1, aft1, aft2, aft3;
  for (auto isf=PattCore.begin();isf!=PattCore.end();){
    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;
    }
      if (isf->second.size()==0) PattCore.erase(isf++);  else  ++isf;
  }
  ///Usefull pattern has at least 4 measurements and has a RPC measurement
  return PattCore.find(POSRPC)!=PattCore.end() && PattCore.size()>4;
}

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;
     if(typeInfos[isf->first].find("RPC")!=std::string::npos){
       RPCDetId rpc(idf->first);
       out<<" ("<<rpc<<") ";
     }
     if(typeInfos[isf->first].find("CSC")!=std::string::npos){
       CSCDetId csc(idf->first);
       out<<" ("<<csc<<") ";
     }
     if(typeInfos[isf->first].find("DT")!=std::string::npos){
       DTChamberId dt(idf->first);
       out<<" ("<<dt<<") ";
     }
     out <<" Value: ";
     for (auto imf = idf->second.cbegin(); imf != idf->second.cend();++imf) 
       { out << imf->first<<":"<<imf->second<<", "; }
     out << std::endl;
   }
 }
 return out;
}
 
Revision:	1.7
Committed:	Fri Jun 28 08:31:12 2013 UTC (11 years, 10 months ago) by akalinow
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Artur_28_06_2013
Changes since 1.6:	+47 -53 lines
Log Message:	*updates 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	akalinow	1.7	float fract = 1;
18	akalinow	1.6
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	akalinow	1.7	theValue = ( nTot > 4) ? pow(fract, 1./((float) nTot)) : 0.;
29	konec	1.1	}
30
31	akalinow	1.7	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	akalinow	1.6	if (normValue > epsilon) ++nMatchedPoints[where];
36	akalinow	1.7	else normValue = 0.05;
37	konec	1.1	return normValue;
38			}
39
40			bool GoldenPattern::Result::operator < (const GoldenPattern::Result &o) const {
41	konec	1.2	if ( *this && o) {
42			if (nMatchedTot() < o.nMatchedTot()) return true;
43			else if (nMatchedTot() == o.nMatchedTot() && value() < o.value()) return true;
44			else return false;
45			}
46			else if (o) {return true; }
47			else if (*this) { return false; }
48			else return false;
49
50			// return (value() < o.value() );
51			// return false;
52	konec	1.1	}
53
54			GoldenPattern::Result::operator bool() const {
55	konec	1.3	return (value() > 0.1); // && hasStation1 && hasStation2);
56	konec	1.1	}
57
58	akalinow	1.7	float GoldenPattern::Result::value() const {
59	konec	1.1	run();
60			return theValue;
61			}
62
63			unsigned int GoldenPattern::Result::nMatchedTot() const {
64			run();
65	akalinow	1.6	unsigned int nTot = 0;
66			for(auto it=nMatchedPoints.cbegin();it!=nMatchedPoints.cend();++it) nTot+=it->second;
67			return nTot;
68	konec	1.1	}
69
70	akalinow	1.5
71	konec	1.1	std::ostream & operator << (std::ostream &out, const GoldenPattern::Result& o)
72			{
73			o.run();
74	akalinow	1.7
75	konec	1.1	out <<"Result: "
76			<< " value: "<<o.theValue
77	akalinow	1.6	<<" nPos+Ben: (";
78	akalinow	1.7	for(auto cit=o.myResults.cbegin();cit!=o.myResults.cend();++cit){
79			out<<o.nMatchedPoints[cit->first]<<"/"<<cit->second.size()<<", ";
80	akalinow	1.6	}
81			out <<", tot:"<<o.nMatchedTot()<<")";
82	akalinow	1.7
83	konec	1.1	return out;
84			}
85
86	akalinow	1.6	void GoldenPattern::add( GoldenPattern::PosBenCase aCase, uint32_t rawId, int posOrBen, unsigned int freq){
87			PattCore[aCase][rawId][posOrBen] += freq;
88			}
89	konec	1.1
90
91	akalinow	1.6	void GoldenPattern::add(const Pattern & p) {
92	konec	1.1
93			const Pattern::DataType & detdigi = p ;
94	akalinow	1.6	for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
95	konec	1.1	uint32_t rawId = is->first;
96			DetId detId(rawId);
97	akalinow	1.6	if (detId.det() != DetId::Muon)
98			std::cout << "PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
99	konec	1.1	switch (detId.subdetId()) {
100			case MuonSubdetId::RPC: {
101			RPCDigiSpec digi(rawId, is->second);
102	akalinow	1.6	PattCore[GoldenPattern::POSRPC][rawId][digi.halfStrip()]++;
103	konec	1.1	break;
104			}
105			case MuonSubdetId::DT: {
106			DTphDigiSpec digi(rawId, is->second);
107	akalinow	1.6	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	konec	1.1	break;
111			}
112			case MuonSubdetId::CSC: {
113			CSCDigiSpec digi(rawId, is->second);
114	akalinow	1.6	PattCore[GoldenPattern::POSCSC][rawId][digi.strip()]++;
115			PattCore[GoldenPattern::BENCSC][rawId][digi.pattern()]++;
116	konec	1.1	break;
117			}
118			default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
119			};
120			}
121			}
122
123	akalinow	1.6	GoldenPattern::Result GoldenPattern::compare(const Pattern &p) const
124	konec	1.1	{
125			Result result;
126	akalinow	1.6	const Pattern::DataType & detdigi = p;
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	akalinow	1.5
139	akalinow	1.6	SystFreq::const_iterator cit;
140			DetFreq::const_iterator idm;
141			PosBenCase mType;
142			for (auto is = detdigi.cbegin(); is != detdigi.cend(); ++is) {
143	konec	1.1	uint32_t rawId = is->first;
144			DetId detId(rawId);
145	akalinow	1.6	if (detId.det() != DetId::Muon){
146	akalinow	1.7	std::cout << "GoldenPattern::compare PROBLEM: hit in unknown Det, detID: "<<detId.det()<<std::endl;
147	akalinow	1.6	return result;
148			}
149	konec	1.1	if (detId.subdetId() == MuonSubdetId::RPC) {
150			RPCDigiSpec digi(rawId, is->second);
151	akalinow	1.6	mType = GoldenPattern::POSRPC;
152			cit = PattCore.find(mType);
153	akalinow	1.7	if(cit==PattCore.cend()) return result; //AK: Ugly, FIX.
154	akalinow	1.6	idm = cit->second.find(rawId);
155			if (idm != cit->second.cend() ) {
156	akalinow	1.7	float f = whereInDistribution(digi.halfStrip(), idm->second);
157	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
158	konec	1.3	RPCDetId rpc(rawId);
159			if(rpc.station()==1) result.hasStation1 = true;
160			if(rpc.station()==2) result.hasStation2 = true;
161	konec	1.1	}
162	akalinow	1.6	}
163			else if (detId.subdetId() == MuonSubdetId::DT) {
164	konec	1.1	DTphDigiSpec digi(rawId, is->second);
165	akalinow	1.6	mType = GoldenPattern::POSDT;
166			cit = PattCore.find(mType);
167	akalinow	1.7	if(cit==PattCore.cend()) return result;
168	akalinow	1.6	idm = cit->second.find(rawId);
169			if (idm != cit->second.cend() ) {
170	akalinow	1.7	float f = whereInDistribution(digi.phi(), idm->second);
171	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
172	konec	1.3	DTChamberId dt(rawId);
173			if(dt.station()==1) result.hasStation1 = true;
174			if(dt.station()==2) result.hasStation2 = true;
175	konec	1.1	}
176	akalinow	1.6	mType = GoldenPattern::BENDT;
177			cit = PattCore.find(mType);
178	akalinow	1.7	if(cit==PattCore.cend()) return result;
179	akalinow	1.6	idm = cit->second.find(rawId);
180			if (idm != cit->second.cend() ) {
181	akalinow	1.7	float f = whereInDistribution(digi.phiB(), idm->second);
182	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
183	konec	1.1	}
184	akalinow	1.6	}
185			else if (detId.subdetId() == MuonSubdetId::CSC) {
186	konec	1.1	CSCDigiSpec digi(rawId, is->second);
187	akalinow	1.6	mType = GoldenPattern::POSCSC;
188			cit = PattCore.find(mType);
189	akalinow	1.7	if(cit==PattCore.cend()) return result;
190	akalinow	1.6	auto idm = cit->second.find(rawId);
191			if (idm != cit->second.cend() ) {
192	akalinow	1.7	float f = whereInDistribution(digi.strip(), idm->second);
193	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
194	konec	1.3	CSCDetId csc(rawId);
195			if (csc.station()==1) result.hasStation1 = true;
196			if (csc.station()==2) result.hasStation1 = true;
197	konec	1.1	}
198	akalinow	1.6	mType = GoldenPattern::BENCSC;
199			cit = PattCore.find(mType);
200	akalinow	1.7	if(cit==PattCore.cend()) return result;
201	akalinow	1.6	idm = cit->second.find(rawId);
202			if (idm != cit->second.cend() ) {
203	akalinow	1.7	float f = whereInDistribution(digi.pattern(), idm->second);
204	akalinow	1.6	result.myResults[mType].push_back(std::make_pair(rawId, f));
205	konec	1.1	}
206			}
207			}
208	akalinow	1.6
209	konec	1.1	return result;
210			}
211
212	akalinow	1.7	float GoldenPattern::whereInDistribution( int obj, const GoldenPattern::MFreq & m) const
213	konec	1.1	{
214	akalinow	1.6
215	akalinow	1.7	float sum_before = 0;
216			float sum_after = 0;
217			float sum_obj = 0;
218	akalinow	1.6	for (MFreq::const_iterator im = m.begin(); im != m.end(); ++im) {
219	konec	1.1	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	akalinow	1.7	float sum = std::max((float)1.,sum_before+sum_after+sum_obj );
224	akalinow	1.6	//return sum_obj/sum; //AK
225	konec	1.1	return (sum_before+sum_obj/2.)/sum;
226			}
227
228	akalinow	1.7	bool GoldenPattern::purge(){
229	akalinow	1.5
230	akalinow	1.6	bool remove = false;
231			int pos;
232			unsigned int bef2, bef1, aft1, aft2, aft3;
233	akalinow	1.7	for (auto isf=PattCore.begin();isf!=PattCore.end();){
234	akalinow	1.6	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	konec	1.1	}
251	akalinow	1.7	if (isf->second.size()==0) PattCore.erase(isf++); else ++isf;
252	konec	1.1	}
253	akalinow	1.7	///Usefull pattern has at least 4 measurements and has a RPC measurement
254			return PattCore.find(POSRPC)!=PattCore.end() && PattCore.size()>4;
255	akalinow	1.6	}
256	konec	1.1
257			std::ostream & operator << (std::ostream &out, const GoldenPattern & o) {
258	akalinow	1.6
259			out <<"GoldenPattern"<< o.theKey <<std::endl;
260
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	akalinow	1.7	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	akalinow	1.6	for (auto imf = idf->second.cbegin(); imf != idf->second.cend();++imf)
280			{ out << imf->first<<":"<<imf->second<<", "; }
281	konec	1.1	out << std::endl;
282			}
283	akalinow	1.6	}
284			return out;
285			}
286