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.;
  //AK theValue = ( nTot > 4) ? -log(fract) : 9999.0;
}

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