L1RpcTriggerAnalysis/src/PatternManager.cc

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

#include "UserCode/L1RpcTriggerAnalysis/interface/EventObj.h"
#include "UserCode/L1RpcTriggerAnalysis/interface/TrackObj.h"
#include "UserCode/L1RpcTriggerAnalysis/interface/HitSpecObj.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"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"

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

#include "TFile.h"
#include "TTree.h"

namespace {
  typedef struct {
    UInt_t key_det;
    unsigned int key_pt;
    unsigned int key_phi;
    int key_ch;
    unsigned int pat_Case;
    UInt_t patDet;
    int posOrBend;
    unsigned int freq; 
  } ENTRY;
}

PatternManager::PatternManager(const edm::ParameterSet &cfg) 
  : theConfig(cfg), theEvForPatCounter(0), theEvUsePatCounter(0) 

{}

PatternManager::~PatternManager()
{
  std::cout <<" Events checked for pattenrs:    " << theEvForPatCounter << std::endl;
  std::cout <<" Events used    for pattenrs:    " << theEvUsePatCounter << std::endl;
  std::cout <<" Size of GoldenPatterns:         " << theGPs.size() << std::endl;
  if (theConfig.getUntrackedParameter<bool>("dump",false)) {
  for (std::map< GoldenPattern::Key, GoldenPattern>::const_iterator 
      im = theGPs.begin(); im!= theGPs.end(); ++im) {
    std::cout <<" GP: "<< (*im).second << std::endl;
  }
  }

}

void PatternManager::run(const EventObj* ev, const TrackObj * simu, const HitSpecObj * hitSpec,  const VDigiSpec & vDigi)
{

  if (theConfig.exists("patternInpFile") && 
      theConfig.exists("patternOutFile")) return;

  if (!hitSpec) return;
  if (hitSpec->rawId() == 0 ) return;

  int stripRef = -1;
  for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {   
    if(is->first!=hitSpec->rawId()) continue;
    RPCDigiSpec digi(is->first, is->second);
    stripRef = digi.halfStrip();
  } 

  ///Use only events with a digi in reference station
  if(stripRef<0) return;

  //float phiref = hitSpec->position().phi();
  float phiref = stripRef;
  float ptref  = simu->pt(); 
  int    chargeref = simu->charge();
  unsigned int detref =  hitSpec->rawId();

  GoldenPattern::Key key( detref, ptref, chargeref, phiref);

  Pattern pattern;
  theEvForPatCounter++;
  for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {
    DetId detId( is->first);    
    bool isOK = pattern.add(*is); 
    if (!isOK) return;
  }
  if (pattern.size() == 0) return; 
  theEvUsePatCounter++;


  if (theGPs.find(key)==theGPs.end()) theGPs[key]=GoldenPattern(key); 
  theGPs[key].add(pattern);
}

