Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:30 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Marina Giunta
0005  */
0006 
0007 #include "CalibMuon/DTCalibration/interface/DTTMax.h"
0008 #include "Geometry/DTGeometry/interface/DTLayer.h"
0009 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0010 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0011 // Declare histograms for debugging purposes
0012 //#include "CalibMuon/DTCalibration/interface/Histogram.h"
0013 
0014 #include <map>
0015 #include <iostream>
0016 
0017 using namespace std;
0018 using namespace dttmaxenums;
0019 using namespace DTEnums;
0020 
0021 DTTMax::InfoLayer::InfoLayer(
0022     const DTRecHit1D& rh_, const DTSuperLayer& isl, GlobalVector dir, GlobalPoint pos, const DTTTrigBaseSync& sync)
0023     : rh(rh_), idWire(rh.wireId()), lr(rh.lrSide()) {
0024   const DTLayer* layer = isl.layer(idWire.layerId());
0025   LocalPoint wirePosInLayer(layer->specificTopology().wirePosition(idWire.wire()), 0, 0);
0026   LocalPoint wirePosInSL = isl.toLocal(layer->toGlobal(wirePosInLayer));
0027   wireX = wirePosInSL.x();
0028 
0029   //-- Correction for signal propagation along wire, t0 subtraction,
0030   LocalVector segDir = layer->toLocal(dir);
0031   LocalPoint segPos = layer->toLocal(pos);
0032   LocalPoint segPosAtLayer = segPos + segDir * (-segPos.z()) / cos(segDir.theta());
0033   LocalPoint hitPos(rh.localPosition().x(), segPosAtLayer.y(), 0.);
0034   time = rh.digiTime() - sync.offset(layer, idWire, layer->toGlobal(hitPos));
0035 
0036   if (time < 0. || time > 415.) {
0037     // FIXME introduce time window to reject "out-of-time" digis
0038     cout << " *** WARNING time = " << time << endl;
0039   }
0040 }
0041 
0042 DTTMax::DTTMax(const vector<DTRecHit1D>& hits,
0043                const DTSuperLayer& isl,
0044                GlobalVector dir,
0045                GlobalPoint pos,
0046                const DTTTrigBaseSync& sync,
0047                dtcalibration::Histograms& hist)
0048     : theInfoLayers(4, (InfoLayer*)nullptr),  //FIXME
0049       theTMaxes(4, (TMax*)nullptr) {
0050   // debug parameter for verbose output
0051   debug = "true";
0052 
0053   // Collect all information using InfoLayer
0054   for (vector<DTRecHit1D>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0055     //     cout << "Hit Pos " << (*hit).localPosition() << endl;
0056 
0057     InfoLayer* layInfo = new InfoLayer((*hit), isl, dir, pos, sync);
0058     int ilay = layInfo->idWire.layer();
0059     if (getInfoLayer(ilay) == nullptr) {
0060       getInfoLayer(ilay) = layInfo;
0061     } else {
0062       // FIXME: in case there is > 1 hit/layer, the first is taken and the others are IGNORED.
0063       delete layInfo;
0064     }
0065   }
0066 
0067   // Get the segment direction
0068   theSegDir = ((isl.toLocal(dir).x() < 0) ? L : R);
0069 
0070   int layersIn = 0;
0071   int nGoodHits = 0;
0072   for (vector<InfoLayer*>::const_iterator ilay = theInfoLayers.begin(); ilay != theInfoLayers.end(); ++ilay) {
0073     if ((*ilay) == nullptr) {
0074       theSegType += "X";
0075       continue;
0076     }
0077     DTEnums::DTCellSide lOrR = (*ilay)->lr;
0078     if (lOrR == Left)
0079       theSegType += "L";
0080     else if (lOrR == Right)
0081       theSegType += "R";
0082     else
0083       theSegType += "X";
0084 
0085     // layersIn : 6 =  layers 1,2,3
0086     //            7 =         1,2,4
0087     //            8 =         1,3,4
0088     //            9 =         2,3,4
0089     //            10=         1,2,3,4
0090     layersIn += (*ilay)->idWire.layer();
0091     nGoodHits++;
0092   }
0093 
0094   if (nGoodHits >= 3 && (theSegType != "RRRR" && theSegType != "LLLL")) {
0095     float t1 = 0.;
0096     float t2 = 0.;
0097     float t3 = 0.;
0098     float t4 = 0.;
0099     float x1 = 0.;
0100     float x2 = 0.;
0101     float x3 = 0.;
0102     float x4 = 0.;
0103 
0104     if (layersIn <= 8 || layersIn == 10) {
0105       t1 = getInfoLayer(1)->time;
0106       x1 = getInfoLayer(1)->wireX;
0107     }
0108     if (layersIn <= 7 || layersIn >= 9) {
0109       t2 = getInfoLayer(2)->time;
0110       x2 = getInfoLayer(2)->wireX;
0111     }
0112     if (layersIn == 6 || layersIn >= 8) {
0113       t3 = getInfoLayer(3)->time;
0114       x3 = getInfoLayer(3)->wireX;
0115     }
0116     if (layersIn >= 7) {
0117       t4 = getInfoLayer(4)->time;
0118       x4 = getInfoLayer(4)->wireX;
0119     }
0120 
0121     float t = 0.;
0122     TMaxCells cGroup = notInit;
0123     string type;
0124     SigmaFactor sigma = noR;  // Return the factor relating the width of the TMax distribution and the cell resolution
0125     float halfCell = 2.1;     // 2.1 is the half cell length in cm
0126     float delta = 0.5;        // (diff. wire pos.) < delta, halfCell+delta, .....
0127     unsigned t0Factor = 99;   // "quantity" of Delta(t0) included in the tmax formula
0128 
0129     //Debug
0130     if (debug) {
0131       cout << "seg. type: " << theSegType << " and dir: " << theSegDir << endl;
0132       cout << "t1, t2, t3, t4: " << t1 << " " << t2 << " " << t3 << " " << t4 << endl;
0133       cout << "x1, x2, x3, x4: " << x1 << " " << x2 << " " << x3 << " " << x4 << endl;
0134     }
0135 
0136     //different t0 hists (if you have at least one hit within a certain distance from the wire)
0137     unsigned hSubGroup = 99;  //
0138     if (t1 == 0. || t2 == 0. || t3 == 0. || t4 == 0.)
0139       hSubGroup = 0;  //if only 3 hits
0140     else if (t1 <= 5. || t2 <= 5. || t3 <= 5. || t4 <= 5.)
0141       hSubGroup = 1;  //if distance of one hit from wire < 275um (v_drift=55um/ns)
0142     else if (t1 <= 10. || t2 <= 10. || t3 <= 10. || t4 <= 10.)
0143       hSubGroup = 2;
0144     else if (t1 <= 20. || t2 <= 20. || t3 <= 20. || t4 <= 20.)
0145       hSubGroup = 3;
0146     else if (t1 <= 50. || t2 <= 50. || t3 <= 50. || t4 <= 50.)
0147       hSubGroup = 4;
0148 
0149     if ((layersIn == 6 || layersIn == 10) && (fabs(x1 - x3) < delta)) {
0150       cGroup = c123;
0151       ((type += theSegType[0]) += theSegType[1]) += theSegType[2];
0152       sigma = r32;
0153       if (type == "LRL" || type == "RLR") {
0154         t0Factor = 2;
0155         t = (t1 + t3) / 2. + t2;
0156         hist.hT123LRL->Fill(t);
0157       } else if ((type == "LLR" && theSegDir == R) || (type == "RRL" && theSegDir == L)) {
0158         t0Factor = 1;
0159         t = (t3 - t1) / 2. + t2;
0160         hist.hT123LLR->Fill(t);
0161       } else if ((type == "LRR" && theSegDir == R) || (type == "RLL" && theSegDir == L)) {
0162         t0Factor = 1;
0163         t = (t1 - t3) / 2. + t2;
0164         hist.hT123LRR->Fill(t);
0165       } else {
0166         t = -1.;
0167         sigma = noR;
0168         hist.hT123Bad->Fill(t);
0169       }
0170       theTMaxes[cGroup] = new TMax(t, cGroup, type, sigma, t0Factor, hSubGroup);
0171       if (debug)
0172         cout << "tmax123 " << t << " " << type << endl;
0173     }
0174     if (layersIn == 7 || layersIn == 10) {
0175       cGroup = c124;
0176       type.clear();
0177       sigma = r72;
0178       ((type += theSegType[0]) += theSegType[1]) += theSegType[3];
0179       if ((theSegType == "LRLR" && type == "LRR" && x1 > x4) || (theSegType == "RLRL" && type == "RLL" && x1 < x4)) {
0180         t0Factor = 2;
0181         t = 1.5 * t2 + t1 - t4 / 2.;
0182         hist.hT124LRR1gt4->Fill(t);
0183       } else if ((type == "LLR" && theSegDir == R && (fabs(x2 - x4) < delta) && x1 < x2) ||
0184                  (type == "RRL" && theSegDir == L && (fabs(x2 - x4) < delta) && x1 > x2)) {
0185         t0Factor = 1;
0186         t = 1.5 * t2 - t1 + t4 / 2.;
0187         hist.hT124LLR->Fill(t);
0188       } else if ((type == "LLL" && theSegDir == R && (fabs(x2 - x4) < delta) && x1 < x2) ||
0189                  (type == "RRR" && theSegDir == L && (fabs(x2 - x4) < delta) && x1 > x2)) {
0190         t0Factor = 0;
0191         t = 1.5 * t2 - t1 - t4 / 2.;
0192         hist.hT124LLLR->Fill(t);
0193       } else if ((type == "LLL" && theSegDir == L && (fabs(x2 - x4) < delta)) ||
0194                  (type == "RRR" && theSegDir == R && (fabs(x2 - x4) < delta))) {
0195         t0Factor = 0;
0196         t = -1.5 * t2 + t1 + t4 / 2.;
0197         hist.hT124LLLL->Fill(t);
0198       } else if ((type == "LRL" && theSegDir == L && (fabs(x2 - x4) < delta)) ||
0199                  (type == "RLR" && theSegDir == R && (fabs(x2 - x4) < delta))) {
0200         t0Factor = 3;
0201         t = 1.5 * t2 + t1 + t4 / 2.;
0202         hist.hT124LRLL->Fill(t);
0203       } else if ((type == "LRL" && theSegDir == R && (fabs(x1 - x4) < (halfCell + delta))) ||
0204                  (type == "RLR" && theSegDir == L && (fabs(x1 - x4) < (halfCell + delta)))) {
0205         t0Factor = 99;  // it's actually 1.5, but this value it's not used
0206         t = 3. / 4. * t2 + t1 / 2. + t4 / 4.;
0207         sigma = r78;
0208         hist.hT124LRLR->Fill(t);
0209       } else if ((type == "LRR" && theSegDir == R && x1 < x4 && (fabs(x1 - x4) < (halfCell + delta))) ||
0210                  (type == "RLL" && theSegDir == L && x1 > x4 && (fabs(x1 - x4) < (halfCell + delta)))) {
0211         t0Factor = 1;
0212         t = 3. / 4. * t2 + t1 / 2. - t4 / 4.;
0213         sigma = r78;
0214         hist.hT124LRR1lt4->Fill(t);
0215       } else {
0216         t = -1.;
0217         sigma = noR;
0218         hist.hT124Bad->Fill(t);
0219       }
0220       theTMaxes[cGroup] = new TMax(t, cGroup, type, sigma, t0Factor, hSubGroup);
0221       if (debug)
0222         cout << "tmax124 " << t << " " << t0Factor << " " << type << endl;
0223     }
0224     if (layersIn == 8 || layersIn == 10) {
0225       cGroup = c134;
0226       type.clear();
0227       ((type += theSegType[0]) += theSegType[2]) += theSegType[3];
0228       sigma = r72;
0229       if ((type == "LLR" && x1 > x4 && theSegType == "LRLR") || (type == "RRL" && x1 < x4 && theSegType == "RLRL")) {
0230         t0Factor = 2;
0231         t = 1.5 * t3 + t4 - t1 / 2.;
0232         hist.hT134LLR1gt4->Fill(t);
0233       } else if ((type == "LLR" && x1 < x4 && (fabs(x1 - x4) < (halfCell + delta))) ||
0234                  (type == "RRL" && x1 > x4 && (fabs(x1 - x4) < (halfCell + delta)))) {
0235         t0Factor = 1;
0236         t = 3. / 4. * t3 + t4 / 2. - t1 / 4.;
0237         sigma = r78;
0238         hist.hT134LLR1lt4->Fill(t);
0239       } else if ((type == "LRR" && theSegDir == R && x1 < x4 && (fabs(x1 - x3) < delta)) ||
0240                  (type == "RLL" && theSegDir == L && x1 > x4 && (fabs(x1 - x3) < delta))) {
0241         t0Factor = 1;
0242         t = 1.5 * t3 - t4 + t1 / 2.;
0243         hist.hT134LRR->Fill(t);
0244       } else if ((type == "LRL" && theSegDir == R && (fabs(x1 - x3) < delta)) ||
0245                  (type == "RLR" && theSegDir == L && (fabs(x1 - x3) < delta))) {
0246         t0Factor = 3;
0247         t = 1.5 * t3 + t4 + t1 / 2.;
0248         hist.hT134LRLR->Fill(t);
0249       } else if ((type == "LRL" && theSegDir == L && (fabs(x1 - x3) < (2. * halfCell + delta))) ||
0250                  (type == "RLR" && theSegDir == R && (fabs(x1 - x3) < (2. * halfCell + delta)))) {
0251         t0Factor = 99;  // it's actually 1.5, but this value it's not used
0252         t = 3. / 4. * t3 + t4 / 2. + t1 / 4.;
0253         sigma = r78;
0254         hist.hT134LRLL->Fill(t);
0255       } else if ((type == "LLL" && theSegDir == L && x1 > x4 && (fabs(x1 - x3) < delta)) ||
0256                  (type == "RRR" && theSegDir == R && x1 < x4 && (fabs(x1 - x3) < delta))) {
0257         t0Factor = 0;
0258         t = 1.5 * t3 - t4 - t1 / 2.;
0259         hist.hT134LLLL->Fill(t);
0260       } else if ((type == "LLL" && theSegDir == R && (fabs(x1 - x3) < delta)) ||
0261                  (type == "RRR" && theSegDir == L && (fabs(x1 - x3) < delta))) {
0262         t0Factor = 0;
0263         t = -1.5 * t3 + t4 + t1 / 2.;
0264         hist.hT134LLLR->Fill(t);
0265       } else {
0266         t = -1;
0267         sigma = noR;
0268         hist.hT134Bad->Fill(t);
0269       }
0270       theTMaxes[cGroup] = new TMax(t, cGroup, type, sigma, t0Factor, hSubGroup);
0271       if (debug)
0272         cout << "tmax134 " << t << " " << t0Factor << " " << type << endl;
0273     }
0274     if ((layersIn == 9 || layersIn == 10) && (fabs(x2 - x4) < delta)) {
0275       cGroup = c234;
0276       type.clear();
0277       ((type += theSegType[1]) += theSegType[2]) += theSegType[3];
0278       sigma = r32;
0279       if ((type == "LRL") || (type == "RLR")) {
0280         t0Factor = 2;
0281         t = (t2 + t4) / 2. + t3;
0282         hist.hT234LRL->Fill(t);
0283       } else if ((type == "LRR" && theSegDir == R) || (type == "RLL" && theSegDir == L)) {
0284         t0Factor = 1;
0285         t = (t2 - t4) / 2. + t3;
0286         hist.hT234LRR->Fill(t);
0287       } else if ((type == "LLR" && theSegDir == R) || (type == "RRL" && theSegDir == L)) {
0288         t0Factor = 1;
0289         t = (t4 - t2) / 2. + t3;
0290         hist.hT234LLR->Fill(t);
0291       } else {
0292         t = -1;
0293         sigma = noR;
0294         hist.hT234Bad->Fill(t);
0295       }
0296       theTMaxes[cGroup] = new TMax(t, cGroup, type, sigma, t0Factor, hSubGroup);
0297       if (debug)
0298         cout << "tmax234 " << t << " " << type << endl;
0299     }
0300   }
0301 }
0302 
0303 vector<const DTTMax::TMax*> DTTMax::getTMax(const DTWireId& idWire) {
0304   vector<const TMax*> v;
0305   if (idWire.layer() == 1) {
0306     v.push_back(getTMax(c123));  //FIXME: needs pointer
0307     v.push_back(getTMax(c124));
0308     v.push_back(getTMax(c134));
0309   } else if (idWire.layer() == 2) {
0310     v.push_back(getTMax(c123));
0311     v.push_back(getTMax(c124));
0312     v.push_back(getTMax(c234));
0313   } else if (idWire.layer() == 3) {
0314     v.push_back(getTMax(c123));
0315     v.push_back(getTMax(c134));
0316     v.push_back(getTMax(c234));
0317   } else {
0318     v.push_back(getTMax(c124));
0319     v.push_back(getTMax(c134));
0320     v.push_back(getTMax(c234));
0321   }
0322   return v;
0323 }
0324 
0325 vector<const DTTMax::TMax*> DTTMax::getTMax(const DTSuperLayerId& isl) {
0326   vector<const TMax*> v;
0327   // add TMax* to the vector only if it really exists
0328   if (getTMax(c123))
0329     v.push_back(getTMax(c123));
0330   if (getTMax(c124))
0331     v.push_back(getTMax(c124));
0332   if (getTMax(c134))
0333     v.push_back(getTMax(c134));
0334   if (getTMax(c234))
0335     v.push_back(getTMax(c234));
0336   return v;
0337 }
0338 
0339 const DTTMax::TMax* DTTMax::getTMax(TMaxCells cCase) { return theTMaxes[cCase]; }
0340 
0341 /* Destructor */
0342 DTTMax::~DTTMax() {
0343   for (vector<InfoLayer*>::const_iterator ilay = theInfoLayers.begin(); ilay != theInfoLayers.end(); ++ilay) {
0344     delete (*ilay);
0345   }
0346 
0347   for (vector<TMax*>::const_iterator iTmax = theTMaxes.begin(); iTmax != theTMaxes.end(); ++iTmax) {
0348     delete (*iTmax);
0349   }
0350 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!