Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-16 00:33:31

0001 // -*- C++ -*-
0002 //
0003 // Package:    StripClusterMCanalysis
0004 // Class:      StripClusterMCanalysis
0005 //
0006 /**\class StripClusterMCanalysis StripClusterMCanalysis.cc UserCode/WTFord/detAnalysis/StripClusterMCanalysis/src/StripClusterMCanalysis.cc
0007 
0008  Description:  Analyze siStrip clusters for multiple track crossings
0009 
0010  The original purpose was to explore the separation via dE/dx between clusters
0011  originating from a single minimum-ionizing charged particle and those from
0012  overlap of multiple particles.  To this end we create a flat root ntuple 
0013  providing for each cluster information extracted from associated simTracks and 
0014  simHits regarding the charge, path length, particle ID, originating physics
0015  process, etc.
0016 
0017  Implementation:
0018      <Notes on implementation>
0019 */
0020 //
0021 //   Original Author:  William T. Ford
0022 //           Created:  Sat Nov 21 18:02:42 MST 2009
0023 // Adapted for CMSSW:  September, 2015
0024 // $Id: StripClusterMCanalysis.cc,v 1.3 2011/06/01 23:22:26 wtford Exp $
0025 //
0026 //
0027 
0028 // this class header
0029 
0030 // system include files
0031 #include <memory>
0032 
0033 // user include files
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 // Data Formats
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 //  Particle interaction process codes are found in
0073 // SimG4Core/Physics/src/G4ProcessTypeEnumerator.cc
0074 // See also https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMCTruth
0075 //
0076 
0077 //
0078 // class declaration
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   // ----------member data ---------------------------
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_;  // tree filled for each cluster
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_;  // tree filled for each pixel vertex
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 // constructors and destructor
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   //now do whatever initialization is needed
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 // member functions
0213 //
0214 
0215 // ------------ method called to for each event  ------------
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   // Retrieve tracker topology from geometry
0226   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0227 
0228   // Get the beam spot
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   // and the primary vertex
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     //    iterate over pixel vertices and fill the pixel vertex Ntuple
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   // Get the simHit info
0281   makeMap(iEvent);
0282 
0283   /////////////////////////
0284   // Cluster collections //
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   // Loop over subdetectors
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     // Find the path length of a straight track from the primary vertex
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     // Traverse the clusters for this subdetector
0336     for (edmNew::DetSet<SiStripCluster>::const_iterator ClusIter = DSViter->begin(); ClusIter != DSViter->end();
0337          ClusIter++) {
0338       clusterCount++;
0339 
0340       // Look for a digisimlink matching this cluster
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           // This digisimlink points to a strip in the current cluster
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           // Find simTracks associated with this cluster
0410           //
0411           unsigned int thisTrackID = linkiter->SimTrackId();
0412           // Does at least one strip have >1 track?
0413           if (stripIdx == prevIdx)
0414             ovlap = 1;
0415           // Have we seen this track yet?
0416           bool newTrack = true;
0417           if (std::find(trackID.begin(), trackID.end(), thisTrackID) != trackID.end())
0418             newTrack = false;
0419           if (newTrack) {
0420             // This is the first time we've encountered this track (linked to this cluster)
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           }  // if newTrack ... else
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           // Find simHits associaed with this cluster
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();  // Add up the contributions of all simHits to this cluster
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           }  // if (newHit)
0482           prevIdx = stripIdx;
0483         }  // DigiSimLink belongs to this cluster
0484       }    // Traverse DigiSimLinks
0485       if (prevIdx == -1)
0486         continue;  // No truth information for this cluster; move on to the next one.
0487 
0488       if (ovlap && trackID[1] < trackID[0])
0489         tkFlip = 1;
0490 
0491       // RecoTracker/DeDx/python/dedxDiscriminator_Prod_cfi.py, line 12 -- MeVperADCStrip = cms.double(3.61e-06*250)
0492 
0493       // Fill the cluster Ntuple
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     }  // traverse clusters in subdetector
0537   }    // traverse subdetectors
0538   if (printOut_ > 0)
0539     std::cout << clusterCount << " total clusters; " << clustersNoDigiSimLink << " without digiSimLinks" << std::endl;
0540 }
0541 
0542 // ------------ method called once each job just before starting event loop  ------------
0543 void StripClusterMCanalysis::beginJob() {
0544   int bufsize = 64000;
0545   edm::Service<TFileService> fs;
0546   // Create the cluster tree
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   // Create the vertex tree
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   // Book some histograms
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   //  The simHit collections are specified in trackerContainers, and can
0573   //  be either crossing frames (e.g., mix/g4SimHitsTrackerHitsTIBLowTof)
0574   //  or just PSimHits (e.g., g4SimHits/TrackerHitsTIBLowTof)
0575 
0576   SimHitCollMap_.clear();
0577   for (auto const& cfToken : cfTokens_) {
0578     edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
0579     int Nhits = 0;
0580     // if (theEvent.getByToken(cfToken, cf_simhit)) {
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 //define this as a plug-in
0628 DEFINE_FWK_MODULE(StripClusterMCanalysis);