L1Obj PatternManager::check(const EventObj* ev, const TrackObj * simu, const HitSpecObj * hitSpec,  const VDigiSpec & vDigi)
{
  L1Obj candidate;
  if (theConfig.exists("patternInpFile") && 
      theConfig.exists("patternOutFile")) return candidate;

  if (!hitSpec) return candidate;
  if (hitSpec->rawId() == 0 ) return candidate;
  //float phiref = hitSpec->position().phi();
  float ptref  = simu->pt();
  int    chargeref = simu->charge();
  unsigned int detref =  hitSpec->rawId();

  VDigiSpec::const_iterator is=vDigi.end();
  for (is = vDigi.begin(); is!=vDigi.end(); ++is) {   
    if(is->first==hitSpec->rawId()) break;
  }
  if(is==vDigi.end()){return candidate;} //Require a digi in reference station. Creates inefficiency!
  RPCDigiSpec digiRef(is->first, is->second);
  int stripRef = digiRef.halfStrip();
  
  std::vector<Pattern> vpattern(1);
  theEvForPatCounter++;
  static bool skipRpcData   = theConfig.getUntrackedParameter<bool>("skipRpcData",  false);
  static bool skipDtCscData = theConfig.getUntrackedParameter<bool>("skipDtCscData",false);
  for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {
    DetId detId( is->first);
    if (skipRpcData   && detId.subdetId()==MuonSubdetId::RPC) continue;
    if (skipDtCscData && (detId.subdetId()==MuonSubdetId::DT || detId.subdetId()==MuonSubdetId::CSC) ) continue;
    if(detId.subdetId()==MuonSubdetId::DT && DTphDigiSpec(is->first, is->second).bxNum()!=0) continue;
    Pattern::add(vpattern,*is);
    if (vpattern.size() > 100) break;
  }
  if (vpattern[0].size() == 0) return candidate;

  GoldenPattern::Key thisKey(detref, ptref, chargeref,stripRef);
  //std::cout << thisKey << std::endl;
  /*
  RPCDetId aId(637599786);
  std::cout<<aId<<std::endl;
  RPCDetId bId(637599914);
  std::cout<<bId<<std::endl;
  */
  GoldenPattern::Result bestMatching;
  GoldenPattern::Key    bestKey;
  for (auto ip=vpattern.cbegin(); ip!= vpattern.cend();++ip) {
    Pattern pattern;
    int phiRotation = -1;
    ///Find rotation degree to match the base set of patterns. This still can be optimised.
    for (std::map< GoldenPattern::Key, GoldenPattern>::iterator igps = theGPs.begin(); igps != theGPs.end(); ++igps) {
      ///Rotate the hit pattern to math the golden pattern      
      for(int iStep=0;iStep<5;++iStep){
        uint32_t rawIdRotated = Pattern::rotateDetId(digiRef.rawId(),iStep);
        if(rawIdRotated==igps->first.theDet && igps->first.thePhiCode==digiRef.halfStrip()){
          pattern = ip->getRotated(iStep);
          phiRotation = iStep;
          break;
        }
      }
      GoldenPattern::Key thisRotatedKey(Pattern::rotateDetId(digiRef.rawId(),phiRotation), ptref, chargeref,stripRef);
      if(!pattern.size()) continue;
      //if(pattern.size()<7) return candidate;
      //if(!(thisRotatedKey==igps->first)) continue;
      
      GoldenPattern & gp = igps->second;
      GoldenPattern::Result result = gp.compare(pattern);  
      /*
      std::cout<<"Step: "<<phiRotation<<std::endl;
      std::cout<<*ip<<std::endl;
      std::cout<<pattern<<std::endl;
      std::cout<<gp<<std::endl;
      std::cout <<"PATT KEY: "<<igps->first<<" "<<result<<std::endl;
      */
      if (bestMatching < result) {
        bestMatching = result;
        bestKey =  igps->first;
        //      std::cout <<" ----"<<" pt: "<< bestKey.ptValue()<<std::endl; 
      } 
      //else std::cout <<std::endl;
    }
  }

  //std::cout <<" ------------------ END COMPARE: " << std::endl;
  //std::cout <<"THIS KEY: "<<thisKey << std::endl;
  //std::cout <<"BEST KEY: "<<bestKey<<" "<< bestMatching<<std::endl;
  //abort();
  if (bestMatching) {
    candidate.pt = bestKey.ptValue();
    candidate.eta = bestKey.etaValue();
    candidate.phi = bestKey.phiValue();
    candidate.charge = bestKey.chargeValue();
    candidate.q   = bestMatching.nMatchedTot();
    candidate.type = L1Obj::OTF;
  }
  return candidate;

}


