ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/L1RpcTriggerAnalysis/src/AnaDigiSpec.cc
Revision: 1.3
Committed: Wed May 29 13:43:33 2013 UTC (11 years, 11 months ago) by konec
Content type: text/plain
Branch: MAIN
CVS Tags: Artur_11_07_2013_B, Artur_11_07_2013_A, Artur_11_07_2013, Artur_28_06_2013, HEAD
Changes since 1.2: +1 -1 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 konec 1.1 #include "UserCode/L1RpcTriggerAnalysis/interface/AnaDigiSpec.h"
2    
3     #include "UserCode/L1RpcTriggerAnalysis/interface/EventObj.h"
4     #include "UserCode/L1RpcTriggerAnalysis/interface/TrackObj.h"
5     #include "UserCode/L1RpcTriggerAnalysis/interface/HitSpecObj.h"
6     #include "UserCode/L1RpcTriggerAnalysis/interface/Utilities.h"
7    
8     #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
9     #include "DataFormats/MuonDetId/interface/RPCDetId.h"
10     #include "DataFormats/MuonDetId/interface/CSCDetId.h"
11     #include "DataFormats/MuonDetId/interface/DTChamberId.h"
12    
13     #include "UserCode/L1RpcTriggerAnalysis/interface/RPCDetIdUtil.h"
14     #include "UserCode/L1RpcTriggerAnalysis/interface/DTphDigiSpec.h"
15     #include "UserCode/L1RpcTriggerAnalysis/interface/CSCDigiSpec.h"
16     #include "UserCode/L1RpcTriggerAnalysis/interface/RPCDigiSpec.h"
17    
18     #include "TObjArray.h"
19     #include "TH1D.h"
20     #include "TH2D.h"
21     #include "TProfile.h"
22     #include "TGraphErrors.h"
23    
24     #include <cmath>
25     #include <algorithm>
26    
27     template <class T> T sqr( T t) {return t*t;}
28     namespace {
29     TH2D *hDigiSpec_PosCorr_DT2, *hDigiSpec_PosCorr_CSC2, *hDigiSpec_PosCorr_RPC2;
30    
31     }
32     void AnaDigiSpec::init(TObjArray& histos)
33     {
34     hDigiSpec_PosCorr_DT2 = new TH2D("hDigiSpec_PosCorr_DT2", "hDigiSpec_PosCorr_DT2", 100,0.8,1.2, 100, -2000.,2000.);
35     histos.Add(hDigiSpec_PosCorr_DT2);
36     hDigiSpec_PosCorr_RPC2 = new TH2D("hDigiSpec_PosCorr_RPC2","hDigiSpec_PosCorr_RPC2",100,0.8,1.2, 191, 2.,193.);
37     histos.Add(hDigiSpec_PosCorr_RPC2);
38     hDigiSpec_PosCorr_CSC2 = new TH2D("hDigiSpec_PosCorr_CSC2","hDigiSpec_PosCorr_CSC2",100,0.8,1.2, 160, 0.,159.);
39     histos.Add(hDigiSpec_PosCorr_CSC2);
40     }
41    
42     void AnaDigiSpec::run(const EventObj* ev, const TrackObj * simu, const HitSpecObj * hitSpec, const VDigiSpec & vdigSpec)
43     {
44     //
45     // position correlation at station 2 wrt RPC SimHit
46     // skip events with not RPC simhit in proper place
47     //
48     if (!hitSpec) return;
49     uint32_t hitRawId = hitSpec->rawId();
50     if (!hitRawId) return;
51     // if (hitSpec->position().phi() < 0.6 && hitSpec->position().phi() > 1.20) return;
52     if (hitRawId != 637602109 && hitRawId != 637634877 &&
53     // hitRawId != 637620266 && hitRawId != 637653034 &&
54     hitRawId != 637599914 && hitRawId != 637632682 ) return;
55     // if (hitRawId != 637599914) return;
56     // if (hitRawId != 637632682 ) return;
57     RPCDetIdUtil idHit( hitRawId );
58 konec 1.2 bool barrelRefHit = idHit.isBarrel();
59 konec 1.1
60     //
61     // correlation betweem sim hit position@RPC and digi
62     //
63     for (VDigiSpec::const_iterator is= vdigSpec.begin(); is!=vdigSpec.end(); is++) {
64    
65     DetId detId(is->first);
66     if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
67     switch (detId.subdetId()) {
68     case MuonSubdetId::RPC: {
69     RPCDigiSpec digi(is->first, is->second);
70     if (digi.rawId() != hitRawId) continue;
71     // std::vector<uint32_t> allDetsInDigis;
72     // for (VDigiSpec::const_iterator is= vdigSpec.begin(); is!=vdigSpec.end(); is++) allDetsInDigis.push_back(is->first);
73     // unsigned int ndet = std::count( allDetsInDigis.begin(), allDetsInDigis.end(), hitRawId);
74     // if (digi.fromStrip() == digi.toStrip() && ndet == 1 )
75     hDigiSpec_PosCorr_RPC2->Fill(hitSpec->position().phi(), digi.halfStrip());
76     break;
77     }
78     case MuonSubdetId::DT: {
79     // if (!idHit.isBarrel()) break;
80     DTphDigiSpec digi(is->first, is->second);
81     DTChamberId dt(digi.rawId());
82     if (dt.station() != 2) break;
83     hDigiSpec_PosCorr_DT2->Fill(hitSpec->position().phi(), digi.phi());
84     break;
85     }
86     case MuonSubdetId::CSC: {
87     // if (idHit.isBarrel()) break;
88     CSCDigiSpec digi(is->first, is->second);
89     CSCDetId csc(digi.rawId());
90     // std::cout << csc << std::endl;
91     if (csc.station() != 2) break;
92     hDigiSpec_PosCorr_CSC2->Fill(hitSpec->position().phi(), digi.strip());
93     break;
94     }
95     default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
96     };
97     }
98    
99     //
100     // for given fixed hit postion in RPC check distributions
101     //
102 konec 1.2 bool precisePos = ( fabs(hitSpec->position().phi()-1.025) < 0.001);
103 konec 1.1 // bool precisePos = (
104     // fabs(hitSpec->position().phi()-0.99) < 0.001
105     // || fabs(hitSpec->position().phi()-0.89) < 0.001 );
106     for (VDigiSpec::const_iterator is= vdigSpec.begin(); is!=vdigSpec.end(); is++) {
107    
108     DetId detId(is->first);
109     if (detId.det() != DetId::Muon) std::cout << "PROBLEM ;;;"<<std::endl;
110     switch (detId.subdetId()) {
111    
112     case MuonSubdetId::RPC: {
113     RPCDigiSpec digi(is->first, is->second);
114     RPCDetId rpcId(digi.rawId());
115     RPCDetIdUtil util(digi.rawId());
116     if (util.isBarrel()) {
117 konec 1.2 if (barrelRefHit && (theRpcBPosB.find(digi.rawId())==theRpcBPosB.end()) ){
118 konec 1.1 std::stringstream str;
119 konec 1.2 str<<"hDigiSpec_RpcBVsPt_PosB_"<<digi.rawId();
120 konec 1.1 std::string hName = str.str();
121     str<<"_la"<<util.layer()<<"sc"<<rpcId.sector()<<"ro"<<rpcId.roll();
122     std::string hTitle = str.str();
123 konec 1.2 theRpcBPosB[digi.rawId()] = new TH2D(hName.c_str(), hTitle.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 191,2.,193.);
124 konec 1.1 }
125 konec 1.2 if (!barrelRefHit && (theRpcBPosE.find(digi.rawId())==theRpcBPosE.end()) ) {
126     std::stringstream str;
127     str<<"hDigiSpec_RpcBVsPt_PosE_"<<digi.rawId();
128     std::string hName = str.str();
129     str<<"_la"<<util.layer()<<"sc"<<rpcId.sector()<<"ro"<<rpcId.roll();
130     std::string hTitle = str.str();
131     theRpcBPosE[digi.rawId()] = new TH2D(hName.c_str(), hTitle.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 191,2.,193.);
132     }
133     if (precisePos && barrelRefHit) theRpcBPosB[digi.rawId()]->Fill(simu->pt(), digi.halfStrip());
134     if (precisePos && !barrelRefHit) theRpcBPosE[digi.rawId()]->Fill(simu->pt(), digi.halfStrip());
135 konec 1.1 } else {
136 konec 1.2 if (barrelRefHit && (theRpcEPosB.find(digi.rawId())==theRpcEPosB.end()) ) {
137 konec 1.1 std::stringstream str;
138 konec 1.2 str<<"hDigiSpec_RpcEVsPt_PosB_"<<digi.rawId();
139 konec 1.1 std::string hName = str.str();
140 konec 1.3 str<<"_la"<<util.layer()<<"ch"<<(rpcId.sector()-1)*6+rpcId.subsector()<<"rn"<<rpcId.ring()<<"ro"<<rpcId.roll();
141 konec 1.1 std::string hTitle = str.str();
142 konec 1.2 theRpcEPosB[digi.rawId()] = new TH2D(hName.c_str(), hTitle.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 191,2.,193.);
143 konec 1.1 }
144 konec 1.2 if (!barrelRefHit && (theRpcEPosE.find(digi.rawId())==theRpcEPosE.end()) ) {
145     std::stringstream str;
146     str<<"hDigiSpec_RpcEVsPt_PosE_"<<digi.rawId();
147     std::string hName = str.str();
148     str<<"_la"<<util.layer()<<"ch"<<(rpcId.sector()-1)*6+rpcId.subsector()<<"rn"<<rpcId.ring()<<"ro"<<rpcId.roll();
149     std::string hTitle = str.str();
150     theRpcEPosE[digi.rawId()] = new TH2D(hName.c_str(), hTitle.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 191,2.,193.);
151     }
152     if (precisePos && barrelRefHit) theRpcEPosB[digi.rawId()]->Fill(simu->pt(), digi.halfStrip());
153     if (precisePos && !barrelRefHit) theRpcEPosE[digi.rawId()]->Fill(simu->pt(), digi.halfStrip());
154 konec 1.1 }
155     break;
156     }
157    
158     case MuonSubdetId::DT: {
159     DTphDigiSpec digi(is->first, is->second);
160     DTChamberId dt(digi.rawId());
161 konec 1.2 if (barrelRefHit && (theDtPosB.find(digi.rawId())==theDtPosB.end()) ) {
162     std::stringstream str1, str2;
163     str1<<"hDigiSpec_DtVsPt_PosB_"<<digi.rawId();
164     str2<<"hDigiSpec_DtVsPt_BendB_"<<digi.rawId();
165     std::string hName1 = str1.str();
166     std::string hName2 = str2.str();
167     str1<<"_s"<<dt.station()<<"sc"<<dt.sector();
168     str2<<"_s"<<dt.station()<<"sc"<<dt.sector();
169     std::string hTitle1 = str1.str();
170     std::string hTitle2 = str2.str();
171     theDtPosB[digi.rawId()] = new TH2D(hName1.c_str(), hTitle1.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 100,-2048.,2048.);
172     theDtBendB[digi.rawId()] = new TH2D(hName2.c_str(), hTitle2.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 1100,-550.,550.);
173     }
174     if (!barrelRefHit && (theDtPosE.find(digi.rawId())==theDtPosE.end()) ) {
175 konec 1.1 std::stringstream str1, str2;
176 konec 1.2 str1<<"hDigiSpec_DtVsPt_PosE_"<<digi.rawId();
177     str2<<"hDigiSpec_DtVsPt_BendE_"<<digi.rawId();
178 konec 1.1 std::string hName1 = str1.str();
179     std::string hName2 = str2.str();
180     str1<<"_s"<<dt.station()<<"sc"<<dt.sector();
181     str2<<"_s"<<dt.station()<<"sc"<<dt.sector();
182     std::string hTitle1 = str1.str();
183     std::string hTitle2 = str2.str();
184 konec 1.2 theDtPosE[digi.rawId()] = new TH2D(hName1.c_str(), hTitle1.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 100,-2048.,2048.);
185     theDtBendE[digi.rawId()] = new TH2D(hName2.c_str(), hTitle2.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 1100,-550.,550.);
186 konec 1.1 }
187 konec 1.2 if (precisePos && barrelRefHit) theDtPosB[digi.rawId()]->Fill(simu->pt(), digi.phi());
188     if (precisePos && barrelRefHit) theDtBendB[digi.rawId()]->Fill(simu->pt(), digi.phiB());
189     if (precisePos && !barrelRefHit) theDtPosE[digi.rawId()]->Fill(simu->pt(), digi.phi());
190     if (precisePos && !barrelRefHit) theDtBendE[digi.rawId()]->Fill(simu->pt(), digi.phiB());
191 konec 1.1 break;
192     }
193    
194     case MuonSubdetId::CSC: {
195     CSCDigiSpec digi(is->first, is->second);
196     CSCDetId csc(digi.rawId());
197 konec 1.2 if ( barrelRefHit && (theCscPosB.find(digi.rawId())==theCscPosB.end()) ) {
198     std::stringstream str1, str2;
199     str1<<"hDigiSpec_CscVsPt_PosB_"<<digi.rawId();
200     str2<<"hDigiSpec_CscVsPt_BendB_"<<digi.rawId();
201     std::string hName1 = str1.str();
202     std::string hName2 = str2.str();
203     str1<<"_s"<<csc.station()<<"ch"<<csc.chamber()<<"rn"<<csc.ring();
204     str2<<"_s"<<csc.station()<<"ch"<<csc.chamber()<<"rn"<<csc.ring();
205     std::string hTitle1 = str1.str();
206     std::string hTitle2 = str2.str();
207     theCscPosB[digi.rawId()] = new TH2D(hName1.c_str(), hTitle1.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 100,0.,160.);
208     theCscBendB[digi.rawId()] = new TH2D(hName2.c_str(), hTitle2.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 11,0.5,11.5);
209     }
210     if (!barrelRefHit && (theCscPosE.find(digi.rawId())==theCscPosE.end()) ) {
211 konec 1.1 std::stringstream str1, str2;
212 konec 1.2 str1<<"hDigiSpec_CscVsPt_PosE_"<<digi.rawId();
213     str2<<"hDigiSpec_CscVsPt_BendE_"<<digi.rawId();
214 konec 1.1 std::string hName1 = str1.str();
215     std::string hName2 = str2.str();
216 konec 1.2 str1<<"_s"<<csc.station()<<"ch"<<csc.chamber()<<"rn"<<csc.ring();
217     str2<<"_s"<<csc.station()<<"ch"<<csc.chamber()<<"rn"<<csc.ring();
218 konec 1.1 std::string hTitle1 = str1.str();
219     std::string hTitle2 = str2.str();
220 konec 1.2 theCscPosE[digi.rawId()] = new TH2D(hName1.c_str(), hTitle1.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 100,0.,160.);
221     theCscBendE[digi.rawId()] = new TH2D(hName2.c_str(), hTitle2.c_str(), L1PtScale::nPtBins, L1PtScale::ptBins, 11,0.5,11.5);
222 konec 1.1 }
223 konec 1.2 if ( barrelRefHit && precisePos) theCscPosB[digi.rawId()]->Fill(simu->pt(), digi.strip());
224     if ( barrelRefHit && precisePos) theCscBendB[digi.rawId()]->Fill(simu->pt(), digi.pattern());
225     if (!barrelRefHit && precisePos) theCscPosE[digi.rawId()]->Fill(simu->pt(), digi.strip());
226     if (!barrelRefHit && precisePos) theCscBendE[digi.rawId()]->Fill(simu->pt(), digi.pattern());
227 konec 1.1 break;
228     }
229    
230     default: {std::cout <<" Unexpeced sebdet case, id ="<<is->first <<std::endl; return; }
231     };
232     }
233     }
234    
235     void AnaDigiSpec::resume(TObjArray& histos)
236     {
237 konec 1.2 std::cout <<" AnaDigiSpec, size of HMaps (barrel ref hit): CscPosB="<<theCscPosB.size()
238     <<", CscBendB="<<theCscBendB.size()
239     <<", DtPosB="<<theDtPosB.size()
240     <<", DtBendB="<<theDtBendB.size()
241     <<", RpcBPosB="<<theRpcBPosB.size()
242     <<", RpcEPosB="<<theRpcEPosB.size()<<std::endl;
243     std::cout <<" AnaDigiSpec, size of HMaps (endcal ref hit): CscPosE="<<theCscPosE.size()
244     <<", CscBendE="<<theCscBendE.size()
245     <<", DtPosE="<<theDtPosE.size()
246     <<", DtBendE="<<theDtBendE.size()
247     <<", RpcBPosE="<<theRpcBPosE.size()
248     <<", RpcEPosE="<<theRpcEPosE.size()<<std::endl;
249     VHisto aCscPos = topN(theCscPosB, 9);
250     VHisto aRpcBPos = topN(theRpcBPosB, 9);
251     VHisto aRpcEPos = topN(theRpcEPosB, 9);
252     VHisto aDtPos = topN(theDtPosB, 9);
253     VHisto aCscBend = topN(theCscBendB, 9);
254     VHisto aDtBend = topN(theDtBendB, 9);
255 konec 1.1 typedef VHisto::const_iterator IV;
256     for (IV iv = aCscPos.begin(); iv != aCscPos.end(); ++iv) histos.Add(*iv);
257     for (IV iv = aRpcBPos.begin(); iv != aRpcBPos.end(); ++iv) histos.Add(*iv);
258     for (IV iv = aRpcEPos.begin(); iv != aRpcEPos.end(); ++iv) histos.Add(*iv);
259     for (IV iv = aDtPos.begin(); iv != aDtPos.end(); ++iv) histos.Add(*iv);
260     for (IV iv = aDtBend.begin(); iv != aDtBend.end(); ++iv) histos.Add(*iv);
261     for (IV iv = aCscBend.begin(); iv != aCscBend.end(); ++iv) histos.Add(*iv);
262    
263 konec 1.2 aCscPos = topN(theCscPosE, 9);
264     aRpcBPos = topN(theRpcBPosE, 9);
265     aRpcEPos = topN(theRpcEPosE, 9);
266     aDtPos = topN(theDtPosE, 9);
267     aCscBend = topN(theCscBendE, 9);
268     aDtBend = topN(theDtBendE, 9);
269     for (IV iv = aCscPos.begin(); iv != aCscPos.end(); ++iv) histos.Add(*iv);
270     for (IV iv = aRpcBPos.begin(); iv != aRpcBPos.end(); ++iv) histos.Add(*iv);
271     for (IV iv = aRpcEPos.begin(); iv != aRpcEPos.end(); ++iv) histos.Add(*iv);
272     for (IV iv = aDtPos.begin(); iv != aDtPos.end(); ++iv) histos.Add(*iv);
273     for (IV iv = aDtBend.begin(); iv != aDtBend.end(); ++iv) histos.Add(*iv);
274     for (IV iv = aCscBend.begin(); iv != aCscBend.end(); ++iv) histos.Add(*iv);
275    
276 konec 1.1 }
277    
278     AnaDigiSpec::VHisto AnaDigiSpec::topN(const HMap & hmap, unsigned int nElements) const
279     {
280     VHisto result;
281     typedef std::map<unsigned int, VHisto> StatMap;
282     StatMap statMap;
283    
284     for (HMap::const_iterator im = hmap.begin(); im != hmap.end(); ++im) {
285     statMap[im->second->Integral()].push_back(im->second);
286     }
287    
288     unsigned int nCurrElements = 0;
289     for (StatMap::const_reverse_iterator im=statMap.rbegin(); im != statMap.rend(); im++) {
290     for (VHisto::const_iterator ih=im->second.begin(); ih != im->second.end(); ih++) {
291     std::string name( (*ih)->GetName() );
292     std::size_t pos = name.find_last_of('_');
293     std::string newName = name.substr(0,pos)+"_"+std::to_string(++nCurrElements);
294     // std::cout<<" name: "<< name <<" NEW NAME: " << newName << std::endl;
295     result.push_back( static_cast<TH2D*>( (*ih)->Clone(newName.c_str()) ) );
296     }
297     if (result.size() >= nElements) return result;
298     }
299     return result;
300     }