145 cout <<
"[DTEfficiencyTask] Analyze #Run: " <<
event.id().run() <<
" #Event: " <<
event.id().event() << endl;
161 for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
166 const vector<const DTSuperLayer*>& SLayers =
chamber->superLayers();
167 map<DTWireId, int> wireAnd1DRecHits;
168 for (vector<const DTSuperLayer*>::const_iterator superlayer = SLayers.begin(); superlayer != SLayers.end();
173 wireAnd1DRecHits[(*rechit).wireId()] = (*rechit).wireId().wire();
181 cout <<
" Chamber: " << *chamberId <<
" has " << nsegm <<
" 4D segments" << endl;
186 cout <<
" == RecSegment dimension: " << (*segment4D).dimension() << endl;
190 if ((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
192 cout <<
"[DTEfficiencyTask]***Warning: RecSegment dimension is not 4 but " << (*segment4D).dimension()
193 <<
", skipping!" << endl;
195 }
else if ((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
197 cout <<
"[DTEfficiencyTask]***Warning: RecSegment dimension is not 2 but " << (*segment4D).dimension()
198 <<
", skipping!" << endl;
202 vector<DTRecHit1D> recHits1D;
210 if (phiRecHits.size() < 7 || phiRecHits.size() > 8) {
212 cout <<
"[DTEfficiencyTask] Phi segments has: " << phiRecHits.size() <<
" hits, skipping"
216 copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D));
218 if ((*segment4D).dimension() == 4) {
220 zSeg = (*segment4D).zSegment();
222 if (zRecHits.size() < 3 || zRecHits.size() > 4) {
224 cout <<
"[DTEfficiencyTask] Theta segments has: " << zRecHits.size() <<
" hits, skipping"
228 copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D));
232 vector<DTWireId> wireMap;
233 for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D.begin(); recHit1D != recHits1D.end(); recHit1D++) {
234 wireMap.push_back((*recHit1D).wireId());
237 bool hitsOnSameLayer =
false;
238 for (vector<DTWireId>::const_iterator channelId = wireMap.begin(); channelId != wireMap.end(); channelId++) {
239 vector<DTWireId>::const_iterator
next = channelId;
241 for (vector<DTWireId>::const_iterator ite =
next; ite != wireMap.end(); ite++) {
242 if ((*channelId).layerId() == (*ite).layerId()) {
243 hitsOnSameLayer =
true;
247 if (hitsOnSameLayer) {
249 cout <<
"[DTEfficiencyTask] This RecHit has 2 hits on the same layer, skipping!" << endl;
254 LocalVector phiDirectionInChamber = (*phiSeg).localDirection();
255 if (rPhi && fabs(phiDirectionInChamber.
x() / phiDirectionInChamber.
z()) > 1) {
257 cout <<
" RPhi segment has angle > 45 deg, skipping! " << endl;
258 cout <<
" Theta = " << phiDirectionInChamber.
theta() << endl;
263 LocalVector zDirectionInChamber = (*zSeg).localDirection();
264 if (fabs(zDirectionInChamber.
y() / zDirectionInChamber.
z()) > 1) {
266 cout <<
" RZ segment has angle > 45 deg, skipping! " << endl;
267 cout <<
" Theta = " << zDirectionInChamber.
theta() << endl;
274 if (recHits1D.size() == 10) {
276 cout <<
"[DTEfficiencyTask] 4D Segment with only 10 hits, skipping!" << endl;
281 if ((rPhi && recHits1D.size() == 7) || (rZ && recHits1D.size() == 11)) {
283 if (rPhi && recHits1D.size() == 7)
284 cout <<
"[DTEfficiencyTask] MB4 Segment with only 7 hits!" << endl;
285 if (rZ && recHits1D.size() == 11)
286 cout <<
"[DTEfficiencyTask] 4D Segment with only 11 hits!" << endl;
290 const vector<const DTSuperLayer*>& SupLayers =
chamber->superLayers();
291 map<DTLayerId, bool> layerMap;
292 map<DTWireId, float> wireAndPosInChamberAtLayerZ;
294 for (vector<const DTSuperLayer*>::const_iterator superlayer = SupLayers.begin(); superlayer != SupLayers.end();
296 const vector<const DTLayer*> Layers = (*superlayer)->layers();
297 for (vector<const DTLayer*>::const_iterator
layer = Layers.begin();
layer != Layers.end();
layer++) {
298 layerMap.insert(make_pair((*layer)->id(),
false));
299 const int firstWire = dtGeom->
layer((*layer)->id())->specificTopology().firstChannel();
300 const int lastWire = dtGeom->
layer((*layer)->id())->specificTopology().lastChannel();
301 for (
int i = firstWire;
i - lastWire <= 0;
i++) {
303 float wireX = (*layer)->specificTopology().wirePosition(wireId.wire());
305 GlobalPoint wirePosGlob = (*layer)->toGlobal(wirePosInLay);
307 if ((*superlayer)->id().superlayer() == 1 || (*superlayer)->id().superlayer() == 3) {
308 wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.
x()));
310 wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.
y()));
317 map<DTLayerId, int> NumWireMap;
318 for (vector<DTRecHit1D>::const_iterator
recHit = recHits1D.begin();
recHit != recHits1D.end();
recHit++) {
319 layerMap[(*recHit).wireId().layerId()] =
true;
320 NumWireMap[(*recHit).wireId().layerId()] = (*recHit).wireId().wire();
325 for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
327 missLayerId = (*iter).first;
330 cout <<
"[DTEfficiencyTask] Layer without recHits is: " << missLayerId << endl;
338 LocalPoint segPosAtZLayer = (*segment4D).localPosition() + (*segment4D).localDirection() *
339 missLayerPosInChamber.
z() /
340 cos((*segment4D).localDirection().theta());
345 for (map<DTWireId, float>::const_iterator wireAndPos = wireAndPosInChamberAtLayerZ.begin();
346 wireAndPos != wireAndPosInChamberAtLayerZ.end();
348 DTWireId wireId = (*wireAndPos).first;
349 if (wireId.
layerId() == missLayerId) {
351 if (fabs(segPosAtZLayer.
x() - (*wireAndPos).second) < 2.1)
354 if (fabs(segPosAtZLayer.
y() - (*wireAndPos).second) < 2.1)
360 cout <<
"[DTEfficiencyTask] Cell without hit is: " << missWireId << endl;
363 bool foundUnAssRechit =
false;
366 if (wireAnd1DRecHits.find(missWireId) != wireAnd1DRecHits.end()) {
368 cout <<
"[DTEfficiencyTask] Unassociated Hit found!" << endl;
369 foundUnAssRechit =
true;
372 for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
377 NumWireMap[(*iter).first]);
388 if ((rPhi && recHits1D.size() == 8) || (rZ && recHits1D.size() == 12)) {
389 map<DTLayerId, int> NumWireMap;
391 for (vector<DTRecHit1D>::const_iterator
recHit = recHits1D.begin();
recHit != recHits1D.end();
recHit++) {
392 LayerID = (*recHit).wireId().layerId();
393 NumWireMap[LayerID] = (*recHit).wireId().wire();
395 for (map<DTLayerId, int>::const_iterator iter = NumWireMap.begin(); iter != NumWireMap.end(); iter++) {
399 NumWireMap[(*iter).first]);