void PatternManager::beginJob()
{
  if ( !theConfig.exists("patternInpFile") ) return;
  std::string patternInpFileName = theConfig.getParameter<std::string>("patternInpFile");
  TFile patternInpFile( patternInpFileName.c_str());
  TTree * tree = (TTree*) patternInpFile.Get("FlatPatterns");
  static ENTRY entry;
  tree->SetBranchAddress("entry",&entry); 
  Int_t nentries = (Int_t) tree->GetEntries();
  for (Int_t i=0; i<nentries; ++i) {
    tree->GetEntry(i);
    GoldenPattern::Key key;
    key.theDet =     entry.key_det;
    key.thePtCode =  entry.key_pt;
    key.thePhiCode = entry.key_phi;
    key.theCharge =  entry.key_ch;

    GoldenPattern::PosBenCase pat_Case = static_cast<GoldenPattern::PosBenCase>(entry.pat_Case);
    if (theGPs.find(key)==theGPs.end()) theGPs[key]=GoldenPattern(key); 
    theGPs[key].add(pat_Case, entry.patDet, entry.posOrBend, entry.freq);
  }
  delete tree;
  patternInpFile.Close();
}

void PatternManager::endJob(){

  if ( !theConfig.exists("patternOutFile") ) return;
  std::string patternOutFileName = theConfig.getParameter<std::string>("patternOutFile");
  TFile patternOutFile( patternOutFileName.c_str(),"RECREATE"); 
  
  if (theConfig.exists("patternInpFile") && theConfig.exists("patternOutFile")){
    for (auto igps = theGPs.begin(); igps != theGPs.end();) {
      GoldenPattern & gp = igps->second;      
      if(!gp.purge()) {theGPs.erase(igps++); } else { ++igps; } 
    }
  }
  
  static ENTRY entry;
  TTree *tree = new TTree("FlatPatterns","FlatPatterns");
  tree->Branch("entry",&entry,"key_det/i:key_pt/i:key_phi/i:key_ch/I:pat_Case/i:patDet/i:posOrBend/I:freq/i");

  for (std::map< GoldenPattern::Key, GoldenPattern>::const_iterator igps = theGPs.begin(); igps != theGPs.end(); ++igps) {
    const GoldenPattern & gp = igps->second; 
    entry.key_det =  gp.theKey.theDet;
    entry.key_pt =  gp.theKey.thePtCode;
    entry.key_phi =  gp.theKey.thePhiCode;
    entry.key_ch =  gp.theKey.theCharge;
    
    for (auto it=gp.PattCore.cbegin();it!=gp.PattCore.cend();++it){
      entry.pat_Case = it->first; 
      for (GoldenPattern::DetFreq::const_iterator idf = it->second.cbegin(); idf !=it->second.cend(); ++idf) {
        entry.patDet = idf->first;
        for (GoldenPattern::MFreq::const_iterator imf = idf->second.cbegin(); imf != idf->second.cend(); ++imf) {
          entry.posOrBend = imf->first;
          entry.freq = imf->second; 
          tree->Fill();
        }
      }
    }
  } 
  patternOutFile.Write();
  patternOutFile.Close();
//  delete tree;
}

