File indexing completed on 2024-04-06 12:33:49
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 #include <memory>
0032
0033
0034 #include "FWCore/Framework/interface/Frameworkfwd.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Utilities/interface/EDGetToken.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "FWCore/ServiceRegistry/interface/Service.h"
0041
0042 #include "TH1F.h"
0043 #include "TH2F.h"
0044 #include "TROOT.h"
0045
0046 #include "TObject.h"
0047 #include "TH1D.h"
0048 #include "TTree.h"
0049 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0050
0051
0052 #include "DataFormats/Common/interface/DetSet.h"
0053 #include "DataFormats/Common/interface/DetSetVector.h"
0054 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0055 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0056 #include "DataFormats/VertexReco/interface/Vertex.h"
0057 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0058 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0059 #include "DataFormats/DetId/interface/DetId.h"
0060 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0061 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0062 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0063 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0064 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0065 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0066 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0067 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0068 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0069 #include "Geometry/Records/interface/TrackerTopologyRcd.h" // PR#7802
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081 class StripClusterMCanalysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0082 public:
0083 explicit StripClusterMCanalysis(const edm::ParameterSet&);
0084 ~StripClusterMCanalysis() override;
0085
0086 private:
0087 void beginJob() override;
0088 void analyze(const edm::Event&, const edm::EventSetup&) override;
0089
0090 typedef std::pair<unsigned int, unsigned int> simHitCollectionID;
0091 typedef std::pair<simHitCollectionID, unsigned int> simhitAddr;
0092 typedef std::map<simHitCollectionID, std::vector<PSimHit> > simhit_collectionMap;
0093 simhit_collectionMap SimHitCollMap_;
0094
0095 void makeMap(const edm::Event& theEvent);
0096
0097
0098
0099 private:
0100 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0101 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0102 edm::ParameterSet conf_;
0103 edm::InputTag beamSpotLabel_;
0104 edm::InputTag primaryVertexLabel_;
0105 edm::InputTag stripClusterSourceLabel_;
0106 edm::InputTag stripSimLinkLabel_;
0107 int printOut_;
0108 edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0109 edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken_;
0110 edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > clusterToken_;
0111 edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink> > stripSimlinkToken_;
0112 std::vector<edm::EDGetTokenT<CrossingFrame<PSimHit> > > cfTokens_;
0113 std::vector<edm::EDGetTokenT<std::vector<PSimHit> > > simHitTokens_;
0114
0115 edm::Handle<edm::DetSetVector<StripDigiSimLink> > stripdigisimlink;
0116
0117 int evCnt_;
0118
0119 TTree* clusTree_;
0120 struct ClusterStruct {
0121 int subDet;
0122 float thickness;
0123 int width;
0124 int NsimHits;
0125 int firstProcess;
0126 int secondProcess;
0127 int thirdProcess;
0128 int fourthProcess;
0129 int firstPID;
0130 int secondPID;
0131 int thirdPID;
0132 int fourthPID;
0133 int Ntp;
0134 float firstTkChg;
0135 float secondTkChg;
0136 float thirdTkChg;
0137 float fourthTkChg;
0138 float charge;
0139 float firstPmag;
0140 float secondPmag;
0141 float thirdPmag;
0142 float fourthPmag;
0143 float firstPathLength;
0144 float secondPathLength;
0145 float thirdPathLength;
0146 float fourthPathLength;
0147 float pathLstraight;
0148 float allHtPathLength;
0149 float Eloss;
0150 int sat;
0151 int tkFlip;
0152 int ovlap;
0153 int layer;
0154 int stereo;
0155 } clusNtp_;
0156
0157 TTree* pvTree_;
0158 struct pvStruct {
0159 int isValid;
0160 int isFake;
0161 int Ntrks;
0162 int nDoF;
0163 float chisq;
0164 float xV;
0165 float yV;
0166 float zV;
0167 float xVsig;
0168 float yVsig;
0169 float zVsig;
0170 } pvNtp_;
0171
0172 TH1F* hNpv;
0173 TH2F* hzV_Iev;
0174 TH2F* hNtrk_zVerrPri;
0175 TH2F* hNtrk_zVerrSec;
0176 TH1F* hZvPri;
0177 TH1F* hZvSec;
0178 };
0179
0180
0181
0182
0183 StripClusterMCanalysis::StripClusterMCanalysis(const edm::ParameterSet& iConfig)
0184 : topoToken_(esConsumes()),
0185 tkGeomToken_(esConsumes()),
0186 conf_(iConfig),
0187 beamSpotLabel_(iConfig.getParameter<edm::InputTag>("beamSpot")),
0188 primaryVertexLabel_(iConfig.getParameter<edm::InputTag>("primaryVertex")),
0189 stripClusterSourceLabel_(iConfig.getParameter<edm::InputTag>("stripClusters")),
0190 stripSimLinkLabel_(iConfig.getParameter<edm::InputTag>("stripSimLinks")),
0191 printOut_(iConfig.getUntrackedParameter<int>("printOut")) {
0192
0193 usesResource(TFileService::kSharedResource);
0194
0195 beamSpotToken_ = consumes<reco::BeamSpot>(beamSpotLabel_);
0196 primaryVertexToken_ = consumes<reco::VertexCollection>(primaryVertexLabel_);
0197 clusterToken_ = consumes<edmNew::DetSetVector<SiStripCluster> >(stripClusterSourceLabel_);
0198 stripSimlinkToken_ = consumes<edm::DetSetVector<StripDigiSimLink> >(stripSimLinkLabel_);
0199
0200 std::vector<std::string> trackerContainers(iConfig.getParameter<std::vector<std::string> >("ROUList"));
0201 cfTokens_.reserve(trackerContainers.size());
0202 simHitTokens_.reserve(trackerContainers.size());
0203 for (auto const& trackerContainer : trackerContainers) {
0204 cfTokens_.push_back(consumes<CrossingFrame<PSimHit> >(edm::InputTag("mix", trackerContainer)));
0205 simHitTokens_.push_back(consumes<std::vector<PSimHit> >(edm::InputTag("g4SimHits", trackerContainer)));
0206 }
0207 }
0208
0209 StripClusterMCanalysis::~StripClusterMCanalysis() = default;
0210
0211
0212
0213
0214
0215
0216 void StripClusterMCanalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0217 using namespace edm;
0218 using namespace std;
0219
0220 if (printOut_ > 0)
0221 std::cout << std::endl;
0222
0223 evCnt_++;
0224
0225
0226 const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0227
0228
0229 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0230 iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0231 const reco::BeamSpot& bs = *recoBeamSpotHandle;
0232 reco::Vertex::Point estPriVtx(bs.position());
0233
0234
0235 edm::Handle<reco::VertexCollection> pixelVertexCollectionHandle;
0236 if (iEvent.getByToken(primaryVertexToken_, pixelVertexCollectionHandle)) {
0237 const reco::VertexCollection pixelVertexColl = *(pixelVertexCollectionHandle.product());
0238 hNpv->Fill(int(pixelVertexColl.size()));
0239 float zv = 0, zverr = 0;
0240 reco::VertexCollection::const_iterator vi = pixelVertexColl.begin();
0241 if (vi != pixelVertexColl.end() && vi->isValid()) {
0242 estPriVtx = vi->position();
0243 zv = vi->z();
0244 zverr = vi->zError();
0245 }
0246 hZvPri->Fill(zv);
0247
0248 int iVtx = 0;
0249 for (reco::VertexCollection::const_iterator vi = pixelVertexColl.begin(); vi != pixelVertexColl.end(); ++vi) {
0250 if (printOut_ > 0)
0251 std::cout << " " << vi->tracksSize() << " " << vi->z() << "+/-" << vi->zError() << std::endl;
0252 pvNtp_.isValid = int(vi->isValid());
0253 pvNtp_.isFake = int(vi->isFake());
0254 pvNtp_.Ntrks = vi->tracksSize();
0255 pvNtp_.nDoF = vi->ndof();
0256 pvNtp_.chisq = vi->chi2();
0257 pvNtp_.xV = vi->x();
0258 pvNtp_.yV = vi->y();
0259 pvNtp_.zV = vi->z();
0260 pvNtp_.xVsig = vi->xError();
0261 pvNtp_.yVsig = vi->yError();
0262 pvNtp_.zVsig = vi->zError();
0263 pvTree_->Fill();
0264 hzV_Iev->Fill(vi->z(), evCnt_);
0265 if (iVtx == 0) {
0266 hNtrk_zVerrPri->Fill(vi->zError(), vi->tracksSize());
0267 } else {
0268 hNtrk_zVerrSec->Fill(vi->zError(), vi->tracksSize());
0269 hZvSec->Fill(vi->z() - zv);
0270 }
0271 ++iVtx;
0272 }
0273 if (printOut_ > 0)
0274 std::cout << " zv = " << zv << " +/- " << zverr << std::endl;
0275 } else {
0276 if (printOut_ > 0)
0277 std::cout << "No vertices found." << std::endl;
0278 }
0279
0280
0281 makeMap(iEvent);
0282
0283
0284
0285
0286 Handle<edmNew::DetSetVector<SiStripCluster> > dsv_SiStripCluster;
0287 iEvent.getByToken(clusterToken_, dsv_SiStripCluster);
0288
0289 iEvent.getByToken(stripSimlinkToken_, stripdigisimlink);
0290 edm::ESHandle<TrackerGeometry> tkgeom = iSetup.getHandle(tkGeomToken_);
0291 if (!tkgeom.isValid()) {
0292 std::cout << "Unable to find TrackerDigiGeometryRecord in event!";
0293 return;
0294 }
0295 const TrackerGeometry& tracker(*tkgeom);
0296
0297
0298 int clusterCount = 0;
0299 int clustersNoDigiSimLink = 0;
0300 for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = dsv_SiStripCluster->begin();
0301 DSViter != dsv_SiStripCluster->end();
0302 DSViter++) {
0303 uint32_t detID = DSViter->id();
0304 DetId theDet(detID);
0305 int subDet = theDet.subdetId();
0306
0307
0308 const StripGeomDetUnit* DetUnit = nullptr;
0309 DetUnit = (const StripGeomDetUnit*)tracker.idToDetUnit(theDet);
0310 float thickness = 0;
0311 if (DetUnit != nullptr)
0312 thickness = DetUnit->surface().bounds().thickness();
0313 int layer = 0;
0314 int stereo = 0;
0315 if (subDet == int(StripSubdetector::TIB)) {
0316 layer = tTopo->tibLayer(detID);
0317 if (tTopo->tibIsStereo(detID))
0318 stereo = 1;
0319 }
0320 if (printOut_ > 0) {
0321 std::cout << tTopo->print(detID);
0322 std::cout << " thickness " << thickness << std::endl;
0323 }
0324
0325 LocalPoint locVtx = DetUnit->toLocal(GlobalPoint(estPriVtx.X(), estPriVtx.Y(), estPriVtx.Z()));
0326 float modPathLength = fabs(thickness * locVtx.mag() / locVtx.z());
0327 if (printOut_ > 0) {
0328 std::cout << " Module at " << DetUnit->position() << std::endl;
0329 std::cout << " PriVtx at " << locVtx;
0330 printf("%s%7.4f%s%4d%s\n", " path ", modPathLength, ", ", DSViter->size(), " clusters");
0331 }
0332
0333 edm::DetSetVector<StripDigiSimLink>::const_iterator isearch = stripdigisimlink->find(detID);
0334
0335
0336 for (edmNew::DetSet<SiStripCluster>::const_iterator ClusIter = DSViter->begin(); ClusIter != DSViter->end();
0337 ClusIter++) {
0338 clusterCount++;
0339
0340
0341 if (isearch == stripdigisimlink->end()) {
0342 clustersNoDigiSimLink++;
0343 if (printOut_ > 0)
0344 std::cout << "No digisimlinks for this module" << std::endl;
0345 continue;
0346 }
0347
0348 const SiStripCluster* clust = ClusIter;
0349 auto const& amp = clust->amplitudes();
0350 int clusiz = amp.size();
0351 int first = clust->firstStrip();
0352 int last = first + clusiz;
0353 float charge = 0;
0354 bool saturates = false;
0355 for (int i = 0; i < clusiz; ++i) {
0356 charge += amp[i];
0357 if (amp[i] == 254 || amp[i] == 255)
0358 saturates = true;
0359 }
0360 if (printOut_ > 0) {
0361 std::cout << "Cluster (width, first, last) = (" << clusiz << ", " << first << ", " << last - 1
0362 << ") charge = " << charge;
0363 if (saturates)
0364 std::cout << " (saturates)";
0365 std::cout << endl;
0366 }
0367 float clusEloss = 0;
0368 int tkFlip = 0;
0369 int ovlap = 0;
0370 std::vector<unsigned int> trackID;
0371 std::vector<simhitAddr> CFaddr;
0372 std::vector<unsigned int> hitProcess;
0373 std::vector<int> hitPID;
0374 std::vector<float> trackCharge;
0375 std::vector<float> hitPmag;
0376 std::vector<float> hitPathLength;
0377
0378 int prevIdx = -1;
0379 edm::DetSet<StripDigiSimLink> link_detset = (*isearch);
0380 for (edm::DetSet<StripDigiSimLink>::const_iterator linkiter = link_detset.data.begin(),
0381 linkEnd = link_detset.data.end();
0382 linkiter != linkEnd;
0383 ++linkiter) {
0384 int theChannel = linkiter->channel();
0385 if (printOut_ > 3)
0386 printf(" %s%4d%s\n", "channel = ", theChannel, " before matching to cluster");
0387 if (theChannel >= first && theChannel < last) {
0388
0389 int stripIdx = theChannel - first;
0390 uint16_t rawAmpl = amp[stripIdx];
0391 if (printOut_ > 1)
0392 printf(" %s%4d%s%8d%s%8d%s%2d%s%8d%s%8.4f%s%5d\n",
0393 "channel = ",
0394 linkiter->channel(),
0395 " TrackID = ",
0396 linkiter->SimTrackId(),
0397 " EventID = ",
0398 linkiter->eventId().rawId(),
0399 " TofBin = ",
0400 linkiter->TofBin(),
0401 " CFPos = ",
0402 linkiter->CFposition(),
0403 " fraction = ",
0404 linkiter->fraction(),
0405 " amp = ",
0406 rawAmpl);
0407
0408
0409
0410
0411 unsigned int thisTrackID = linkiter->SimTrackId();
0412
0413 if (stripIdx == prevIdx)
0414 ovlap = 1;
0415
0416 bool newTrack = true;
0417 if (std::find(trackID.begin(), trackID.end(), thisTrackID) != trackID.end())
0418 newTrack = false;
0419 if (newTrack) {
0420
0421 trackID.push_back(thisTrackID);
0422
0423 trackCharge.push_back(linkiter->fraction() * rawAmpl);
0424 } else {
0425 for (unsigned int i = 0; i < trackID.size(); ++i)
0426 if (trackID[i] == thisTrackID)
0427 trackCharge[i] += linkiter->fraction() * rawAmpl;
0428 }
0429 if (printOut_ > 2) {
0430 std::cout << " Track charge accumulator = ";
0431 for (unsigned int it = 0; it < trackCharge.size(); ++it)
0432 printf("%7.1f ", trackCharge[it]);
0433 std::cout << std::endl;
0434 }
0435
0436
0437
0438
0439 unsigned int currentCFPos = linkiter->CFposition();
0440 unsigned int tofBin = linkiter->TofBin();
0441 simHitCollectionID theSimHitCollID = std::make_pair(theDet.subdetId(), tofBin);
0442 simhitAddr currentAddr = std::make_pair(theSimHitCollID, currentCFPos);
0443 bool newHit = true;
0444 if (std::find(CFaddr.begin(), CFaddr.end(), currentAddr) != CFaddr.end())
0445 newHit = false;
0446 if (newHit) {
0447 simhit_collectionMap::const_iterator it = SimHitCollMap_.find(theSimHitCollID);
0448 if (it != SimHitCollMap_.end()) {
0449 if (currentCFPos < (it->second).size()) {
0450 const PSimHit& theSimHit = (it->second)[currentCFPos];
0451 CFaddr.push_back(currentAddr);
0452 hitProcess.push_back(theSimHit.processType());
0453 hitPID.push_back(theSimHit.particleType());
0454 hitPmag.push_back(theSimHit.pabs());
0455 Local3DPoint entry = theSimHit.entryPoint();
0456 Local3DPoint exit = theSimHit.exitPoint();
0457 Local3DVector segment = exit - entry;
0458 hitPathLength.push_back(segment.mag());
0459 clusEloss += theSimHit.energyLoss();
0460 if (printOut_ > 1 && theSimHit.pabs() >= 1.0)
0461 printf(" %s%3d%s%3d%s%5d%s%7.2f%s%10.6f%s%8.4f%s%8.4f\n",
0462 "SimHit ",
0463 int(CFaddr.size()),
0464 ", process = ",
0465 theSimHit.processType(),
0466 ", PID = ",
0467 theSimHit.particleType(),
0468 ", p = ",
0469 theSimHit.pabs(),
0470 ", Eloss = ",
0471 theSimHit.energyLoss(),
0472 ", segment = ",
0473 segment.mag(),
0474 ", str segment = ",
0475 modPathLength);
0476 } else {
0477 std::cout << "currentCFPos " << currentCFPos << " is out of range for " << (it->second).size()
0478 << std::endl;
0479 }
0480 }
0481 }
0482 prevIdx = stripIdx;
0483 }
0484 }
0485 if (prevIdx == -1)
0486 continue;
0487
0488 if (ovlap && trackID[1] < trackID[0])
0489 tkFlip = 1;
0490
0491
0492
0493
0494 clusNtp_.subDet = subDet;
0495 clusNtp_.thickness = thickness;
0496 clusNtp_.width = amp.size();
0497 clusNtp_.NsimHits = CFaddr.size();
0498 clusNtp_.firstProcess = !hitProcess.empty() ? hitProcess[0] : -1;
0499 clusNtp_.secondProcess = hitProcess.size() > 1 ? hitProcess[1] : -1;
0500 clusNtp_.thirdProcess = hitProcess.size() > 2 ? hitProcess[2] : -1;
0501 clusNtp_.fourthProcess = hitProcess.size() > 3 ? hitProcess[3] : -1;
0502 clusNtp_.firstPID = !hitPID.empty() ? hitPID[0] : 0;
0503 clusNtp_.secondPID = hitPID.size() > 1 ? hitPID[1] : 0;
0504 clusNtp_.thirdPID = hitPID.size() > 2 ? hitPID[2] : 0;
0505 clusNtp_.fourthPID = hitPID.size() > 3 ? hitPID[3] : 0;
0506 clusNtp_.firstPmag = !hitPmag.empty() ? hitPmag[0] : 0;
0507 clusNtp_.secondPmag = hitPmag.size() > 1 ? hitPmag[1] : 0;
0508 clusNtp_.thirdPmag = hitPmag.size() > 2 ? hitPmag[2] : 0;
0509 clusNtp_.fourthPmag = hitPmag.size() > 3 ? hitPmag[3] : 0;
0510 clusNtp_.firstPathLength = !hitPathLength.empty() ? hitPathLength[0] : 0;
0511 clusNtp_.secondPathLength = hitPathLength.size() > 1 ? hitPathLength[1] : 0;
0512 clusNtp_.thirdPathLength = hitPathLength.size() > 2 ? hitPathLength[2] : 0;
0513 clusNtp_.fourthPathLength = hitPathLength.size() > 3 ? hitPathLength[3] : 0;
0514 clusNtp_.pathLstraight = modPathLength;
0515 float allHtPathLength = 0;
0516 for (unsigned int ih = 0; ih < hitPathLength.size(); ++ih)
0517 allHtPathLength += hitPathLength[ih];
0518 clusNtp_.allHtPathLength = allHtPathLength;
0519 clusNtp_.Ntp = trackID.size();
0520 if (printOut_ > 0 && trackCharge.size() == 2)
0521 cout << " charge 1st, 2nd, dE/dx 1st, 2nd, asymmetry = " << trackCharge[0] << " " << trackCharge[1] << " "
0522 << 3.36e-4 * trackCharge[0] / modPathLength << " " << 3.36e-4 * trackCharge[1] / modPathLength << " "
0523 << (trackCharge[0] - trackCharge[1]) / (trackCharge[0] + trackCharge[1]) << " " << tkFlip << endl;
0524 clusNtp_.firstTkChg = !trackCharge.empty() ? trackCharge[0] : 0;
0525 clusNtp_.secondTkChg = trackCharge.size() > 1 ? trackCharge[1] : 0;
0526 clusNtp_.thirdTkChg = trackCharge.size() > 2 ? trackCharge[2] : 0;
0527 clusNtp_.fourthTkChg = trackCharge.size() > 3 ? trackCharge[3] : 0;
0528 clusNtp_.charge = charge;
0529 clusNtp_.Eloss = clusEloss;
0530 clusNtp_.sat = saturates ? 1 : 0;
0531 clusNtp_.tkFlip = tkFlip;
0532 clusNtp_.ovlap = ovlap;
0533 clusNtp_.layer = layer;
0534 clusNtp_.stereo = stereo;
0535 clusTree_->Fill();
0536 }
0537 }
0538 if (printOut_ > 0)
0539 std::cout << clusterCount << " total clusters; " << clustersNoDigiSimLink << " without digiSimLinks" << std::endl;
0540 }
0541
0542
0543 void StripClusterMCanalysis::beginJob() {
0544 int bufsize = 64000;
0545 edm::Service<TFileService> fs;
0546
0547 clusTree_ = fs->make<TTree>("ClusterNtuple", "Cluster analyzer ntuple");
0548 clusTree_->Branch(
0549 "cluster",
0550 &clusNtp_,
0551 "subDet/I:thickness/F:width/"
0552 "I:NsimHits:firstProcess:secondProcess:thirdProcess:fourthProcess:firstPID:secondPID:thirdPID:fourthPID:Ntp:"
0553 "firstTkChg/F:secondTkChg:thirdTkChg/"
0554 "F:fourthTkChg:charge:firstPmag:secondPmag:thirdPmag:fourthPmag:firstPathLength:secondPathLength:thirdPathLength:"
0555 "fourthPathLength:pathLstraight:allHtPathLength:Eloss:sat/I:tkFlip:ovlap:layer:stereo",
0556 bufsize);
0557
0558 pvTree_ = fs->make<TTree>("pVertexNtuple", "Pixel vertex ntuple");
0559 pvTree_->Branch("pVertex", &pvNtp_, "isValid/I:isFake:Ntrks:nDoF:chisq/F:xV:yV:zV:xVsig:yVsig:zVsig", bufsize);
0560
0561
0562 hNpv = fs->make<TH1F>("hNpv", "No. of pixel vertices", 40, 0, 40);
0563 hzV_Iev = fs->make<TH2F>("hzV_Iev", "Zvertex vs event index", 20, -10, 10, 100, 0, 100);
0564 hNtrk_zVerrPri = fs->make<TH2F>("hzVerr_NtrkPri", "Ntracks vs Zvertex error, primary", 50, 0, 0.025, 30, 0, 150);
0565 hNtrk_zVerrSec = fs->make<TH2F>("hzVerr_NtrkSec", "Ntracks vs Zvertex error, secondary", 50, 0, 0.025, 30, 0, 150);
0566 hZvPri = fs->make<TH1F>("hZvPri", "Zvertex, primary", 48, -15, 15);
0567 hZvSec = fs->make<TH1F>("hZvSec", "Zvertex, secondary", 48, -15, 15);
0568 evCnt_ = 0;
0569 }
0570
0571 void StripClusterMCanalysis::makeMap(const edm::Event& theEvent) {
0572
0573
0574
0575
0576 SimHitCollMap_.clear();
0577 for (auto const& cfToken : cfTokens_) {
0578 edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
0579 int Nhits = 0;
0580
0581 if (theEvent.getByToken(cfToken, cf_simhit)) {
0582 std::unique_ptr<MixCollection<PSimHit> > thisContainerHits(new MixCollection<PSimHit>(cf_simhit.product()));
0583 for (auto const& isim : *thisContainerHits) {
0584 DetId theDet(isim.detUnitId());
0585 edm::EDConsumerBase::Labels labels;
0586 theEvent.labelsForToken(cfToken, labels);
0587 std::string trackerContainer(labels.productInstance);
0588 if (printOut_ && Nhits == 0)
0589 std::cout << " trackerContainer " << trackerContainer << std::endl;
0590 unsigned int tofBin = StripDigiSimLink::LowTof;
0591 if (trackerContainer.find(std::string("HighTof")) != std::string::npos)
0592 tofBin = StripDigiSimLink::HighTof;
0593 simHitCollectionID theSimHitCollID = std::make_pair(theDet.subdetId(), tofBin);
0594 SimHitCollMap_[theSimHitCollID].push_back(isim);
0595 ++Nhits;
0596 }
0597 if (printOut_ > 0)
0598 std::cout << "simHits from crossing frames; map size = " << SimHitCollMap_.size() << ", Hit count = " << Nhits
0599 << ", " << sizeof(SimHitCollMap_) << " bytes" << std::endl;
0600 }
0601 }
0602 for (auto const& simHitToken : simHitTokens_) {
0603 edm::Handle<std::vector<PSimHit> > simHits;
0604 int Nhits = 0;
0605 if (theEvent.getByToken(simHitToken, simHits)) {
0606 for (auto const& isim : *simHits) {
0607 DetId theDet(isim.detUnitId());
0608 edm::EDConsumerBase::Labels labels;
0609 theEvent.labelsForToken(simHitToken, labels);
0610 std::string trackerContainer(labels.productInstance);
0611 if (printOut_ && Nhits == 0)
0612 std::cout << " trackerContainer " << trackerContainer << std::endl;
0613 unsigned int tofBin = StripDigiSimLink::LowTof;
0614 if (trackerContainer.find(std::string("HighTof")) != std::string::npos)
0615 tofBin = StripDigiSimLink::HighTof;
0616 simHitCollectionID theSimHitCollID = std::make_pair(theDet.subdetId(), tofBin);
0617 SimHitCollMap_[theSimHitCollID].push_back(isim);
0618 ++Nhits;
0619 }
0620 if (printOut_ > 0)
0621 std::cout << "simHits from hard-scatter collection; map size = " << SimHitCollMap_.size()
0622 << ", Hit count = " << Nhits << ", " << sizeof(SimHitCollMap_) << " bytes" << std::endl;
0623 }
0624 }
0625 }
0626
0627
0628 DEFINE_FWK_MODULE(StripClusterMCanalysis);
