File indexing completed on 2024-09-07 04:36:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h"
0014 #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h"
0015 #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h"
0016 #include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h"
0017
0018
0019 #include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h"
0020
0021
0022
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0025 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0026 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0027
0028
0029 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0030 #include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h"
0031
0032
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h"
0036
0037
0038 #include <TFile.h>
0039 #include <TH1F.h>
0040 #include <TH2F.h>
0041
0042 using namespace std;
0043
0044 const double microntocm = 0.0001;
0045
0046 PixelTemplateSmearerBase::PixelTemplateSmearerBase(const std::string& name,
0047 const edm::ParameterSet& config,
0048 edm::ConsumesCollector& consumesCollector)
0049 : TrackingRecHitAlgorithm(name, config, consumesCollector) {
0050
0051 mergeHitsOn = config.getParameter<bool>("MergeHitsOn");
0052 isBarrel = config.getParameter<bool>("isBarrel");
0053 int detType = (isBarrel) ? 1 : 0;
0054
0055
0056 theBigPixelResolutionFileName = config.getParameter<string>("BigPixelResolutionFile");
0057 theBigPixelResolutionFileName = edm::FileInPath(theBigPixelResolutionFileName).fullPath();
0058
0059 theEdgePixelResolutionFileName = config.getParameter<string>("EdgePixelResolutionFile");
0060 theEdgePixelResolutionFileName = edm::FileInPath(theEdgePixelResolutionFileName).fullPath();
0061
0062 theRegularPixelResolutionFileName = config.getParameter<string>("RegularPixelResolutionFile");
0063 theRegularPixelResolutionFileName = edm::FileInPath(theRegularPixelResolutionFileName).fullPath();
0064
0065
0066
0067
0068 int status = 0;
0069 theRegularPixelResolutions = std::make_shared<PixelResolutionHistograms>(theRegularPixelResolutionFileName, "");
0070 if ((status = theRegularPixelResolutions->status()) != 0) {
0071 throw cms::Exception("PixelTemplateSmearerBase:")
0072 << " constructing PixelResolutionHistograms file " << theRegularPixelResolutionFileName
0073 << " failed with status = " << status << std::endl;
0074 }
0075
0076 theBigPixelResolutions = std::make_shared<PixelResolutionHistograms>(
0077 theBigPixelResolutionFileName, "", detType, (!isBarrel), false, true);
0078 if ((status = theBigPixelResolutions->status()) != 0) {
0079 throw cms::Exception("PixelTemplateSmearerBase:")
0080 << " constructing PixelResolutionHistograms file " << theBigPixelResolutionFileName
0081 << " failed with status = " << status << std::endl;
0082 }
0083
0084 theEdgePixelResolutions = std::make_shared<PixelResolutionHistograms>(
0085 theEdgePixelResolutionFileName, "", detType, false, true, true);
0086 if ((status = theEdgePixelResolutions->status()) != 0) {
0087 throw cms::Exception("PixelTemplateSmearerBase:")
0088 << " constructing PixelResolutionHistograms file " << theEdgePixelResolutionFileName
0089 << " failed with status = " << status << std::endl;
0090 }
0091
0092
0093 theMergingProbabilityFileName = config.getParameter<string>("MergingProbabilityFile");
0094 theMergingProbabilityFileName = edm::FileInPath(theMergingProbabilityFileName).fullPath();
0095 theMergingProbabilityFile = std::make_unique<TFile>(theMergingProbabilityFileName.c_str(), "READ");
0096
0097 theMergedPixelResolutionXFileName = config.getParameter<string>("MergedPixelResolutionXFile");
0098 theMergedPixelResolutionXFileName = edm::FileInPath(theMergedPixelResolutionXFileName).fullPath();
0099 theMergedPixelResolutionXFile = std::make_unique<TFile>(theMergedPixelResolutionXFileName.c_str(), "READ");
0100
0101 theMergedPixelResolutionYFileName = config.getParameter<string>("MergedPixelResolutionYFile");
0102 theMergedPixelResolutionYFileName = edm::FileInPath(theMergedPixelResolutionYFileName).fullPath();
0103 theMergedPixelResolutionYFile = std::make_unique<TFile>(theMergedPixelResolutionYFileName.c_str(), "READ");
0104
0105
0106
0107
0108
0109 if (config.exists("templateId")) {
0110
0111 templateId = config.getParameter<int>("templateId");
0112 if (templateId > 0) {
0113 if (!SiPixelTemplate::pushfile(templateId, thePixelTemp_)) {
0114 throw cms::Exception("PixelTemplateSmearerBase:")
0115 << "SiPixel Template " << templateId << " Not Loaded Correctly!" << std::endl;
0116 }
0117 }
0118 }
0119
0120
0121
0122
0123
0124
0125
0126 }
0127
0128 PixelTemplateSmearerBase::~PixelTemplateSmearerBase() {}
0129
0130
0131
0132
0133
0134
0135
0136 void PixelTemplateSmearerBase::beginRun(edm::Run const& run,
0137 const edm::EventSetup& eventSetup,
0138 const SiPixelTemplateDBObject* pixelTemplateDBObjectPtr,
0139 const std::vector<SiPixelTemplateStore>& tempStoreRef) {
0140
0141
0142
0143
0144
0145 if (templateId == -1) {
0146 thePixelTempRef = &tempStoreRef;
0147 pixelTemplateDBObject_ = pixelTemplateDBObjectPtr;
0148 }
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162 }
0163
0164
0165
0166
0167
0168
0169
0170
0171 TrackingRecHitProductPtr PixelTemplateSmearerBase::process(TrackingRecHitProductPtr product) const {
0172 std::vector<std::pair<unsigned int, const PSimHit*>>& simHitIdPairs = product->getSimHitIdPairs();
0173 std::vector<const PSimHit*> simHits(simHitIdPairs.size());
0174 for (unsigned int ihit = 0; ihit < simHitIdPairs.size(); ++ihit) {
0175 simHits[ihit] = simHitIdPairs[ihit].second;
0176 }
0177
0178 RandomEngineAndDistribution const& randomEngine = getRandomEngine();
0179
0180 const GeomDet* geomDet = getTrackerGeometry().idToDetUnit(product->getDetId());
0181 const PixelGeomDetUnit* pixelGeomDet = dynamic_cast<const PixelGeomDetUnit*>(geomDet);
0182 if (pixelGeomDet == nullptr) {
0183 throw cms::Exception("FastSimulation/TrackingRecHitProducer")
0184 << "The GeomDetUnit is not a PixelGeomDetUnit. This should never happen!";
0185 }
0186 const BoundPlane& theDetPlane = pixelGeomDet->surface();
0187 const Bounds& theBounds = theDetPlane.bounds();
0188 const double boundX = theBounds.width() / 2.;
0189 const double boundY = theBounds.length() / 2.;
0190
0191 std::vector<TrackingRecHitProduct::SimHitIdPair> listOfUnmergedHits;
0192 std::vector<MergeGroup*> listOfMergeGroups;
0193 int nHits = simHits.size();
0194
0195
0196 MergeGroup* mergeGroupByHit[nHits];
0197
0198 if (nHits == 0) {
0199 return product;
0200 } else if (nHits == 1) {
0201 listOfUnmergedHits.push_back(simHitIdPairs[0]);
0202 } else {
0203 if (mergeHitsOn) {
0204 for (int i = 0; i < nHits; ++i) {
0205
0206 mergeGroupByHit[i] = nullptr;
0207 }
0208 for (int i = 0; i < nHits - 1; ++i) {
0209 for (int j = i + 1; j < nHits; ++j) {
0210
0211 bool merged = hitsMerge(*simHitIdPairs[i].second, *simHitIdPairs[j].second);
0212
0213 if (merged) {
0214
0215 if (mergeGroupByHit[j] != nullptr) {
0216 if (mergeGroupByHit[i] == nullptr) {
0217 mergeGroupByHit[i] = mergeGroupByHit[j];
0218 mergeGroupByHit[i]->group.push_back(simHitIdPairs[i]);
0219 mergeGroupByHit[i]->smearIt = true;
0220 } else {
0221 if (mergeGroupByHit[i] != mergeGroupByHit[j]) {
0222 for (auto hit_it = mergeGroupByHit[j]->group.begin(); hit_it != mergeGroupByHit[j]->group.end();
0223 ++hit_it) {
0224 mergeGroupByHit[i]->group.push_back(*hit_it);
0225 mergeGroupByHit[i]->smearIt = true;
0226 }
0227
0228
0229 MergeGroup* mgbhj = mergeGroupByHit[j];
0230 for (int k = 0; k < nHits; ++k) {
0231 if (mgbhj == mergeGroupByHit[k]) {
0232
0233 mergeGroupByHit[k] = mergeGroupByHit[i];
0234 }
0235 }
0236 mgbhj->smearIt = false;
0237 mergeGroupByHit[i]->smearIt = true;
0238
0239
0240
0241
0242
0243 }
0244 }
0245 } else {
0246
0247
0248 if (mergeGroupByHit[i] == nullptr) {
0249
0250
0251 mergeGroupByHit[i] = new MergeGroup();
0252 listOfMergeGroups.push_back(mergeGroupByHit[i]);
0253
0254
0255
0256 mergeGroupByHit[i]->group.push_back(simHitIdPairs[i]);
0257 mergeGroupByHit[i]->smearIt = true;
0258 }
0259
0260 mergeGroupByHit[i]->group.push_back(simHitIdPairs[j]);
0261 mergeGroupByHit[i]->smearIt = true;
0262
0263 mergeGroupByHit[j] = mergeGroupByHit[i];
0264
0265 }
0266
0267 }
0268
0269 }
0270
0271
0272
0273
0274
0275
0276
0277 if (mergeGroupByHit[i] == nullptr) {
0278
0279 listOfUnmergedHits.push_back(simHitIdPairs[i]);
0280 }
0281 }
0282 }
0283 else {
0284
0285
0286 for (int i = 0; i < nHits; ++i) {
0287 listOfUnmergedHits.push_back(simHitIdPairs[i]);
0288 }
0289 }
0290 }
0291
0292
0293
0294
0295 product = processUnmergedHits(listOfUnmergedHits, product, pixelGeomDet, boundX, boundY, &randomEngine);
0296
0297 product = processMergeGroups(listOfMergeGroups, product, pixelGeomDet, boundX, boundY, &randomEngine);
0298
0299
0300
0301
0302
0303
0304 for (auto mg_it = listOfMergeGroups.begin(); mg_it != listOfMergeGroups.end(); ++mg_it) {
0305 delete *mg_it;
0306 }
0307
0308 return product;
0309 }
0310
0311
0312
0313
0314 FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit(const PSimHit& simHit,
0315 const PixelGeomDetUnit* detUnit,
0316 const double boundX,
0317 const double boundY,
0318 RandomEngineAndDistribution const* random) const {
0319
0320
0321
0322 LocalVector localDir = simHit.momentumAtEntry();
0323 float locx = localDir.x();
0324 float locy = localDir.y();
0325 float locz = localDir.z();
0326
0327
0328
0329
0330
0331 float cotalpha = locx / locz;
0332 float cotbeta = locy / locz;
0333
0334
0335 int signOfCotalpha = (cotalpha < 0) ? -1 : 1;
0336 int signOfCotbeta = (cotbeta < 0) ? -1 : 1;
0337
0338
0339 cotalpha *= signOfCotalpha;
0340 cotbeta *= signOfCotbeta;
0341
0342 LogDebug("SmearHit") << "LocalVector=" << locx << "," << locy << "," << locz << " momentum=" << localDir.mag()
0343 << " cotalpha=" << cotalpha << ", cotbeta=" << cotbeta;
0344
0345 const PixelTopology* theSpecificTopology = &(detUnit->specificType().specificTopology());
0346 const RectangularPixelTopology* rectPixelTopology = static_cast<const RectangularPixelTopology*>(theSpecificTopology);
0347
0348 const int nrows = theSpecificTopology->nrows();
0349 const int ncolumns = theSpecificTopology->ncolumns();
0350
0351 const Local3DPoint lp = simHit.localPosition();
0352
0353 const MeasurementPoint mp = rectPixelTopology->measurementPosition(lp);
0354 float mpy = mp.y();
0355 float mpx = mp.x();
0356
0357 float pixelCenterY = 0.5 + (int)mpy;
0358 float pixelCenterX = 0.5 + (int)mpx;
0359
0360 const MeasurementPoint mpCenter(pixelCenterX, pixelCenterY);
0361
0362 const Local3DPoint lpCenter = rectPixelTopology->localPosition(mpCenter);
0363
0364
0365 float xtrk = lp.x() - lpCenter.x();
0366 float ytrk = lp.y() - lpCenter.y();
0367
0368 const float ysize = {0.015}, xsize = {0.01};
0369
0370 float yhit = 20. + 8. * (ytrk / ysize);
0371 float xhit = 20. + 8. * (xtrk / xsize);
0372 int ybin = (int)yhit;
0373 int xbin = (int)xhit;
0374 float yfrac = yhit - (float)ybin;
0375 float xfrac = xhit - (float)xbin;
0376
0377 if (ybin < 0)
0378 ybin = 0;
0379 if (ybin > 39)
0380 ybin = 39;
0381 if (xbin < 0)
0382 xbin = 0;
0383 if (xbin > 39)
0384 xbin = 39;
0385
0386 int ID = templateId;
0387 if (templateId == -1) {
0388
0389
0390 ID = pixelTemplateDBObject_->getTemplateID(detUnit->geographicalId());
0391
0392 }
0393
0394
0395 SiPixelTemplate templ(*thePixelTempRef);
0396
0397
0398 templ.interpolate(ID, cotalpha, cotbeta);
0399
0400
0401
0402 float qbin_frac[4];
0403
0404 float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
0405 bool singlex = false, singley = false;
0406 templ.qbin_dist(ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac);
0407 int nqbin;
0408
0409 double xsizeProbability = random->flatShoot();
0410 double ysizeProbability = random->flatShoot();
0411 bool hitbigx = rectPixelTopology->isItBigPixelInX((int)mpx);
0412 bool hitbigy = rectPixelTopology->isItBigPixelInY((int)mpy);
0413
0414 if (hitbigx)
0415 if (xsizeProbability < nx2_frac)
0416 singlex = true;
0417 else
0418 singlex = false;
0419 else if (xsizeProbability < nx1_frac)
0420 singlex = true;
0421 else
0422 singlex = false;
0423
0424 if (hitbigy)
0425 if (ysizeProbability < ny2_frac)
0426 singley = true;
0427 else
0428 singley = false;
0429 else if (ysizeProbability < ny1_frac)
0430 singley = true;
0431 else
0432 singley = false;
0433
0434
0435 double qbinProbability = random->flatShoot();
0436 for (int i = 0; i < 4; ++i) {
0437 nqbin = i;
0438 if (qbinProbability < qbin_frac[i]) {
0439 break;
0440 }
0441 }
0442
0443
0444
0445 float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}};
0446 templ.ytemp(0, 40, ytempl);
0447 templ.xtemp(0, 40, xtempl);
0448
0449 std::vector<double> ytemp(BYSIZE);
0450 for (int i = 0; i < BYSIZE; ++i) {
0451 ytemp[i] = (1. - yfrac) * ytempl[ybin][i] + yfrac * ytempl[ybin + 1][i];
0452 }
0453
0454 std::vector<double> xtemp(BXSIZE);
0455 for (int i = 0; i < BXSIZE; ++i) {
0456 xtemp[i] = (1. - xfrac) * xtempl[xbin][i] + xfrac * xtempl[xbin + 1][i];
0457 }
0458
0459
0460 const float qThreshold = templ.s50() * 2.0;
0461
0462
0463
0464 int offsetX1 = 0, offsetX2 = 0, offsetY1 = 0, offsetY2 = 0;
0465 int firstY, lastY, firstX, lastX;
0466 for (firstY = 0; firstY < BYSIZE; ++firstY) {
0467 bool yCluster = ytemp[firstY] > qThreshold;
0468 if (yCluster) {
0469 offsetY1 = BHY - firstY;
0470 break;
0471 }
0472 }
0473 for (lastY = firstY; lastY < BYSIZE; ++lastY) {
0474 bool yCluster = ytemp[lastY] > qThreshold;
0475 if (!yCluster) {
0476 lastY = lastY - 1;
0477 offsetY2 = lastY - BHY;
0478 break;
0479 }
0480 }
0481
0482 for (firstX = 0; firstX < BXSIZE; ++firstX) {
0483 bool xCluster = xtemp[firstX] > qThreshold;
0484 if (xCluster) {
0485 offsetX1 = BHX - firstX;
0486 break;
0487 }
0488 }
0489 for (lastX = firstX; lastX < BXSIZE; ++lastX) {
0490 bool xCluster = xtemp[lastX] > qThreshold;
0491 if (!xCluster) {
0492 lastX = lastX - 1;
0493 offsetX2 = lastX - BHX;
0494 break;
0495 }
0496 }
0497
0498
0499 Local3DPoint thePosition;
0500 double theShiftInX;
0501 double theShiftInY;
0502 double theShiftInZ;
0503 LocalError theError;
0504 double theErrorX;
0505 double theErrorY;
0506
0507
0508
0509
0510
0511
0512 bool edge, edgex, edgey;
0513
0514
0515 int firstPixelInX = (int)mpx - offsetX1;
0516 int firstPixelInY = (int)mpy - offsetY1;
0517 int lastPixelInX = (int)mpx + offsetX2;
0518 int lastPixelInY = (int)mpy + offsetY2;
0519 firstPixelInX = (firstPixelInX >= 0) ? firstPixelInX : 0;
0520 firstPixelInY = (firstPixelInY >= 0) ? firstPixelInY : 0;
0521 lastPixelInX = (lastPixelInX < nrows) ? lastPixelInX : nrows - 1;
0522 lastPixelInY = (lastPixelInY < ncolumns) ? lastPixelInY : ncolumns - 1;
0523
0524 edgex = rectPixelTopology->isItEdgePixelInX(firstPixelInX) || rectPixelTopology->isItEdgePixelInX(lastPixelInX);
0525 edgey = rectPixelTopology->isItEdgePixelInY(firstPixelInY) || rectPixelTopology->isItEdgePixelInY(lastPixelInY);
0526 edge = edgex || edgey;
0527
0528
0529
0530 bool hasBigPixelInX = rectPixelTopology->containsBigPixelInX(firstPixelInX, lastPixelInX);
0531 bool hasBigPixelInY = rectPixelTopology->containsBigPixelInY(firstPixelInY, lastPixelInY);
0532
0533
0534 float sigmay, sigmax, sy1, sy2, sx1, sx2;
0535 templ.temperrors(ID,
0536 cotalpha,
0537 cotbeta,
0538 nqbin,
0539 sigmay,
0540 sigmax,
0541 sy1,
0542 sy2,
0543 sx1,
0544 sx2
0545 );
0546
0547 if (edge) {
0548 if (edgex && !edgey) {
0549 theErrorX = 23.0 * microntocm;
0550 theErrorY = 39.0 * microntocm;
0551 } else if (!edgex && edgey) {
0552 theErrorX = 24.0 * microntocm;
0553 theErrorY = 96.0 * microntocm;
0554 } else {
0555 theErrorX = 31.0 * microntocm;
0556 theErrorY = 90.0 * microntocm;
0557 }
0558 } else {
0559 if (singlex) {
0560 if (hitbigx) {
0561 theErrorX = sx2 * microntocm;
0562 } else {
0563 theErrorX = sx1 * microntocm;
0564 }
0565 } else {
0566 theErrorX = sigmax * microntocm;
0567 }
0568 if (singley) {
0569 if (hitbigy) {
0570 theErrorY = sy2 * microntocm;
0571 } else {
0572 theErrorY = sy1 * microntocm;
0573 }
0574 } else {
0575 theErrorY = sigmay * microntocm;
0576 }
0577 }
0578
0579
0580 const TrackerGeomDet* misalignmentDetUnit = getMisalignedGeometry().idToDet(detUnit->geographicalId());
0581 const LocalError& alignmentError = misalignmentDetUnit->localAlignmentError();
0582 if (alignmentError.valid()) {
0583 theError = LocalError(
0584 theErrorX * theErrorX + alignmentError.xx(), alignmentError.xy(), theErrorY * theErrorY + alignmentError.yy());
0585 } else {
0586 theError = LocalError(theErrorX * theErrorX, 0.0, theErrorY * theErrorY);
0587 }
0588
0589
0590
0591
0592
0593
0594
0595
0596 shared_ptr<PixelResolutionHistograms> resHistsX = nullptr;
0597 shared_ptr<PixelResolutionHistograms> resHistsY = nullptr;
0598
0599 if (edge) {
0600 resHistsX = resHistsY = theEdgePixelResolutions;
0601 singlex = singley = false;
0602 } else {
0603
0604 if ((singlex && hitbigx) || (isBarrel && hasBigPixelInX)) {
0605 resHistsX = theBigPixelResolutions;
0606 } else {
0607 resHistsX = theRegularPixelResolutions;
0608 }
0609
0610 if ((singley && hitbigy) || (isBarrel && hasBigPixelInY)) {
0611 resHistsY = theBigPixelResolutions;
0612 } else {
0613 resHistsY = theRegularPixelResolutions;
0614 }
0615 }
0616
0617
0618 const SimpleHistogramGenerator* xgen = resHistsX->getGeneratorX(cotalpha, cotbeta, nqbin, singlex);
0619 const SimpleHistogramGenerator* ygen = resHistsY->getGeneratorY(cotalpha, cotbeta, nqbin, singley);
0620
0621
0622 if (!xgen || !ygen) {
0623 throw cms::Exception("FastSimulation/TrackingRecHitProducer")
0624 << "Histogram (cot\alpha=" << cotalpha << ", cot\beta=" << cotbeta << ", nQbin=" << nqbin
0625 << ") was not found for PixelTemplateSmearer. Check if the smearing resolution histogram exists.";
0626 }
0627
0628
0629
0630
0631 unsigned int retry = 0;
0632
0633 do {
0634
0635 theShiftInX = xgen->generate(random);
0636 theShiftInY = ygen->generate(random);
0637
0638
0639 theShiftInX *= signOfCotalpha;
0640 theShiftInY *= signOfCotbeta;
0641
0642 theShiftInZ = 0.0;
0643
0644 thePosition = Local3DPoint(simHit.localPosition().x() + theShiftInX,
0645 simHit.localPosition().y() + theShiftInY,
0646 simHit.localPosition().z() + theShiftInZ);
0647 retry++;
0648 if (retry > 10) {
0649
0650
0651 thePosition = Local3DPoint(simHit.localPosition().x(), simHit.localPosition().y(), simHit.localPosition().z());
0652 break;
0653 }
0654 } while (fabs(thePosition.x()) > boundX || fabs(thePosition.y()) > boundY);
0655
0656 FastSingleTrackerRecHit recHit(thePosition, theError, *detUnit, fastTrackerRecHitType::siPixel);
0657 return recHit;
0658 }
0659
0660
0661
0662
0663 TrackingRecHitProductPtr PixelTemplateSmearerBase::processUnmergedHits(
0664 std::vector<TrackingRecHitProduct::SimHitIdPair>& unmergedHits,
0665 TrackingRecHitProductPtr product,
0666 const PixelGeomDetUnit* detUnit,
0667 const double boundX,
0668 const double boundY,
0669 RandomEngineAndDistribution const* random) const {
0670 for (auto simHitIdPair : unmergedHits) {
0671 FastSingleTrackerRecHit recHit = smearHit(*simHitIdPair.second, detUnit, boundX, boundY, random);
0672 product->addRecHit(recHit, {simHitIdPair});
0673 }
0674 return product;
0675 }
0676
0677
0678
0679
0680 TrackingRecHitProductPtr PixelTemplateSmearerBase::processMergeGroups(std::vector<MergeGroup*>& mergeGroups,
0681 TrackingRecHitProductPtr product,
0682 const PixelGeomDetUnit* detUnit,
0683 const double boundX,
0684 const double boundY,
0685 RandomEngineAndDistribution const* random) const {
0686 for (auto mg_it = mergeGroups.begin(); mg_it != mergeGroups.end(); ++mg_it) {
0687 if ((*mg_it)->smearIt) {
0688 FastSingleTrackerRecHit recHit = smearMergeGroup(*mg_it, detUnit, boundX, boundY, random);
0689 product->addRecHit(recHit, (*mg_it)->group);
0690 }
0691 }
0692 return product;
0693 }
0694
0695
0696
0697
0698 FastSingleTrackerRecHit PixelTemplateSmearerBase::smearMergeGroup(MergeGroup* mg,
0699 const PixelGeomDetUnit* detUnit,
0700 const double boundX,
0701 const double boundY,
0702 RandomEngineAndDistribution const* random) const {
0703 float loccx = 0;
0704 float loccy = 0;
0705 float loccz = 0;
0706 float nHit = 0;
0707 float locpx = 0;
0708 float locpy = 0;
0709 float locpz = 0;
0710
0711 for (auto hit_it = mg->group.begin(); hit_it != mg->group.end(); ++hit_it) {
0712 const PSimHit simHit = *hit_it->second;
0713
0714 LocalVector localDir = simHit.momentumAtEntry().unit();
0715 loccx += localDir.x();
0716 loccy += localDir.y();
0717 loccz += localDir.z();
0718
0719 const Local3DPoint lpos = simHit.localPosition();
0720 locpx += lpos.x();
0721 locpy += lpos.y();
0722 locpz += lpos.z();
0723
0724 nHit += 1;
0725 }
0726
0727 float locx = loccx / nHit;
0728 float locy = loccy / nHit;
0729 float locz = loccz / nHit;
0730
0731
0732
0733
0734
0735 float cotalpha = locx / locz;
0736 float cotbeta = locy / locz;
0737
0738
0739 int signOfCotalpha = (cotalpha < 0) ? -1 : 1;
0740 int signOfCotbeta = (cotbeta < 0) ? -1 : 1;
0741
0742
0743 cotalpha *= signOfCotalpha;
0744 cotbeta *= signOfCotbeta;
0745
0746 float lpx = locpx / nHit;
0747 float lpy = locpy / nHit;
0748 float lpz = locpz / nHit;
0749
0750
0751 float xtrk = lpx;
0752 float ytrk = lpy;
0753
0754 const float ysize = {0.015}, xsize = {0.01};
0755
0756 float yhit = 20. + 8. * (ytrk / ysize);
0757 float xhit = 20. + 8. * (xtrk / xsize);
0758 int ybin = (int)yhit;
0759 int xbin = (int)xhit;
0760 float yfrac = yhit - (float)ybin;
0761 float xfrac = xhit - (float)xbin;
0762
0763 if (ybin < 0)
0764 ybin = 0;
0765 if (ybin > 39)
0766 ybin = 39;
0767 if (xbin < 0)
0768 xbin = 0;
0769 if (xbin > 39)
0770 xbin = 39;
0771
0772 int ID = templateId;
0773 if (templateId == -1) {
0774
0775
0776 ID = pixelTemplateDBObject_->getTemplateID(detUnit->geographicalId());
0777
0778 }
0779
0780
0781 SiPixelTemplate templ(*thePixelTempRef);
0782
0783
0784 templ.interpolate(ID, cotalpha, cotbeta);
0785
0786
0787
0788 float qbin_frac[4];
0789
0790 float ny1_frac, ny2_frac, nx1_frac, nx2_frac;
0791 bool singlex = false, singley = false;
0792 templ.qbin_dist(ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac);
0793 int nqbin;
0794
0795
0796
0797 bool hitbigx = false;
0798 bool hitbigy = false;
0799
0800
0801
0802 double qbinProbability = random->flatShoot();
0803 for (int i = 0; i < 4; ++i) {
0804 nqbin = i;
0805 if (qbinProbability < qbin_frac[i])
0806 break;
0807 }
0808
0809
0810
0811 float ytempl[41][BYSIZE] = {{0}}, xtempl[41][BXSIZE] = {{0}};
0812 templ.ytemp(0, 40, ytempl);
0813 templ.xtemp(0, 40, xtempl);
0814
0815 std::vector<double> ytemp(BYSIZE);
0816 for (int i = 0; i < BYSIZE; ++i) {
0817 ytemp[i] = (1. - yfrac) * ytempl[ybin][i] + yfrac * ytempl[ybin + 1][i];
0818 }
0819
0820 std::vector<double> xtemp(BXSIZE);
0821 for (int i = 0; i < BXSIZE; ++i) {
0822 xtemp[i] = (1. - xfrac) * xtempl[xbin][i] + xfrac * xtempl[xbin + 1][i];
0823 }
0824
0825
0826 Local3DPoint thePosition;
0827 double theShiftInX;
0828 double theShiftInY;
0829 double theShiftInZ;
0830 LocalError theError;
0831 double theErrorX;
0832 double theErrorY;
0833
0834
0835
0836
0837
0838
0839 bool edge = false;
0840 bool edgex = false;
0841 bool edgey = false;
0842
0843
0844 float sigmay, sigmax, sy1, sy2, sx1, sx2;
0845 templ.temperrors(ID,
0846 cotalpha,
0847 cotbeta,
0848 nqbin,
0849 sigmay,
0850 sigmax,
0851 sy1,
0852 sy2,
0853 sx1,
0854 sx2);
0855
0856
0857 if (edge) {
0858 if (edgex && !edgey) {
0859 theErrorX = 23.0 * microntocm;
0860 theErrorY = 39.0 * microntocm;
0861 } else if (!edgex && edgey) {
0862 theErrorX = 24.0 * microntocm;
0863 theErrorY = 96.0 * microntocm;
0864 } else {
0865 theErrorX = 31.0 * microntocm;
0866 theErrorY = 90.0 * microntocm;
0867 }
0868
0869 } else {
0870 if (singlex) {
0871 if (hitbigx) {
0872 theErrorX = sx2 * microntocm;
0873 } else {
0874 theErrorX = sx1 * microntocm;
0875 }
0876 } else {
0877 theErrorX = sigmax * microntocm;
0878 }
0879
0880 if (singley) {
0881 if (hitbigy) {
0882 theErrorY = sy2 * microntocm;
0883 } else {
0884 theErrorY = sy1 * microntocm;
0885 }
0886 } else {
0887 theErrorY = sigmay * microntocm;
0888 }
0889 }
0890
0891 theError = LocalError(theErrorX * theErrorX, 0., theErrorY * theErrorY);
0892
0893 unsigned int retry = 0;
0894 do {
0895 const SimpleHistogramGenerator* xgen =
0896 new SimpleHistogramGenerator((TH1F*)theMergedPixelResolutionXFile->Get("th1x"));
0897 const SimpleHistogramGenerator* ygen =
0898 new SimpleHistogramGenerator((TH1F*)theMergedPixelResolutionYFile->Get("th1y"));
0899
0900
0901 theShiftInX = xgen->generate(random);
0902 theShiftInY = ygen->generate(random);
0903
0904
0905 theShiftInX *= signOfCotalpha;
0906 theShiftInY *= signOfCotbeta;
0907
0908 theShiftInZ = 0.0;
0909
0910 thePosition = Local3DPoint(lpx + theShiftInX, lpy + theShiftInY, lpz + theShiftInZ);
0911
0912 retry++;
0913 if (retry > 10) {
0914
0915
0916 thePosition = Local3DPoint(lpx, lpy, lpz);
0917 break;
0918 }
0919 } while (fabs(thePosition.x()) > boundX || fabs(thePosition.y()) > boundY);
0920
0921 FastSingleTrackerRecHit recHit(thePosition, theError, *detUnit, fastTrackerRecHitType::siPixel);
0922 return recHit;
0923 }
0924
0925 bool PixelTemplateSmearerBase::hitsMerge(const PSimHit& simHit1, const PSimHit& simHit2) const {
0926 LocalVector localDir = simHit1.momentumAtEntry().unit();
0927 float locy1 = localDir.y();
0928 float locz1 = localDir.z();
0929 float cotbeta = locy1 / locz1;
0930 float loceta = fabs(-log((double)(-cotbeta + sqrt((double)(1. + cotbeta * cotbeta)))));
0931
0932 const Local3DPoint lp1 = simHit1.localPosition();
0933 const Local3DPoint lp2 = simHit2.localPosition();
0934 float lpy1 = lp1.y();
0935 float lpx1 = lp1.x();
0936 float lpy2 = lp2.y();
0937 float lpx2 = lp2.x();
0938 float locdis = 10000. * sqrt(pow(lpx1 - lpx2, 2) + pow(lpy1 - lpy2, 2));
0939 TH2F* probhisto = (TH2F*)theMergingProbabilityFile->Get("h2bc");
0940 float prob =
0941 probhisto->GetBinContent(probhisto->GetXaxis()->FindFixBin(locdis), probhisto->GetYaxis()->FindFixBin(loceta));
0942 return prob > 0;
0943 }