Revision:	1.13
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.12:	+48 -35 lines
Log Message:	*last commit before migration to Git.
#	User	Rev	Content
1	konec	1.1	#include "UserCode/L1RpcTriggerAnalysis/interface/PatternManager.h"
2
3			#include "UserCode/L1RpcTriggerAnalysis/interface/EventObj.h"
4			#include "UserCode/L1RpcTriggerAnalysis/interface/TrackObj.h"
5			#include "UserCode/L1RpcTriggerAnalysis/interface/HitSpecObj.h"
6
7			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
8			#include "UserCode/L1RpcTriggerAnalysis/interface/DTphDigiSpec.h"
9			#include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
10			#include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"
11	konec	1.8	#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
12	konec	1.1
13			#include "UserCode/L1RpcTriggerAnalysis/interface/Pattern.h"
14			#include "UserCode/L1RpcTriggerAnalysis/interface/GoldenPattern.h"
15
16			#include "TFile.h"
17			#include "TTree.h"
18
19			namespace {
20			typedef struct {
21			UInt_t key_det;
22			unsigned int key_pt;
23			unsigned int key_phi;
24			int key_ch;
25			unsigned int pat_Case;
26			UInt_t patDet;
27			int posOrBend;
28			unsigned int freq;
29			} ENTRY;
30			}
31
32			PatternManager::PatternManager(const edm::ParameterSet &cfg)
33			: theConfig(cfg), theEvForPatCounter(0), theEvUsePatCounter(0)
34
35			{}
36
37			PatternManager::~PatternManager()
38			{
39			std::cout <<" Events checked for pattenrs: " << theEvForPatCounter << std::endl;
40	konec	1.5	std::cout <<" Events used for pattenrs: " << theEvUsePatCounter << std::endl;
41	konec	1.1	std::cout <<" Size of GoldenPatterns: " << theGPs.size() << std::endl;
42			if (theConfig.getUntrackedParameter<bool>("dump",false)) {
43			for (std::map< GoldenPattern::Key, GoldenPattern>::const_iterator
44			im = theGPs.begin(); im!= theGPs.end(); ++im) {
45			std::cout <<" GP: "<< (*im).second << std::endl;
46			}
47			}
48
49			}
50
51			void PatternManager::run(const EventObj* ev, const TrackObj * simu, const HitSpecObj * hitSpec, const VDigiSpec & vDigi)
52			{
53	akalinow	1.12
54			if (theConfig.exists("patternInpFile") &&
55			theConfig.exists("patternOutFile")) return;
56
57	konec	1.1	if (!hitSpec) return;
58			if (hitSpec->rawId() == 0 ) return;
59	akalinow	1.12
60			int stripRef = -1;
61			for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {
62			if(is->first!=hitSpec->rawId()) continue;
63			RPCDigiSpec digi(is->first, is->second);
64			stripRef = digi.halfStrip();
65			}
66
67			///Use only events with a digi in reference station
68			if(stripRef<0) return;
69
70			//float phiref = hitSpec->position().phi();
71			float phiref = stripRef;
72			float ptref = simu->pt();
73	konec	1.1	int chargeref = simu->charge();
74			unsigned int detref = hitSpec->rawId();
75
76			GoldenPattern::Key key( detref, ptref, chargeref, phiref);
77
78			Pattern pattern;
79			theEvForPatCounter++;
80	akalinow	1.11	for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {
81	akalinow	1.13	DetId detId( is->first);
82	konec	1.1	bool isOK = pattern.add(*is);
83			if (!isOK) return;
84			}
85			if (pattern.size() == 0) return;
86			theEvUsePatCounter++;
87
88	akalinow	1.13
89	konec	1.1	if (theGPs.find(key)==theGPs.end()) theGPs[key]=GoldenPattern(key);
90			theGPs[key].add(pattern);
91			}
92
93			L1Obj PatternManager::check(const EventObj* ev, const TrackObj * simu, const HitSpecObj * hitSpec, const VDigiSpec & vDigi)
94			{
95			L1Obj candidate;
96	akalinow	1.12	if (theConfig.exists("patternInpFile") &&
97			theConfig.exists("patternOutFile")) return candidate;
98
99	konec	1.1	if (!hitSpec) return candidate;
100			if (hitSpec->rawId() == 0 ) return candidate;
101	akalinow	1.13	//float phiref = hitSpec->position().phi();
102	akalinow	1.12	float ptref = simu->pt();
103	konec	1.1	int chargeref = simu->charge();
104			unsigned int detref = hitSpec->rawId();
105
106	akalinow	1.13	VDigiSpec::const_iterator is=vDigi.end();
107			for (is = vDigi.begin(); is!=vDigi.end(); ++is) {
108			if(is->first==hitSpec->rawId()) break;
109			}
110			if(is==vDigi.end()){return candidate;} //Require a digi in reference station. Creates inefficiency!
111			RPCDigiSpec digiRef(is->first, is->second);
112			int stripRef = digiRef.halfStrip();
113
114	konec	1.6	std::vector<Pattern> vpattern(1);
115	konec	1.1	theEvForPatCounter++;
116	konec	1.8	static bool skipRpcData = theConfig.getUntrackedParameter<bool>("skipRpcData", false);
117			static bool skipDtCscData = theConfig.getUntrackedParameter<bool>("skipDtCscData",false);
118	akalinow	1.11	for (VDigiSpec::const_iterator is= vDigi.begin(); is!=vDigi.end(); ++is) {
119	konec	1.8	DetId detId( is->first);
120			if (skipRpcData && detId.subdetId()==MuonSubdetId::RPC) continue;
121			if (skipDtCscData && (detId.subdetId()==MuonSubdetId::DT \|\| detId.subdetId()==MuonSubdetId::CSC) ) continue;
122	akalinow	1.13	if(detId.subdetId()==MuonSubdetId::DT && DTphDigiSpec(is->first, is->second).bxNum()!=0) continue;
123	konec	1.8	Pattern::add(vpattern,*is);
124	konec	1.9	if (vpattern.size() > 100) break;
125	konec	1.8	}
126	konec	1.6	if (vpattern[0].size() == 0) return candidate;
127	konec	1.1
128	akalinow	1.13	GoldenPattern::Key thisKey(detref, ptref, chargeref,stripRef);
129	akalinow	1.12	//std::cout << thisKey << std::endl;
130	akalinow	1.13	/*
131			RPCDetId aId(637599786);
132			std::cout<<aId<<std::endl;
133			RPCDetId bId(637599914);
134			std::cout<<bId<<std::endl;
135			*/
136	konec	1.1	GoldenPattern::Result bestMatching;
137			GoldenPattern::Key bestKey;
138	akalinow	1.11	for (auto ip=vpattern.cbegin(); ip!= vpattern.cend();++ip) {
139	akalinow	1.13	Pattern pattern;
140			int phiRotation = -1;
141			///Find rotation degree to match the base set of patterns. This still can be optimised.
142	akalinow	1.11	for (std::map< GoldenPattern::Key, GoldenPattern>::iterator igps = theGPs.begin(); igps != theGPs.end(); ++igps) {
143	akalinow	1.13	///Rotate the hit pattern to math the golden pattern
144			for(int iStep=0;iStep<5;++iStep){
145			uint32_t rawIdRotated = Pattern::rotateDetId(digiRef.rawId(),iStep);
146			if(rawIdRotated==igps->first.theDet && igps->first.thePhiCode==digiRef.halfStrip()){
147			pattern = ip->getRotated(iStep);
148			phiRotation = iStep;
149			break;
150			}
151			}
152			GoldenPattern::Key thisRotatedKey(Pattern::rotateDetId(digiRef.rawId(),phiRotation), ptref, chargeref,stripRef);
153			if(!pattern.size()) continue;
154			//if(pattern.size()<7) return candidate;
155			//if(!(thisRotatedKey==igps->first)) continue;
156
157	akalinow	1.11	GoldenPattern & gp = igps->second;
158	akalinow	1.13	GoldenPattern::Result result = gp.compare(pattern);
159			/*
160			std::cout<<"Step: "<<phiRotation<<std::endl;
161			std::cout<<*ip<<std::endl;
162			std::cout<<pattern<<std::endl;
163			std::cout<<gp<<std::endl;
164	akalinow	1.12	std::cout <<"PATT KEY: "<<igps->first<<" "<<result<<std::endl;
165	akalinow	1.13	*/
166	akalinow	1.11	if (bestMatching < result) {
167			bestMatching = result;
168			bestKey = igps->first;
169			// std::cout <<" ----"<<" pt: "<< bestKey.ptValue()<<std::endl;
170			}
171			//else std::cout <<std::endl;
172			}
173	konec	1.6	}
174	konec	1.5
175	akalinow	1.13	//std::cout <<" ------------------ END COMPARE: " << std::endl;
176			//std::cout <<"THIS KEY: "<<thisKey << std::endl;
177			//std::cout <<"BEST KEY: "<<bestKey<<" "<< bestMatching<<std::endl;
178			//abort();
179	konec	1.1	if (bestMatching) {
180			candidate.pt = bestKey.ptValue();
181	akalinow	1.11	candidate.eta = bestKey.etaValue();
182	konec	1.1	candidate.phi = bestKey.phiValue();
183	akalinow	1.13	candidate.charge = bestKey.chargeValue();
184	konec	1.1	candidate.q = bestMatching.nMatchedTot();
185			candidate.type = L1Obj::OTF;
186			}
187			return candidate;
188
189			}
190
191
192			void PatternManager::beginJob()
193			{
194			if ( !theConfig.exists("patternInpFile") ) return;
195			std::string patternInpFileName = theConfig.getParameter<std::string>("patternInpFile");
196			TFile patternInpFile( patternInpFileName.c_str());
197			TTree * tree = (TTree*) patternInpFile.Get("FlatPatterns");
198			static ENTRY entry;
199			tree->SetBranchAddress("entry",&entry);
200			Int_t nentries = (Int_t) tree->GetEntries();
201	akalinow	1.11	for (Int_t i=0; i<nentries; ++i) {
202	konec	1.1	tree->GetEntry(i);
203			GoldenPattern::Key key;
204			key.theDet = entry.key_det;
205			key.thePtCode = entry.key_pt;
206			key.thePhiCode = entry.key_phi;
207			key.theCharge = entry.key_ch;
208	akalinow	1.12
209	konec	1.1	GoldenPattern::PosBenCase pat_Case = static_cast<GoldenPattern::PosBenCase>(entry.pat_Case);
210			if (theGPs.find(key)==theGPs.end()) theGPs[key]=GoldenPattern(key);
211			theGPs[key].add(pat_Case, entry.patDet, entry.posOrBend, entry.freq);
212			}
213			delete tree;
214			patternInpFile.Close();
215			}
216
217	akalinow	1.11	void PatternManager::endJob(){
218	konec	1.1
219			if ( !theConfig.exists("patternOutFile") ) return;
220			std::string patternOutFileName = theConfig.getParameter<std::string>("patternOutFile");
221			TFile patternOutFile( patternOutFileName.c_str(),"RECREATE");
222	akalinow	1.12
223			if (theConfig.exists("patternInpFile") && theConfig.exists("patternOutFile")){
224			for (auto igps = theGPs.begin(); igps != theGPs.end();) {
225			GoldenPattern & gp = igps->second;
226			if(!gp.purge()) {theGPs.erase(igps++); } else { ++igps; }
227			}
228	akalinow	1.11	}
229	akalinow	1.12
230	konec	1.1	static ENTRY entry;
231			TTree *tree = new TTree("FlatPatterns","FlatPatterns");
232			tree->Branch("entry",&entry,"key_det/i:key_pt/i:key_phi/i:key_ch/I:pat_Case/i:patDet/i:posOrBend/I:freq/i");
233
234	akalinow	1.11	for (std::map< GoldenPattern::Key, GoldenPattern>::const_iterator igps = theGPs.begin(); igps != theGPs.end(); ++igps) {
235	konec	1.1	const GoldenPattern & gp = igps->second;
236			entry.key_det = gp.theKey.theDet;
237			entry.key_pt = gp.theKey.thePtCode;
238			entry.key_phi = gp.theKey.thePhiCode;
239			entry.key_ch = gp.theKey.theCharge;
240	akalinow	1.11
241			for (auto it=gp.PattCore.cbegin();it!=gp.PattCore.cend();++it){
242			entry.pat_Case = it->first;
243			for (GoldenPattern::DetFreq::const_iterator idf = it->second.cbegin(); idf !=it->second.cend(); ++idf) {
244	konec	1.1	entry.patDet = idf->first;
245	akalinow	1.11	for (GoldenPattern::MFreq::const_iterator imf = idf->second.cbegin(); imf != idf->second.cend(); ++imf) {
246	konec	1.1	entry.posOrBend = imf->first;
247			entry.freq = imf->second;
248			tree->Fill();
249			}
250			}
251			}
252			}
253			patternOutFile.Write();
254			patternOutFile.Close();
255			// delete tree;
256			}
257