Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-01-17 23:31:49 

0001 // File: ReadPixClusters.cc
0002 // Description: TO test the pixel clusters with tracks (full)
0003 // Author: Danek Kotlinski
0004 // Creation Date:  Initial version. 3/06
0005 //
0006 //--------------------------------------------
0007 #include <memory>
0008 #include <string>
0009 #include <iostream>
0010 
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017 
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 //#include "DataFormats/Common/interface/Handle.h"
0021 
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 
0026 //#include "DataFormats/SiPixelCluster/interface/SiPixelClusterCollection.h"
0027 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0028 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0029 #include "DataFormats/Common/interface/DetSetVector.h"
0030 #include "DataFormats/Common/interface/Ref.h"
0031 #include "DataFormats/DetId/interface/DetId.h"
0032 
0033 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0034 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0035 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0036 
0037 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0038 
0039 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0040 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0043 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 //#include "Geometry/CommonTopologies/interface/PixelTopology.h"
0046 
0047 // For L1
0048 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0049 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0050 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0051 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0052 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0053 
0054 // For HLT
0055 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0056 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0057 #include "DataFormats/Common/interface/TriggerResults.h"
0058 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0059 #include "FWCore/Common/interface/TriggerNames.h"
0060 
0061 // For tracks
0062 #include "DataFormats/TrackReco/interface/Track.h"
0063 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0064 
0065 //#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0066 
0067 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0068 
0069 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0070 //#include "TrackingTools/PatternTools/interface/TrajectoryFitter.h"
0071 
0072 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0073 
0074 //#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0075 //#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0076 
0077 //#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
0078 
0079 //#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0080 //#include "TrackingTools/PatternTools/interface/Trajectory.h"
0081 //#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0082 //#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0083 //#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0084 
0085 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0086 #include <DataFormats/VertexReco/interface/VertexFwd.h>
0087 
0088 // For luminisoty
0089 #include "DataFormats/Luminosity/interface/LumiSummary.h"
0090 #include "DataFormats/Common/interface/ConditionsInEdm.h"
0091 
0092 // To use root histos
0093 #include "FWCore/ServiceRegistry/interface/Service.h"
0094 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0095 
0096 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0097 
0098 // For ROOT
0099 #include <TROOT.h>
0100 #include <TChain.h>
0101 #include <TFile.h>
0102 #include <TF1.h>
0103 #include <TH2F.h>
0104 #include <TH1F.h>
0105 #include <TProfile.h>
0106 #include <TVector3.h>
0107 
0108 #define HISTOS
0109 //#define L1
0110 //#define HLT
0111 
0112 #define VDM_STUDIES
0113 
0114 const bool isData = false;  // set false for MC
0115 
0116 using namespace std;
0117 
0118 class TestWithTracks : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0119 public:
0120   explicit TestWithTracks(const edm::ParameterSet &conf);
0121   virtual ~TestWithTracks();
0122   virtual void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0123   virtual void beginJob() override;
0124   virtual void endJob() override;
0125 
0126 private:
0127   edm::EDGetTokenT<LumiSummary> lumiToken_;
0128   edm::EDGetTokenT<edm::ConditionsInLumiBlock> condToken_;
0129   edm::EDGetTokenT<L1GlobalTriggerReadoutRecord> l1gtrrToken_;
0130   edm::EDGetTokenT<edm::TriggerResults> hltToken_;
0131   edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0132   edm::EDGetTokenT<reco::TrackCollection> srcToken_;
0133   edm::EDGetTokenT<TrajTrackAssociationCollection> trackAssocToken_;
0134 
0135   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0136   //const static bool PRINT = false;
0137   bool PRINT;
0138   float countTracks, countGoodTracks, countTracksInPix, countPVs, countEvents, countLumi;
0139 
0140   TH1D *hcharge1, *hcharge2, *hcharge3, *hcharge4, *hcharge5;  // ADD FPIX
0141   TH1D *hsize1, *hsize2, *hsize3, *hsizex1, *hsizex2, *hsizex3, *hsizey1, *hsizey2, *hsizey3;
0142   TH1D *hsize4, *hsize5,  // ADD FPIX
0143       *hsizex4, *hsizex5, *hsizey4, *hsizey5;
0144 
0145   TH1D *hclusPerTrk1, *hclusPerTrk2, *hclusPerTrk3, *hclusPerTrk4, *hclusPerTrk5;
0146   TH1D *hclusPerLay1, *hclusPerLay2, *hclusPerLay3;
0147   TH1D *hclusPerDisk1, *hclusPerDisk2, *hclusPerDisk3, *hclusPerDisk4;
0148   TH1D *hclusPerTrk, *hclusPerTrkB, *hclusPerTrkF;
0149 
0150   TH2F *hDetMap1, *hDetMap2, *hDetMap3;           // clusters
0151   TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3;  // MODULE PROJECTION
0152 
0153   TH2F *hpvxy, *hclusMap1, *hclusMap2, *hclusMap3;  // Z vs PHI
0154 
0155   TH1D *hpvz, *hpvr, *hNumPv, *hNumPvClean;
0156   TH1D *hPt, *hEta, *hDz, *hD0, *hzdiff;
0157 
0158   TH1D *hl1a, *hl1t, *hlt1;
0159 
0160   TH1D *hclusBpix, *hpixBpix;
0161   TH1D *htracks, *htracksGood, *htracksGoodInPix;
0162 
0163   TProfile *hclumult1, *hclumult2, *hclumult3;
0164   TProfile *hclumultx1, *hclumultx2, *hclumultx3;
0165   TProfile *hclumulty1, *hclumulty2, *hclumulty3;
0166   TProfile *hcluchar1, *hcluchar2, *hcluchar3;
0167 
0168   TProfile *hclumultld1, *hclumultld2, *hclumultld3;
0169   TProfile *hclumultxld1, *hclumultxld2, *hclumultxld3;
0170   TProfile *hclumultyld1, *hclumultyld2, *hclumultyld3;
0171   TProfile *hclucharld1, *hclucharld2, *hclucharld3;
0172 
0173   TProfile *htracksls, *hpvsls, *htrackslsn, *hpvslsn, *hintgl, *hinstl, *hbeam1, *hbeam2;
0174   TProfile *hmult1, *hmult2, *hmult3;
0175   TProfile *hclusPerTrkVsEta, *hclusPerTrkVsPt, *hclusPerTrkVsls, *hclusPerTrkVsEtaB, *hclusPerTrkVsEtaF;
0176 
0177   TH1D *hlumi, *hlumi0, *hbx, *hbx0;
0178 
0179   TH1D *recHitXError1, *recHitXError2, *recHitXError3;
0180   TH1D *recHitYError1, *recHitYError2, *recHitYError3;
0181   TH1D *recHitXAlignError1, *recHitXAlignError2, *recHitXAlignError3;
0182   TH1D *recHitYAlignError1, *recHitYAlignError2, *recHitYAlignError3;
0183   TH1D *recHitXError4, *recHitXError5, *recHitXError6, *recHitXError7;
0184   TH1D *recHitYError4, *recHitYError5, *recHitYError6, *recHitYError7;
0185   TH1D *recHitXAlignError4, *recHitXAlignError5, *recHitXAlignError6, *recHitXAlignError7;
0186   TH1D *recHitYAlignError4, *recHitYAlignError5, *recHitYAlignError6, *recHitYAlignError7;
0187   TProfile *hErrorXB, *hErrorYB, *hErrorXF, *hErrorYF;
0188   TProfile *hAErrorXB, *hAErrorYB, *hAErrorXF, *hAErrorYF;
0189 
0190 #ifdef VDM_STUDIES
0191   TProfile *hcharCluls, *hcharPixls, *hsizeCluls, *hsizeXCluls;
0192   TProfile *hcharCluls1, *hcharPixls1, *hsizeCluls1, *hsizeXCluls1;
0193   TProfile *hcharCluls2, *hcharPixls2, *hsizeCluls2, *hsizeXCluls2;
0194   TProfile *hcharCluls3, *hcharPixls3, *hsizeCluls3, *hsizeXCluls3;
0195   TProfile *hclusls;  //   *hpixls;
0196   TProfile *hclusls1, *hclusls2, *hclusls3;
0197   TProfile *hclubx, *hpvbx, *htrackbx;  // *hcharClubx, *hcharPixbx,*hsizeClubx, *hsizeYClubx;
0198 #endif
0199 };
0200 /////////////////////////////////////////////////////////////////
0201 // Contructor,
0202 TestWithTracks::TestWithTracks(edm::ParameterSet const &conf) {
0203   usesResource(TFileService::kSharedResource);
0204   PRINT = conf.getUntrackedParameter<bool>("Verbosity", false);
0205   lumiToken_ = consumes<LumiSummary>(edm::InputTag("lumiProducer"));
0206   condToken_ = consumes<edm::ConditionsInLumiBlock>(edm::InputTag("conditionsInEdm"));
0207   l1gtrrToken_ = consumes<L1GlobalTriggerReadoutRecord>(edm::InputTag("gtDigis"));
0208   hltToken_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", "HLT"));
0209   vtxToken_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0210   srcToken_ = consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("src"));
0211   trackAssocToken_ =
0212       consumes<TrajTrackAssociationCollection>(edm::InputTag(conf.getParameter<std::string>("trajectoryInput")));
0213   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0214   //if(PRINT) edm::LogPrint("TestWithTracks")<<" Construct "<<endl;
0215 }
0216 
0217 // Virtual destructor needed.
0218 TestWithTracks::~TestWithTracks() = default;
0219 
0220 // ------------ method called at the begining   ------------
0221 void TestWithTracks::beginJob() {
0222   edm::LogPrint("TestWithTracks") << "BeginJob, Verbosity " << PRINT;
0223 
0224   countTracks = 0.;
0225   countGoodTracks = 0.;
0226   countTracksInPix = 0.;
0227   countPVs = 0.;
0228   countEvents = 0.;
0229   countLumi = 0.;
0230 
0231 #ifdef HISTOS
0232   edm::Service<TFileService> fs;
0233 
0234   int sizeH = 20;
0235   float lowH = -0.5;
0236   float highH = 19.5;
0237   hclusPerTrk1 = fs->make<TH1D>("hclusPerTrk1", "Clus per track l1", sizeH, lowH, highH);
0238   hclusPerTrk2 = fs->make<TH1D>("hclusPerTrk2", "Clus per track l2", sizeH, lowH, highH);
0239   hclusPerTrk3 = fs->make<TH1D>("hclusPerTrk3", "Clus per track l3", sizeH, lowH, highH);
0240   hclusPerTrk4 = fs->make<TH1D>("hclusPerTrk4", "Clus per track d1", sizeH, lowH, highH);
0241   hclusPerTrk5 = fs->make<TH1D>("hclusPerTrk5", "Clus per track d2", sizeH, lowH, highH);
0242   hclusPerTrk = fs->make<TH1D>("hclusPerTrk", "Clus per track", sizeH, lowH, highH);
0243   hclusPerTrkB = fs->make<TH1D>("hclusPerTrkB", "B Clus per track", sizeH, lowH, highH);
0244   hclusPerTrkF = fs->make<TH1D>("hclusPerTrkF", "F Clus per track", sizeH, lowH, highH);
0245 
0246   sizeH = 2000;
0247   highH = 1999.5;
0248   hclusPerLay1 = fs->make<TH1D>("hclusPerLay1", "Clus per layer l1", sizeH, lowH, highH);
0249   hclusPerLay2 = fs->make<TH1D>("hclusPerLay2", "Clus per layer l2", sizeH, lowH, highH);
0250   hclusPerLay3 = fs->make<TH1D>("hclusPerLay3", "Clus per layer l3", sizeH, lowH, highH);
0251   hclusPerDisk1 = fs->make<TH1D>("hclusPerDisk1", "Clus per disk 1", sizeH, lowH, highH);
0252   hclusPerDisk2 = fs->make<TH1D>("hclusPerDisk2", "Clus per disk 2", sizeH, lowH, highH);
0253   hclusPerDisk3 = fs->make<TH1D>("hclusPerDisk3", "Clus per disk 3", sizeH, lowH, highH);
0254   hclusPerDisk4 = fs->make<TH1D>("hclusPerDisk4", "Clus per disk 4", sizeH, lowH, highH);
0255 
0256   hcharge1 = fs->make<TH1D>("hcharge1", "Clu charge l1", 400, 0., 200.);  //in ke
0257   hcharge2 = fs->make<TH1D>("hcharge2", "Clu charge l2", 400, 0., 200.);
0258   hcharge3 = fs->make<TH1D>("hcharge3", "Clu charge l3", 400, 0., 200.);
0259   hcharge4 = fs->make<TH1D>("hcharge4", "Clu charge d1", 400, 0., 200.);
0260   hcharge5 = fs->make<TH1D>("hcharge5", "Clu charge d2", 400, 0., 200.);
0261 
0262   hsize1 = fs->make<TH1D>("hsize1", "layer 1 clu size", 100, -0.5, 99.5);
0263   hsize2 = fs->make<TH1D>("hsize2", "layer 2 clu size", 100, -0.5, 99.5);
0264   hsize3 = fs->make<TH1D>("hsize3", "layer 3 clu size", 100, -0.5, 99.5);
0265   hsizex1 = fs->make<TH1D>("hsizex1", "lay1 clu size in x", 20, -0.5, 19.5);
0266   hsizex2 = fs->make<TH1D>("hsizex2", "lay2 clu size in x", 20, -0.5, 19.5);
0267   hsizex3 = fs->make<TH1D>("hsizex3", "lay3 clu size in x", 20, -0.5, 19.5);
0268   hsizey1 = fs->make<TH1D>("hsizey1", "lay1 clu size in y", 30, -0.5, 29.5);
0269   hsizey2 = fs->make<TH1D>("hsizey2", "lay2 clu size in y", 30, -0.5, 29.5);
0270   hsizey3 = fs->make<TH1D>("hsizey3", "lay3 clu size in y", 30, -0.5, 29.5);
0271 
0272   hsize4 = fs->make<TH1D>("hsize4", "disk 1 clu size", 100, -0.5, 99.5);
0273   hsize5 = fs->make<TH1D>("hsize5", "disk 2 clu size", 100, -0.5, 99.5);
0274   hsizex4 = fs->make<TH1D>("hsizex4", "d1 clu size in x", 20, -0.5, 19.5);
0275   hsizex5 = fs->make<TH1D>("hsizex5", "d2 clu size in x", 20, -0.5, 19.5);
0276   hsizey4 = fs->make<TH1D>("hsizey4", "d1 clu size in y", 30, -0.5, 29.5);
0277   hsizey5 = fs->make<TH1D>("hsizey5", "d2 clu size in y", 30, -0.5, 29.5);
0278 
0279   hDetMap1 = fs->make<TH2F>("hDetMap1", "layer 1 clus map", 9, 0., 9., 23, 0., 23.);
0280   hDetMap2 = fs->make<TH2F>("hDetMap2", "layer 2 clus map", 9, 0., 9., 33, 0., 33.);
0281   hDetMap3 = fs->make<TH2F>("hDetMap3", "layer 3 clus map", 9, 0., 9., 45, 0., 45.);
0282   //hpixDetMap1 = fs->make<TH2F>( "hpixDetMap1", "pix det layer 1",
0283   //          416,0.,416.,160,0.,160.);
0284   //hpixDetMap2 = fs->make<TH2F>( "hpixDetMap2", "pix det layer 2",
0285   //          416,0.,416.,160,0.,160.);
0286   //hpixDetMap3 = fs->make<TH2F>( "hpixDetMap3", "pix det layer 3",
0287   //          416,0.,416.,160,0.,160.);
0288   hcluDetMap1 = fs->make<TH2F>("hcluDetMap1", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0289   hcluDetMap2 = fs->make<TH2F>("hcluDetMap2", "clu det layer 2", 416, 0., 416., 160, 0., 160.);
0290   hcluDetMap3 = fs->make<TH2F>("hcluDetMap3", "clu det layer 3", 416, 0., 416., 160, 0., 160.);
0291 
0292   htracksGoodInPix = fs->make<TH1D>("htracksGoodInPix", "count good tracks in pix", 2000, -0.5, 1999.5);
0293   htracksGood = fs->make<TH1D>("htracksGood", "count good tracks", 2000, -0.5, 1999.5);
0294   htracks = fs->make<TH1D>("htracks", "count tracks", 2000, -0.5, 1999.5);
0295   hclusBpix = fs->make<TH1D>("hclusBpix", "count clus in bpix", 200, -0.5, 1999.5);
0296   hpixBpix = fs->make<TH1D>("hpixBpix", "count pixels", 200, -0.5, 1999.5);
0297 
0298   hpvxy = fs->make<TH2F>("hpvxy", "pv xy", 100, -1., 1., 100, -1., 1.);
0299   hpvz = fs->make<TH1D>("hpvz", "pv z", 1000, -50., 50.);
0300   hpvr = fs->make<TH1D>("hpvr", "pv r", 100, 0., 1.);
0301   hNumPv = fs->make<TH1D>("hNumPv", "num of pv", 100, 0., 100.);
0302   hNumPvClean = fs->make<TH1D>("hNumPvClean", "num of pv clean", 100, 0., 100.);
0303 
0304   hPt = fs->make<TH1D>("hPt", "pt", 120, 0., 120.);
0305   hEta = fs->make<TH1D>("hEta", "eta", 50, -2.5, 2.5);
0306   hD0 = fs->make<TH1D>("hD0", "d0", 500, 0., 5.);
0307   hDz = fs->make<TH1D>("hDz", "pt", 250, -25., 25.);
0308   hzdiff = fs->make<TH1D>("hzdiff", "PVz-Trackz", 200, -10., 10.);
0309 
0310   hl1a = fs->make<TH1D>("hl1a", "l1a", 128, -0.5, 127.5);
0311   hl1t = fs->make<TH1D>("hl1t", "l1t", 128, -0.5, 127.5);
0312   hlt1 = fs->make<TH1D>("hlt1", "hlt1", 256, -0.5, 255.5);
0313 
0314   hclumult1 = fs->make<TProfile>("hclumult1", "cluster size layer 1", 60, -3., 3., 0.0, 100.);
0315   hclumult2 = fs->make<TProfile>("hclumult2", "cluster size layer 2", 60, -3., 3., 0.0, 100.);
0316   hclumult3 = fs->make<TProfile>("hclumult3", "cluster size layer 3", 60, -3., 3., 0.0, 100.);
0317 
0318   hclumultx1 = fs->make<TProfile>("hclumultx1", "cluster x-size layer 1", 60, -3., 3., 0.0, 100.);
0319   hclumultx2 = fs->make<TProfile>("hclumultx2", "cluster x-size layer 2", 60, -3., 3., 0.0, 100.);
0320   hclumultx3 = fs->make<TProfile>("hclumultx3", "cluster x-size layer 3", 60, -3., 3., 0.0, 100.);
0321 
0322   hclumulty1 = fs->make<TProfile>("hclumulty1", "cluster y-size layer 1", 60, -3., 3., 0.0, 100.);
0323   hclumulty2 = fs->make<TProfile>("hclumulty2", "cluster y-size layer 2", 60, -3., 3., 0.0, 100.);
0324   hclumulty3 = fs->make<TProfile>("hclumulty3", "cluster y-size layer 3", 60, -3., 3., 0.0, 100.);
0325 
0326   hcluchar1 = fs->make<TProfile>("hcluchar1", "cluster char layer 1", 60, -3., 3., 0.0, 1000.);
0327   hcluchar2 = fs->make<TProfile>("hcluchar2", "cluster char layer 2", 60, -3., 3., 0.0, 1000.);
0328   hcluchar3 = fs->make<TProfile>("hcluchar3", "cluster char layer 3", 60, -3., 3., 0.0, 1000.);
0329 
0330   // profiles versus ladder
0331   hclumultld1 = fs->make<TProfile>("hclumultld1", "cluster size layer 1", 23, -11.5, 11.5, 0.0, 100.);
0332   hclumultld2 = fs->make<TProfile>("hclumultld2", "cluster size layer 2", 35, -17.5, 17.5, 0.0, 100.);
0333   hclumultld3 = fs->make<TProfile>("hclumultld3", "cluster size layer 3", 47, -23.5, 23.5, 0.0, 100.);
0334 
0335   hclumultxld1 = fs->make<TProfile>("hclumultxld1", "cluster x-size layer 1", 23, -11.5, 11.5, 0.0, 100.);
0336   hclumultxld2 = fs->make<TProfile>("hclumultxld2", "cluster x-size layer 2", 35, -17.5, 17.5, 0.0, 100.);
0337   hclumultxld3 = fs->make<TProfile>("hclumultxld3", "cluster x-size layer 3", 47, -23.5, 23.5, 0.0, 100.);
0338 
0339   hclumultyld1 = fs->make<TProfile>("hclumultyld1", "cluster y-size layer 1", 23, -11.5, 11.5, 0.0, 100.);
0340   hclumultyld2 = fs->make<TProfile>("hclumultyld2", "cluster y-size layer 2", 35, -17.5, 17.5, 0.0, 100.);
0341   hclumultyld3 = fs->make<TProfile>("hclumultyld3", "cluster y-size layer 3", 47, -23.5, 23.5, 0.0, 100.);
0342 
0343   hclucharld1 = fs->make<TProfile>("hclucharld1", "cluster char layer 1", 23, -11.5, 11.5, 0.0, 1000.);
0344   hclucharld2 = fs->make<TProfile>("hclucharld2", "cluster char layer 2", 35, -17.5, 17.5, 0.0, 1000.);
0345   hclucharld3 = fs->make<TProfile>("hclucharld3", "cluster char layer 3", 47, -23.5, 23.5, 0.0, 1000.);
0346 
0347   hintgl = fs->make<TProfile>("hintgl", "inst lumi vs ls ", 1000, 0., 3000., 0.0, 10000.);
0348   hinstl = fs->make<TProfile>("hinstl", "intg lumi vs ls ", 1000, 0., 3000., 0.0, 100.);
0349   hbeam1 = fs->make<TProfile>("hbeam1", "beam1 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0350   hbeam2 = fs->make<TProfile>("hbeam2", "beam2 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0351 
0352   htracksls = fs->make<TProfile>("htracksls", "tracks with pix hits  vs ls", 1000, 0., 3000., 0.0, 10000.);
0353   hpvsls = fs->make<TProfile>("hpvsls", "pvs  vs ls", 1000, 0., 3000., 0.0, 1000.);
0354   htrackslsn = fs->make<TProfile>("htrackslsn", "tracks with pix hits/lumi  vs ls", 1000, 0., 3000., 0.0, 10000.);
0355   hpvslsn = fs->make<TProfile>("hpvslsn", "pvs/lumi  vs ls", 1000, 0., 3000., 0.0, 1000.);
0356 
0357   hmult1 = fs->make<TProfile>("hmult1", "clu mult layer 1", 10, 0., 10., 0.0, 1000.);
0358   hmult2 = fs->make<TProfile>("hmult2", "clu mult layer 2", 10, 0., 10., 0.0, 1000.);
0359   hmult3 = fs->make<TProfile>("hmult3", "clu mult layer 3", 10, 0., 10., 0.0, 1000.);
0360 
0361   hclusPerTrkVsEta = fs->make<TProfile>("hclusPerTrkVsEta", "clus per trk vs.eta", 60, -3., 3., 0.0, 100.);
0362   hclusPerTrkVsPt = fs->make<TProfile>("hclusPerTrkVsPt", "clus per trk vs.pt", 120, 0., 120., 0.0, 100.);
0363   hclusPerTrkVsls = fs->make<TProfile>("hclusPerTrkVsls", "clus per trk vs.ls", 300, 0., 3000., 0.0, 100.);
0364   hclusPerTrkVsEtaF = fs->make<TProfile>("hclusPerTrkVsEtaF", "F clus per trk vs.eta", 60, -3., 3., 0.0, 100.);
0365   hclusPerTrkVsEtaB = fs->make<TProfile>("hclusPerTrkVsEtaB", "B clus per trk vs.eta", 60, -3., 3., 0.0, 100.);
0366 
0367   hlumi0 = fs->make<TH1D>("hlumi0", "lumi", 2000, 0, 2000.);
0368   hlumi = fs->make<TH1D>("hlumi", "lumi", 2000, 0, 2000.);
0369   hbx0 = fs->make<TH1D>("hbx0", "bx", 4000, 0, 4000.);
0370   hbx = fs->make<TH1D>("hbx", "bx", 4000, 0, 4000.);
0371 
0372   hclusMap1 = fs->make<TH2F>("hclusMap1", "clus - lay1", 260, -26., 26., 350, -3.5, 3.5);
0373   hclusMap2 = fs->make<TH2F>("hclusMap2", "clus - lay2", 260, -26., 26., 350, -3.5, 3.5);
0374   hclusMap3 = fs->make<TH2F>("hclusMap3", "clus - lay3", 260, -26., 26., 350, -3.5, 3.5);
0375 
0376   // RecHit errors
0377   // alignment errors
0378   recHitXAlignError1 = fs->make<TH1D>("recHitXAlignError1", "RecHit X Alignment errors bpix 1", 100, 0., 100.);
0379   recHitYAlignError1 = fs->make<TH1D>("recHitYAlignError1", "RecHit Y Alignment errors bpix 1", 100, 0., 100.);
0380   recHitXAlignError2 = fs->make<TH1D>("recHitXAlignError2", "RecHit X Alignment errors bpix 2", 100, 0., 100.);
0381   recHitYAlignError2 = fs->make<TH1D>("recHitYAlignError2", "RecHit Y Alignment errors bpix 2", 100, 0., 100.);
0382   recHitXAlignError3 = fs->make<TH1D>("recHitXAlignError3", "RecHit X Alignment errors bpix 3", 100, 0., 100.);
0383   recHitYAlignError3 = fs->make<TH1D>("recHitYAlignError3", "RecHit Y Alignment errors bpix 3", 100, 0., 100.);
0384   recHitXAlignError4 = fs->make<TH1D>("recHitXAlignError4", "RecHit X Alignment errors fpix -d2", 100, 0., 100.);
0385   recHitYAlignError4 = fs->make<TH1D>("recHitYAlignError4", "RecHit Y Alignment errors fpix -d2", 100, 0., 100.);
0386   recHitXAlignError5 = fs->make<TH1D>("recHitXAlignError5", "RecHit X Alignment errors fpix -d1", 100, 0., 100.);
0387   recHitYAlignError5 = fs->make<TH1D>("recHitYAlignError5", "RecHit Y Alignment errors fpix -d1", 100, 0., 100.);
0388   recHitXAlignError6 = fs->make<TH1D>("recHitXAlignError6", "RecHit X Alignment errors fpix +d1", 100, 0., 100.);
0389   recHitYAlignError6 = fs->make<TH1D>("recHitYAlignError6", "RecHit Y Alignment errors fpix +d1", 100, 0., 100.);
0390   recHitXAlignError7 = fs->make<TH1D>("recHitXAlignError7", "RecHit X Alignment errors fpix +d2", 100, 0., 100.);
0391   recHitYAlignError7 = fs->make<TH1D>("recHitYAlignError7", "RecHit Y Alignment errors fpix +d2", 100, 0., 100.);
0392 
0393   recHitXError1 = fs->make<TH1D>("recHitXError1", "RecHit X errors bpix 1", 100, 0., 100.);
0394   recHitYError1 = fs->make<TH1D>("recHitYError1", "RecHit Y errors bpix 1", 100, 0., 100.);
0395   recHitXError2 = fs->make<TH1D>("recHitXError2", "RecHit X errors bpix 2", 100, 0., 100.);
0396   recHitYError2 = fs->make<TH1D>("recHitYError2", "RecHit Y errors bpix 2", 100, 0., 100.);
0397   recHitXError3 = fs->make<TH1D>("recHitXError3", "RecHit X errors bpix 3", 100, 0., 100.);
0398   recHitYError3 = fs->make<TH1D>("recHitYError3", "RecHit Y errors bpix 3", 100, 0., 100.);
0399   recHitXError4 = fs->make<TH1D>("recHitXError4", "RecHit X errors fpix -d2", 100, 0., 100.);
0400   recHitYError4 = fs->make<TH1D>("recHitYError4", "RecHit Y errors fpix -d2", 100, 0., 100.);
0401   recHitXError5 = fs->make<TH1D>("recHitXError5", "RecHit X errors fpix -d1", 100, 0., 100.);
0402   recHitYError5 = fs->make<TH1D>("recHitYError5", "RecHit Y errors fpix -d1", 100, 0., 100.);
0403   recHitXError6 = fs->make<TH1D>("recHitXError6", "RecHit X errors fpix +d1", 100, 0., 100.);
0404   recHitYError6 = fs->make<TH1D>("recHitYError6", "RecHit Y errors fpix +d1", 100, 0., 100.);
0405   recHitXError7 = fs->make<TH1D>("recHitXError7", "RecHit X errors fpix +d2", 100, 0., 100.);
0406   recHitYError7 = fs->make<TH1D>("recHitYError7", "RecHit Y errors fpix +d2", 100, 0., 100.);
0407 
0408   hErrorXB = fs->make<TProfile>("hErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.);
0409   hErrorXF = fs->make<TProfile>("hErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.);
0410   hErrorYB = fs->make<TProfile>("hErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.);
0411   hErrorYF = fs->make<TProfile>("hErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.);
0412 
0413   hAErrorXB = fs->make<TProfile>("hAErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.);
0414   hAErrorXF = fs->make<TProfile>("hAErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.);
0415   hAErrorYB = fs->make<TProfile>("hAErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.);
0416   hAErrorYF = fs->make<TProfile>("hAErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.);
0417 
0418 #ifdef VDM_STUDIES
0419   highH = 3000.;
0420   sizeH = 1000;
0421 
0422   hclusls = fs->make<TProfile>("hclusls", "clus vs ls", sizeH, 0., highH, 0.0, 30000.);
0423   //hpixls  = fs->make<TProfile>("hpixls", "pix vs ls ",sizeH,0.,highH,0.0,100000.);
0424 
0425   hcharCluls = fs->make<TProfile>("hcharCluls", "clu char vs ls", sizeH, 0., highH, 0.0, 100.);
0426   //hcharPixls = fs->make<TProfile>("hcharPixls","pix char vs ls",sizeH,0.,highH,0.0,100.);
0427   hsizeCluls = fs->make<TProfile>("hsizeCluls", "clu size vs ls", sizeH, 0., highH, 0.0, 1000.);
0428   hsizeXCluls = fs->make<TProfile>("hsizeXCluls", "clu size-x vs ls", sizeH, 0., highH, 0.0, 100.);
0429 
0430   hcharCluls1 = fs->make<TProfile>("hcharCluls1", "clu char vs ls", sizeH, 0., highH, 0.0, 100.);
0431   //hcharPixls1 = fs->make<TProfile>("hcharPixls1","pix char vs ls",sizeH,0.,highH,0.0,100.);
0432   hsizeCluls1 = fs->make<TProfile>("hsizeCluls1", "clu size vs ls", sizeH, 0., highH, 0.0, 1000.);
0433   hsizeXCluls1 = fs->make<TProfile>("hsizeXCluls1", "clu size-x vs ls", sizeH, 0., highH, 0.0, 100.);
0434   hcharCluls2 = fs->make<TProfile>("hcharCluls2", "clu char vs ls", sizeH, 0., highH, 0.0, 100.);
0435   //hcharPixls2 = fs->make<TProfile>("hcharPixls2","pix char vs ls",sizeH,0.,highH,0.0,100.);
0436   hsizeCluls2 = fs->make<TProfile>("hsizeCluls2", "clu size vs ls", sizeH, 0., highH, 0.0, 1000.);
0437   hsizeXCluls2 = fs->make<TProfile>("hsizeXCluls2", "clu size-x vs ls", sizeH, 0., highH, 0.0, 100.);
0438   hcharCluls3 = fs->make<TProfile>("hcharCluls3", "clu char vs ls", sizeH, 0., highH, 0.0, 100.);
0439   //hcharPixls3 = fs->make<TProfile>("hcharPixls3","pix char vs ls",sizeH,0.,highH,0.0,100.);
0440   hsizeCluls3 = fs->make<TProfile>("hsizeCluls3", "clu size vs ls", sizeH, 0., highH, 0.0, 1000.);
0441   hsizeXCluls3 = fs->make<TProfile>("hsizeXCluls3", "clu size-x vs ls", sizeH, 0., highH, 0.0, 100.);
0442   hclusls1 = fs->make<TProfile>("hclusls1", "clus vs ls", sizeH, 0., highH, 0.0, 30000.);
0443   //hpixls1  = fs->make<TProfile>("hpixls1", "pix vs ls ",sizeH,0.,highH,0.0,100000.);
0444   hclusls2 = fs->make<TProfile>("hclusls2", "clus vs ls", sizeH, 0., highH, 0.0, 30000.);
0445   //hpixls2  = fs->make<TProfile>("hpixls2", "pix vs ls ",sizeH,0.,highH,0.0,100000.);
0446   hclusls3 = fs->make<TProfile>("hclusls3", "clus vs ls", sizeH, 0., highH, 0.0, 30000.);
0447   //hpixls3  = fs->make<TProfile>("hpixls3", "pix vs ls ",sizeH,0.,highH,0.0,100000.);
0448 
0449   // Profiles versus bx
0450   //hpixbx  = fs->make<TProfile>("hpixbx", "pixs vs bx ",4000,-0.5,3999.5,0.0,1000000.);
0451   hclubx = fs->make<TProfile>("hclubx", "clus vs bx ", 4000, -0.5, 3999.5, 0.0, 1000000.);
0452   hpvbx = fs->make<TProfile>("hpvbx", "pvs vs bx ", 4000, -0.5, 3999.5, 0.0, 1000000.);
0453   htrackbx = fs->make<TProfile>("htrackbx", "tracks vs bx ", 4000, -0.5, 3999.5, 0.0, 1000000.);
0454 
0455 #endif
0456 
0457 #endif
0458 }
0459 // ------------ method called to at the end of the job  ------------
0460 void TestWithTracks::endJob() {
0461   edm::LogPrint("TestWithTracks") << " End PixelTracksTest, events =  " << countEvents;
0462 
0463   if (countEvents > 0.) {
0464     countTracks /= countEvents;
0465     countGoodTracks /= countEvents;
0466     countTracksInPix /= countEvents;
0467     countPVs /= countEvents;
0468     countLumi /= 1000.;
0469     edm::LogPrint("TestWithTracks") << " Average tracks/event " << countTracks << " good " << countGoodTracks
0470                                     << " in pix " << countTracksInPix << " PVs " << countPVs << " events "
0471                                     << countEvents << " lumi pb-1 " << countLumi << "/10, bug!";
0472   }
0473 }
0474 //////////////////////////////////////////////////////////////////
0475 // Functions that gets called by framework every event
0476 void TestWithTracks::analyze(const edm::Event &e, const edm::EventSetup &es) {
0477   using namespace edm;
0478   using namespace reco;
0479   static LuminosityBlockNumber_t lumiBlockOld = -9999;
0480 
0481   const float CLU_SIZE_PT_CUT = 1.;
0482 
0483   int trackNumber = 0;
0484   int countNiceTracks = 0;
0485   int countPixTracks = 0;
0486 
0487   int numberOfClusters = 0;
0488   int numberOfPixels = 0;
0489 
0490   int numOfClusPerTrk1 = 0;
0491   int numOfClustersPerLay1 = 0;
0492   int numOfPixelsPerLay1 = 0;
0493 
0494   int numOfClusPerTrk2 = 0;
0495   int numOfClustersPerLay2 = 0;
0496   int numOfPixelsPerLay2 = 0;
0497 
0498   int numOfClusPerTrk3 = 0;
0499   int numOfClustersPerLay3 = 0;
0500   int numOfPixelsPerLay3 = 0;
0501 
0502   int numOfClusPerTrk4 = 0;
0503   int numOfClustersPerLay4 = 0;
0504   int numOfPixelsPerLay4 = 0;
0505 
0506   int numOfClusPerTrk5 = 0;
0507   int numOfClustersPerLay5 = 0;
0508   int numOfPixelsPerLay5 = 0;
0509 
0510   int numOfClustersPerDisk1 = 0;
0511   int numOfClustersPerDisk2 = 0;
0512   int numOfClustersPerDisk3 = 0;
0513   int numOfClustersPerDisk4 = 0;
0514 
0515   RunNumber_t const run = e.id().run();
0516   EventNumber_t const event = e.id().event();
0517   LuminosityBlockNumber_t const lumiBlock = e.luminosityBlock();
0518 
0519   int bx = e.bunchCrossing();
0520   //int orbit     = e.orbitNumber(); // unused
0521 
0522   if (PRINT)
0523     edm::LogPrint("TestWithTracks") << "Run " << run << " Event " << event << " LS " << lumiBlock;
0524 
0525   hbx0->Fill(float(bx));
0526   hlumi0->Fill(float(lumiBlock));
0527 
0528   edm::LuminosityBlock const &iLumi = e.getLuminosityBlock();
0529   edm::Handle<LumiSummary> lumi;
0530   iLumi.getByToken(lumiToken_, lumi);
0531   edm::Handle<edm::ConditionsInLumiBlock> cond;
0532   float intlumi = 0, instlumi = 0;
0533   int beamint1 = 0, beamint2 = 0;
0534   iLumi.getByToken(condToken_, cond);
0535   // This will only work when running on RECO until (if) they fix it in the FW
0536   // When running on RAW and reconstructing, the LumiSummary will not appear
0537   // in the event before reaching endLuminosityBlock(). Therefore, it is not
0538   // possible to get this info in the event
0539   if (lumi.isValid()) {
0540     intlumi = (lumi->intgRecLumi()) / 1000.;  // 10^30 -> 10^33/cm2/sec  ->  1/nb/sec
0541     instlumi = (lumi->avgInsDelLumi()) / 1000.;
0542     beamint1 = (cond->totalIntensityBeam1) / 1000;
0543     beamint2 = (cond->totalIntensityBeam2) / 1000;
0544   } else {
0545     //std::edm::LogPrint("TestWithTracks") << "** ERROR: Event does not get lumi info\n";
0546   }
0547 
0548   hinstl->Fill(float(lumiBlock), float(instlumi));
0549   hintgl->Fill(float(lumiBlock), float(intlumi));
0550   hbeam1->Fill(float(lumiBlock), float(beamint1));
0551   hbeam2->Fill(float(lumiBlock), float(beamint2));
0552 
0553 #ifdef L1
0554   // Get L1
0555   Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
0556   e.getByToken(l1gtrrToken_, L1GTRR);
0557 
0558   if (L1GTRR.isValid()) {
0559     //bool l1a = L1GTRR->decision();  // global decission?
0560     //edm::LogPrint("TestWithTracks")<<" L1 status :"<<l1a<<" "<<hex;
0561     for (unsigned int i = 0; i < L1GTRR->decisionWord().size(); ++i) {
0562       int l1flag = L1GTRR->decisionWord()[i];
0563       int t1flag = L1GTRR->technicalTriggerWord()[i];
0564 
0565       if (l1flag > 0)
0566         hl1a->Fill(float(i));
0567       if (t1flag > 0 && i < 64)
0568         hl1t->Fill(float(i));
0569     }  // for loop
0570   }    // if l1a
0571 #endif
0572 
0573 #ifdef HLT
0574 
0575   bool hlt[256];
0576   for (int i = 0; i < 256; ++i)
0577     hlt[i] = false;
0578 
0579   edm::TriggerNames TrigNames;
0580   edm::Handle<edm::TriggerResults> HLTResults;
0581 
0582   // Extract the HLT results
0583   e.getByToken(hltToken_, HLTResults);
0584   if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) {
0585     //TrigNames.init(*HLTResults);
0586     const edm::TriggerNames &TrigNames = e.triggerNames(*HLTResults);
0587 
0588     //edm::LogPrint("TestWithTracks")<<TrigNames.triggerNames().size()<<endl;
0589 
0590     for (unsigned int i = 0; i < TrigNames.triggerNames().size(); i++) {  // loop over trigger
0591       //if(countAllEvents==1) edm::LogPrint("TestWithTracks")<<i<<" "<<TrigNames.triggerName(i)<<endl;
0592 
0593       if ((HLTResults->wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0594           (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0595           (HLTResults->error(TrigNames.triggerIndex(TrigNames.triggerName(i))) == false)) {
0596         hlt[i] = true;
0597         hlt1->Fill(float(i));
0598 
0599       }  // if hlt
0600 
0601     }  // loop
0602   }    // if valid
0603 #endif
0604 
0605   // Get event setup
0606   edm::ESHandle<TrackerGeometry> geom = es.getHandle(trackerGeomToken_);
0607   const TrackerGeometry &theTracker(*geom);
0608 
0609   // -- Primary vertices
0610   // ----------------------------------------------------------------------
0611   edm::Handle<reco::VertexCollection> vertices;
0612   e.getByToken(vtxToken_, vertices);
0613 
0614   if (PRINT)
0615     edm::LogPrint("TestWithTracks") << " PV list " << vertices->size();
0616   int pvNotFake = 0, pvsTrue = 0;
0617   vector<float> pvzVector;
0618   for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) {
0619     if ((iv->isFake()) == 1)
0620       continue;
0621     pvNotFake++;
0622     float pvx = iv->x();
0623     float pvy = iv->y();
0624     float pvz = iv->z();
0625     int numTracksPerPV = iv->tracksSize();
0626     //int numTracksPerPV = iv->nTracks();
0627 
0628     //float xe = iv->xError();
0629     //float ye = iv->yError();
0630     //float ze = iv->zError();
0631     //int chi2 = iv->chi2();
0632     //int dof = iv->ndof();
0633 
0634     if (PRINT)
0635       edm::LogPrint("TestWithTracks") << " PV " << pvNotFake << " pos = " << pvx << "/" << pvy << "/" << pvz
0636                                       << ", Num of tracks " << numTracksPerPV;
0637 
0638     hpvz->Fill(pvz);
0639     if (pvz > -22. && pvz < 22.) {
0640       float pvr = sqrt(pvx * pvx + pvy * pvy);
0641       hpvxy->Fill(pvx, pvy);
0642       hpvr->Fill(pvr);
0643       if (pvr < 0.3) {
0644         pvsTrue++;
0645         pvzVector.push_back(pvz);
0646         //if(PRINT) edm::LogPrint("TestWithTracks")<<"PV "<<pvsTrue<<" "<<pvz<<endl;
0647       }  //pvr
0648     }    // pvz
0649 
0650     //if(pvsTrue<1) continue; // skip events with no PV
0651 
0652   }  // loop pvs
0653   hNumPv->Fill(float(pvNotFake));
0654   hNumPvClean->Fill(float(pvsTrue));
0655 
0656   if (PRINT)
0657     edm::LogPrint("TestWithTracks") << " Not fake PVs = " << pvNotFake << " good position " << pvsTrue;
0658 
0659   // -- Tracks
0660   // ----------------------------------------------------------------------
0661   Handle<reco::TrackCollection> recTracks;
0662   e.getByToken(srcToken_, recTracks);
0663 
0664   if (PRINT)
0665     edm::LogPrint("TestWithTracks") << " Tracks " << recTracks->size();
0666   for (reco::TrackCollection::const_iterator t = recTracks->begin(); t != recTracks->end(); ++t) {
0667     trackNumber++;
0668     numOfClusPerTrk1 = 0;  // this is confusing, it is used as clus per track
0669     numOfClusPerTrk2 = 0;
0670     numOfClusPerTrk3 = 0;
0671     numOfClusPerTrk4 = 0;
0672     numOfClusPerTrk5 = 0;
0673     int pixelHits = 0;
0674 
0675     int size = t->recHitsSize();
0676     float pt = t->pt();
0677     float eta = t->eta();
0678     float phi = t->phi();
0679     //float trackCharge = t->charge(); // unused
0680     float d0 = t->d0();
0681     float dz = t->dz();
0682     //float tkvx = t->vx();  // unused
0683     //float tkvy = t->vy();
0684     //float tkvz = t->vz();
0685 
0686     if (PRINT)
0687       edm::LogPrint("TestWithTracks") << "Track " << trackNumber << " Pt " << pt << " Eta " << eta << " d0/dz " << d0
0688                                       << " " << dz << " Hits " << size;
0689 
0690     hEta->Fill(eta);
0691     hDz->Fill(dz);
0692     if (abs(eta) > 2.8 || abs(dz) > 25.)
0693       continue;  //  skip
0694 
0695     hD0->Fill(d0);
0696     if (d0 > 1.0)
0697       continue;  // skip
0698 
0699     bool goodTrack = false;
0700     for (vector<float>::iterator m = pvzVector.begin(); m != pvzVector.end(); ++m) {
0701       float z = *m;
0702       float tmp = abs(z - dz);
0703       hzdiff->Fill(tmp);
0704       if (tmp < 1.)
0705         goodTrack = true;
0706     }
0707 
0708     if (isData && !goodTrack)
0709       continue;
0710     countNiceTracks++;
0711     hPt->Fill(pt);
0712 
0713     // Loop over rechits
0714     for (trackingRecHit_iterator recHit = t->recHitsBegin(); recHit != t->recHitsEnd(); ++recHit) {
0715       if (!((*recHit)->isValid()))
0716         continue;
0717 
0718       if ((*recHit)->geographicalId().det() != DetId::Tracker)
0719         continue;
0720 
0721       const DetId &hit_detId = (*recHit)->geographicalId();
0722       uint IntSubDetID = (hit_detId.subdetId());
0723 
0724       // Select pixel rechits
0725       if (IntSubDetID == 0)
0726         continue;  // Select ??
0727 
0728       int layer = 0, ladder = 0, zindex = 0, ladderOn = 0, module = 0, shell = 0;
0729       unsigned int disk = 0;     //1,2,3
0730       unsigned int blade = 0;    //1-24
0731       unsigned int zindexF = 0;  //
0732       unsigned int side = 0;     //size=1 for -z, 2 for +z
0733       unsigned int panel = 0;    //panel=1
0734 
0735       if (IntSubDetID == PixelSubdetector::PixelBarrel) {  // bpix
0736 
0737         // Pixel detector
0738         PXBDetId pdetId = PXBDetId(hit_detId);
0739         //unsigned int detTypeP=pdetId.det();
0740         //unsigned int subidP=pdetId.subdetId();
0741         // Barell layer = 1,2,3
0742         layer = pdetId.layer();
0743         // Barrel ladder id 1-20,32,44.
0744         ladder = pdetId.ladder();
0745         // Barrel Z-index=1,8
0746         zindex = pdetId.module();
0747         if (zindex < 5)
0748           side = 1;
0749         else
0750           side = 2;
0751 
0752         // Convert to online
0753         PixelBarrelName pbn(pdetId);
0754         // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0755         PixelBarrelName::Shell sh = pbn.shell();  //enum
0756         //sector = pbn.sectorName();
0757         ladderOn = pbn.ladderName();
0758         //layerOn  = pbn.layerName();
0759         module = pbn.moduleName();  // 1 to 4
0760         //half  = pbn.isHalfModule();
0761         shell = int(sh);
0762         // change the module sign for z<0
0763         if (shell == 1 || shell == 2)
0764           module = -module;  // make -1 to -4 for -z
0765                              // change ladeer sign for Outer )x<0)
0766         if (shell == 1 || shell == 3)
0767           ladderOn = -ladderOn;
0768 
0769         if (PRINT)
0770           edm::LogPrint("TestWithTracks") << "barrel layer/ladder/module: " << layer << "/" << ladder << "/" << zindex;
0771 
0772       } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {  // fpix
0773 
0774         PXFDetId pdetId = PXFDetId(hit_detId);
0775         disk = pdetId.disk();       //1,2,3
0776         blade = pdetId.blade();     //1-24
0777         zindexF = pdetId.module();  //
0778         side = pdetId.side();       //size=1 for -z, 2 for +z
0779         panel = pdetId.panel();     //panel=1
0780 
0781         if (PRINT)
0782           edm::LogPrint("TestWithTracks") << " forward det, disk " << disk << ", blade " << blade << ", module "
0783                                           << zindexF << ", side " << side << ", panel " << panel;
0784 
0785       } else {     // nothings
0786         continue;  // skip
0787       }
0788 
0789       // Get the geom-detector
0790       const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
0791       //double detZ = theGeomDet->surface().position().z();  // unused
0792       //double detR = theGeomDet->surface().position().perp(); // unused
0793       const PixelTopology *topol = &(theGeomDet->specificTopology());  // pixel topology
0794 
0795       //std::vector<SiStripRecHit2D*> output = getRecHitComponents((*recHit).get());
0796       //std::vector<SiPixelRecHit*> TrkComparison::getRecHitComponents(const TrackingRecHit* rechit){
0797 
0798       const SiPixelRecHit *hit = dynamic_cast<const SiPixelRecHit *>((*recHit));
0799       //edm::Ref<edmNew::DetSetVector<SiStripCluster> ,SiStripCluster> cluster = hit->cluster();
0800       // get the edm::Ref to the cluster
0801 
0802       if (hit) {
0803         if (pt > 1.) {  // eliminate low pt tracks
0804           // RecHit (recthits are transient, so not available without track refit)
0805           double xloc = hit->localPosition().x();  // 1st meas coord
0806           double yloc = hit->localPosition().y();  // 2nd meas coord or zero
0807           //double zloc = hit->localPosition().z();// up, always zero
0808           LocalError lerr = hit->localPositionError();
0809           float lerr_x = sqrt(lerr.xx()) * 1E4;
0810           float lerr_y = sqrt(lerr.yy()) * 1E4;
0811 
0812           if (layer == 1) {
0813             recHitXError1->Fill(lerr_x);
0814             recHitYError1->Fill(lerr_y);
0815             hErrorXB->Fill(float(ladder + (110 * (side - 1))), lerr_x);
0816             hErrorYB->Fill(float(ladder + (110 * (side - 1))), lerr_y);
0817           } else if (layer == 2) {
0818             recHitXError2->Fill(lerr_x);
0819             recHitYError2->Fill(lerr_y);
0820             hErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), lerr_x);
0821             hErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), lerr_y);
0822 
0823           } else if (layer == 3) {
0824             recHitXError3->Fill(lerr_x);
0825             recHitYError3->Fill(lerr_y);
0826             hErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), lerr_x);
0827             hErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), lerr_y);
0828           } else if ((disk == 2) && (side == 1)) {
0829             recHitXError4->Fill(lerr_x);
0830             recHitYError4->Fill(lerr_y);
0831             hErrorXF->Fill(float(blade), lerr_x);
0832             hErrorYF->Fill(float(blade), lerr_y);
0833 
0834           } else if ((disk == 1) && (side == 1)) {
0835             recHitXError5->Fill(lerr_x);
0836             recHitYError5->Fill(lerr_y);
0837             hErrorXF->Fill(float(blade + 25), lerr_x);
0838             hErrorYF->Fill(float(blade + 25), lerr_y);
0839 
0840           } else if ((disk == 1) && (side == 2)) {
0841             recHitXError6->Fill(lerr_x);
0842             recHitYError6->Fill(lerr_y);
0843             hErrorXF->Fill(float(blade + 50), lerr_x);
0844             hErrorYF->Fill(float(blade + 50), lerr_y);
0845           } else if ((disk == 2) && (side == 2)) {
0846             recHitXError7->Fill(lerr_x);
0847             recHitYError7->Fill(lerr_y);
0848             hErrorXF->Fill(float(blade + 75), lerr_x);
0849             hErrorYF->Fill(float(blade + 75), lerr_y);
0850           }
0851 
0852           LocalError lape = theGeomDet->localAlignmentError();
0853           if (lape.valid()) {
0854             float tmp11 = 0.;
0855             if (lape.xx() > 0.)
0856               tmp11 = sqrt(lape.xx()) * 1E4;
0857             //float tmp12= sqrt(lape.xy())*1E4;
0858             float tmp13 = 0.;
0859             if (lape.yy() > 0.)
0860               tmp13 = sqrt(lape.yy()) * 1E4;
0861             //bool tmp14=tmp2<tmp1;
0862             if (layer == 1) {
0863               recHitXAlignError1->Fill(tmp11);
0864               recHitYAlignError1->Fill(tmp13);
0865               hAErrorXB->Fill(float(ladder + (110 * (side - 1))), tmp11);
0866               hAErrorYB->Fill(float(ladder + (110 * (side - 1))), tmp13);
0867             } else if (layer == 2) {
0868               recHitXAlignError2->Fill(tmp11);
0869               recHitYAlignError2->Fill(tmp13);
0870               hAErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), tmp11);
0871               hAErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), tmp13);
0872             } else if (layer == 3) {
0873               recHitXAlignError3->Fill(tmp11);
0874               recHitYAlignError3->Fill(tmp13);
0875               hAErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), tmp11);
0876               hAErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), tmp13);
0877 
0878             } else if ((disk == 2) && (side == 1)) {
0879               recHitXAlignError4->Fill(tmp11);
0880               recHitYAlignError4->Fill(tmp13);
0881               hAErrorXF->Fill(float(blade), tmp11);
0882               hAErrorYF->Fill(float(blade), tmp13);
0883 
0884             } else if ((disk == 1) && (side == 1)) {
0885               recHitXAlignError5->Fill(tmp11);
0886               recHitYAlignError5->Fill(tmp13);
0887               hAErrorXF->Fill(float(blade + 25), tmp11);
0888               hAErrorYF->Fill(float(blade + 25), tmp13);
0889             } else if ((disk == 1) && (side == 2)) {
0890               recHitXAlignError6->Fill(tmp11);
0891               recHitYAlignError6->Fill(tmp13);
0892               hAErrorXF->Fill(float(blade + 50), tmp11);
0893               hAErrorYF->Fill(float(blade + 50), tmp13);
0894             } else if ((disk == 2) && (side == 2)) {
0895               recHitXAlignError7->Fill(tmp11);
0896               recHitYAlignError7->Fill(tmp13);
0897               hAErrorXF->Fill(float(blade + 75), tmp11);
0898               hAErrorYF->Fill(float(blade + 75), tmp13);
0899             }
0900 
0901             //edm::LogPrint("TestWithTracks")<<tTopo->pxbLayer(detId)<<" "<<tTopo->pxbModule(detId)<<" "<<rows<<" "<<tmp14<<" "
0902             if (PRINT)
0903               edm::LogPrint("TestWithTracks") << " align error " << layer << tmp11 << " " << tmp13;
0904           } else {
0905             edm::LogPrint("TestWithTracks") << " lape = 0";
0906           }  // if lape
0907 
0908           if (PRINT)
0909             edm::LogPrint("TestWithTracks") << " rechit loc " << xloc << " " << yloc << " " << lerr_x << " " << lerr_y;
0910         }  // limit pt
0911 
0912         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = hit->cluster();
0913         //  check if the ref is not null
0914         if (!clust.isNonnull())
0915           continue;
0916 
0917         numberOfClusters++;
0918         pixelHits++;
0919         float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
0920         int size = clust->size();
0921         int sizeX = clust->sizeX();
0922         int sizeY = clust->sizeY();
0923         float row = clust->x();
0924         float col = clust->y();
0925         numberOfPixels += size;
0926 
0927         //edm::LogPrint("TestWithTracks")<<" clus loc "<<row<<" "<<col<<endl;
0928 
0929         if (PRINT)
0930           edm::LogPrint("TestWithTracks")
0931               << " cluster " << numberOfClusters << " charge = " << charge << " size = " << size;
0932 
0933         LocalPoint lp = topol->localPosition(MeasurementPoint(clust->x(), clust->y()));
0934         //float x = lp.x();
0935         //float y = lp.y();
0936         //edm::LogPrint("TestWithTracks")<<" clu loc "<<x<<" "<<y<<endl;
0937 
0938         GlobalPoint clustgp = theGeomDet->surface().toGlobal(lp);
0939         double gX = clustgp.x();
0940         double gY = clustgp.y();
0941         double gZ = clustgp.z();
0942 
0943         //edm::LogPrint("TestWithTracks")<<" CLU GLOBAL "<<gX<<" "<<gY<<" "<<gZ<<endl;
0944 
0945         TVector3 v(gX, gY, gZ);
0946         //float phi = v.Phi(); // unused
0947 
0948         //int maxPixelCol = clust->maxPixelCol();
0949         //int maxPixelRow = clust->maxPixelRow();
0950         //int minPixelCol = clust->minPixelCol();
0951         //int minPixelRow = clust->minPixelRow();
0952         //int geoId = PixGeom->geographicalId().rawId();
0953         // Replace with the topology methods
0954         // edge method moved to topologi class
0955         //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
0956         //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
0957 
0958         // calculate alpha and beta from cluster position
0959         //LocalTrajectoryParameters ltp = tsos.localParameters();
0960         //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0961         //float locx = localDir.x();
0962         //float locy = localDir.y();
0963         //float locz = localDir.z();
0964         //float loctheta = localDir.theta(); // currently unused
0965         //float alpha = atan2( locz, locx );
0966         //float beta = atan2( locz, locy );
0967 
0968         if (layer == 1) {
0969           hDetMap1->Fill(float(zindex), float(ladder));
0970           hcluDetMap1->Fill(col, row);
0971           hcharge1->Fill(charge);
0972           //hcols1->Fill(col);
0973           //hrows1->Fill(row);
0974 
0975           hclusMap1->Fill(gZ, phi);
0976           hmult1->Fill(zindex, float(size));
0977 
0978           if (pt > CLU_SIZE_PT_CUT) {
0979             hsize1->Fill(float(size));
0980             hsizex1->Fill(float(sizeX));
0981             hsizey1->Fill(float(sizeY));
0982 
0983             hclumult1->Fill(eta, float(size));
0984             hclumultx1->Fill(eta, float(sizeX));
0985             hclumulty1->Fill(eta, float(sizeY));
0986             hcluchar1->Fill(eta, float(charge));
0987 
0988             //edm::LogPrint("TestWithTracks")<<ladder<<" "<<ladderOn<<endl;
0989 
0990             hclumultld1->Fill(float(ladderOn), size);
0991             hclumultxld1->Fill(float(ladderOn), sizeX);
0992             hclumultyld1->Fill(float(ladderOn), sizeY);
0993             hclucharld1->Fill(float(ladderOn), charge);
0994           }
0995 
0996 #ifdef VDM_STUDIES
0997           hcharCluls->Fill(lumiBlock, charge);
0998           hsizeCluls->Fill(lumiBlock, size);
0999           hsizeXCluls->Fill(lumiBlock, sizeX);
1000           hcharCluls1->Fill(lumiBlock, charge);
1001           hsizeCluls1->Fill(lumiBlock, size);
1002           hsizeXCluls1->Fill(lumiBlock, sizeX);
1003 #endif
1004 
1005           numOfClusPerTrk1++;
1006           numOfClustersPerLay1++;
1007           numOfPixelsPerLay1 += size;
1008 
1009         } else if (layer == 2) {
1010           hDetMap2->Fill(float(zindex), float(ladder));
1011           hcluDetMap2->Fill(col, row);
1012           hcharge2->Fill(charge);
1013           //hcols2->Fill(col);
1014           //hrows2->Fill(row);
1015 
1016           hclusMap2->Fill(gZ, phi);
1017           hmult2->Fill(zindex, float(size));
1018 
1019           if (pt > CLU_SIZE_PT_CUT) {
1020             hsize2->Fill(float(size));
1021             hsizex2->Fill(float(sizeX));
1022             hsizey2->Fill(float(sizeY));
1023 
1024             hclumult2->Fill(eta, float(size));
1025             hclumultx2->Fill(eta, float(sizeX));
1026             hclumulty2->Fill(eta, float(sizeY));
1027             hcluchar2->Fill(eta, float(charge));
1028 
1029             hclumultld2->Fill(float(ladderOn), size);
1030             hclumultxld2->Fill(float(ladderOn), sizeX);
1031             hclumultyld2->Fill(float(ladderOn), sizeY);
1032             hclucharld2->Fill(float(ladderOn), charge);
1033           }
1034 
1035 #ifdef VDM_STUDIES
1036           hcharCluls->Fill(lumiBlock, charge);
1037           hsizeCluls->Fill(lumiBlock, size);
1038           hsizeXCluls->Fill(lumiBlock, sizeX);
1039           hcharCluls2->Fill(lumiBlock, charge);
1040           hsizeCluls2->Fill(lumiBlock, size);
1041           hsizeXCluls2->Fill(lumiBlock, sizeX);
1042 #endif
1043 
1044           numOfClusPerTrk2++;
1045           numOfClustersPerLay2++;
1046           numOfPixelsPerLay2 += size;
1047 
1048         } else if (layer == 3) {
1049           hDetMap3->Fill(float(zindex), float(ladder));
1050           hcluDetMap3->Fill(col, row);
1051           hcharge3->Fill(charge);
1052           //hcols3->Fill(col);
1053           //hrows3->Fill(row);
1054 
1055           hclusMap3->Fill(gZ, phi);
1056           hmult3->Fill(zindex, float(size));
1057 
1058           if (pt > CLU_SIZE_PT_CUT) {
1059             hsize3->Fill(float(size));
1060             hsizex3->Fill(float(sizeX));
1061             hsizey3->Fill(float(sizeY));
1062             hclumult3->Fill(eta, float(size));
1063             hclumultx3->Fill(eta, float(sizeX));
1064             hclumulty3->Fill(eta, float(sizeY));
1065             hcluchar3->Fill(eta, float(charge));
1066 
1067             hclumultld3->Fill(float(ladderOn), size);
1068             hclumultxld3->Fill(float(ladderOn), sizeX);
1069             hclumultyld3->Fill(float(ladderOn), sizeY);
1070             hclucharld3->Fill(float(ladderOn), charge);
1071           }
1072 
1073 #ifdef VDM_STUDIES
1074           hcharCluls->Fill(lumiBlock, charge);
1075           hsizeCluls->Fill(lumiBlock, size);
1076           hsizeXCluls->Fill(lumiBlock, sizeX);
1077           hcharCluls3->Fill(lumiBlock, charge);
1078           hsizeCluls3->Fill(lumiBlock, size);
1079           hsizeXCluls3->Fill(lumiBlock, sizeX);
1080 #endif
1081 
1082           numOfClusPerTrk3++;
1083           numOfClustersPerLay3++;
1084           numOfPixelsPerLay3 += size;
1085 
1086         } else if (disk == 1) {
1087           numOfClusPerTrk4++;
1088           numOfClustersPerLay4++;
1089           numOfPixelsPerLay4 += size;
1090 
1091           hcharge4->Fill(charge);
1092           if (pt > CLU_SIZE_PT_CUT) {
1093             hsize4->Fill(float(size));
1094             hsizex4->Fill(float(sizeX));
1095             hsizey4->Fill(float(sizeY));
1096           }
1097 
1098           if (side == 1)
1099             numOfClustersPerDisk2++;  // -z
1100           else if (side == 2)
1101             numOfClustersPerDisk3++;  // +z
1102 
1103         } else if (disk == 2) {
1104           numOfClusPerTrk5++;
1105           numOfClustersPerLay5++;
1106           numOfPixelsPerLay5 += size;
1107 
1108           hcharge5->Fill(charge);
1109           if (pt > CLU_SIZE_PT_CUT) {
1110             hsize5->Fill(float(size));
1111             hsizex5->Fill(float(sizeX));
1112             hsizey5->Fill(float(sizeY));
1113           }
1114 
1115           if (side == 1)
1116             numOfClustersPerDisk1++;  // -z
1117           else if (side == 2)
1118             numOfClustersPerDisk4++;  // +z
1119 
1120         } else {
1121           edm::LogPrint("TestWithTracks") << " which layer is this? " << layer << " " << disk;
1122         }  // if layer
1123 
1124       }  // if valid
1125 
1126     }  // clusters
1127 
1128     if (pixelHits > 0)
1129       countPixTracks++;
1130 
1131     if (PRINT)
1132       edm::LogPrint("TestWithTracks") << " Clusters for track " << trackNumber << " num of clusters "
1133                                       << numberOfClusters << " num of pixels " << pixelHits;
1134 
1135 #ifdef HISTOS
1136     // per track histos
1137     if (numberOfClusters > 0) {
1138       hclusPerTrk1->Fill(float(numOfClusPerTrk1));
1139       if (PRINT)
1140         edm::LogPrint("TestWithTracks") << "Lay1: number of clusters per track = " << numOfClusPerTrk1;
1141       hclusPerTrk2->Fill(float(numOfClusPerTrk2));
1142       if (PRINT)
1143         edm::LogPrint("TestWithTracks") << "Lay2: number of clusters per track = " << numOfClusPerTrk1;
1144       hclusPerTrk3->Fill(float(numOfClusPerTrk3));
1145       if (PRINT)
1146         edm::LogPrint("TestWithTracks") << "Lay3: number of clusters per track = " << numOfClusPerTrk1;
1147       hclusPerTrk4->Fill(float(numOfClusPerTrk4));  // fpix  disk1
1148       hclusPerTrk5->Fill(float(numOfClusPerTrk5));  // fpix disk2
1149 
1150       float clusPerTrkB = numOfClusPerTrk1 + numOfClusPerTrk2 + numOfClusPerTrk3;
1151       float clusPerTrkF = numOfClusPerTrk4 + numOfClusPerTrk5;
1152       float clusPerTrk = clusPerTrkB + clusPerTrkF;
1153 
1154       hclusPerTrkB->Fill(clusPerTrkB);
1155       hclusPerTrkF->Fill(clusPerTrkF);
1156       hclusPerTrk->Fill(clusPerTrk);
1157 
1158       hclusPerTrkVsEta->Fill(eta, clusPerTrk);
1159       hclusPerTrkVsEtaB->Fill(eta, clusPerTrkB);
1160       hclusPerTrkVsEtaF->Fill(eta, clusPerTrkF);
1161       hclusPerTrkVsPt->Fill(pt, clusPerTrk);
1162       hclusPerTrkVsls->Fill(lumiBlock, clusPerTrk);
1163     }
1164 #endif  // HISTOS
1165 
1166   }  // tracks
1167 
1168 #ifdef HISTOS
1169   // total layer histos
1170   if (numberOfClusters > 0) {
1171     hclusPerLay1->Fill(float(numOfClustersPerLay1));
1172     hclusPerLay2->Fill(float(numOfClustersPerLay2));
1173     hclusPerLay3->Fill(float(numOfClustersPerLay3));
1174 
1175     hclusPerDisk1->Fill(float(numOfClustersPerDisk1));
1176     hclusPerDisk2->Fill(float(numOfClustersPerDisk2));
1177     hclusPerDisk3->Fill(float(numOfClustersPerDisk3));
1178     hclusPerDisk4->Fill(float(numOfClustersPerDisk4));
1179 
1180     //hdetsPerLay1->Fill(float(numberOfDetUnits1));
1181     //hdetsPerLay2->Fill(float(numberOfDetUnits2));
1182     //hdetsPerLay3->Fill(float(numberOfDetUnits3));
1183     //int tmp = numberOfDetUnits1+numberOfDetUnits2+numberOfDetUnits3;
1184     //hpixPerLay1->Fill(float(numOfPixelsPerLay1));
1185     //hpixPerLay2->Fill(float(numOfPixelsPerLay2));
1186     //hpixPerLay3->Fill(float(numOfPixelsPerLay3));
1187     //htest7->Fill(float(tmp));
1188     hclusBpix->Fill(float(numberOfClusters));
1189     hpixBpix->Fill(float(numberOfPixels));
1190   }
1191   htracksGood->Fill(float(countNiceTracks));
1192   htracksGoodInPix->Fill(float(countPixTracks));
1193   htracks->Fill(float(trackNumber));
1194 
1195   hbx->Fill(float(bx));
1196   hlumi->Fill(float(lumiBlock));
1197 
1198   htracksls->Fill(float(lumiBlock), float(countPixTracks));
1199   hpvsls->Fill(float(lumiBlock), float(pvsTrue));
1200   if (instlumi > 0.) {
1201     float tmp = float(countPixTracks) / instlumi;
1202     htrackslsn->Fill(float(lumiBlock), tmp);
1203     tmp = float(pvsTrue) / instlumi;
1204     hpvslsn->Fill(float(lumiBlock), tmp);
1205   }
1206 
1207 #ifdef VDM_STUDIES
1208 
1209   hclusls->Fill(float(lumiBlock), float(numberOfClusters));  // clusters fpix+bpix
1210   //hpixls->Fill(float(lumiBlock),float(numberOfPixels)); // pixels fpix+bpix
1211 
1212   hclubx->Fill(float(bx), float(numberOfClusters));  // clusters fpix+bpix
1213   //hpixbx->Fill(float(bx),float(numberOfPixels)); // pixels fpix+bpix
1214   hpvbx->Fill(float(bx), float(pvsTrue));            // pvs
1215   htrackbx->Fill(float(bx), float(countPixTracks));  // tracks
1216 
1217   hclusls1->Fill(float(lumiBlock), float(numOfClustersPerLay1));  // clusters bpix1
1218   //hpixls1->Fill( float(lumiBlock),float(numOfPixPerLay1)); // pixels bpix1
1219   hclusls2->Fill(float(lumiBlock), float(numOfClustersPerLay2));  // clusters bpix2
1220   //hpixls2->Fill( float(lumiBlock),float(numOfPixPerLay2)); // pixels bpix2
1221   hclusls3->Fill(float(lumiBlock), float(numOfClustersPerLay3));  // clusters bpix3
1222   //hpixls3->Fill( float(lumiBlock),float(numOfPixPerLay3)); // pixels bpix3
1223 #endif
1224 
1225 #endif  // HISTOS
1226 
1227   //
1228   countTracks += float(trackNumber);
1229   countGoodTracks += float(countNiceTracks);
1230   countTracksInPix += float(countPixTracks);
1231   countPVs += float(pvsTrue);
1232   countEvents++;
1233   if (lumiBlock != lumiBlockOld) {
1234     countLumi += intlumi;
1235     lumiBlockOld = lumiBlock;
1236   }
1237 
1238   if (PRINT)
1239     edm::LogPrint("TestWithTracks") << " event with tracks = " << trackNumber << " " << countNiceTracks;
1240 
1241   return;
1242 
1243 #ifdef USE_TRAJ
1244   // Not used
1245 
1246   //----------------------------------------------------------------------------
1247   // Use Trajectories
1248 
1249   edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
1250   e.getByToken(trackAssocToken_, trajTrackCollectionHandle);
1251 
1252   TrajectoryStateCombiner tsoscomb;
1253 
1254   int NbrTracks = trajTrackCollectionHandle->size();
1255   std::edm::LogPrint("TestWithTracks") << " track measurements " << trajTrackCollectionHandle->size() << std::endl;
1256 
1257   int trackNumber = 0;
1258   int numberOfClusters = 0;
1259 
1260   for (TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(),
1261                                                       itEnd = trajTrackCollectionHandle->end();
1262        it != itEnd;
1263        ++it) {
1264     int pixelHits = 0;
1265     int stripHits = 0;
1266     const Track &track = *it->val;
1267     const Trajectory &traj = *it->key;
1268 
1269     std::vector<TrajectoryMeasurement> checkColl = traj.measurements();
1270     for (std::vector<TrajectoryMeasurement>::const_iterator checkTraj = checkColl.begin(); checkTraj != checkColl.end();
1271          ++checkTraj) {
1272       if (!checkTraj->updatedState().isValid())
1273         continue;
1274       TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit();
1275       if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
1276         continue;
1277       uint testSubDetID = (testhit->geographicalId().subdetId());
1278       if (testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap)
1279         pixelHits++;
1280       else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB ||
1281                testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC)
1282         stripHits++;
1283     }
1284 
1285     if (pixelHits == 0)
1286       continue;
1287 
1288     trackNumber++;
1289     std::edm::LogPrint("TestWithTracks") << " track " << trackNumber << " has pixelhits " << pixelHits << std::endl;
1290     pixelHits = 0;
1291 
1292     //std::vector<TrajectoryMeasurement> tmColl = traj.measurements();
1293     for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = checkColl.begin(); itTraj != checkColl.end();
1294          ++itTraj) {
1295       if (!itTraj->updatedState().isValid())
1296         continue;
1297 
1298       TrajectoryStateOnSurface tsos = tsoscomb(itTraj->forwardPredictedState(), itTraj->backwardPredictedState());
1299       TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
1300       if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker)
1301         continue;
1302 
1303       const DetId &hit_detId = hit->geographicalId();
1304       uint IntSubDetID = (hit_detId.subdetId());
1305 
1306       if (IntSubDetID == 0)
1307         continue;  // Select ??
1308       if (IntSubDetID != PixelSubdetector::PixelBarrel)
1309         continue;  // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) {
1310 
1311       //       const GeomDetUnit* detUnit = hit->detUnit();
1312       //       if(detUnit) {
1313       //    const Surface& surface = hit->detUnit()->surface();
1314       //    const TrackerGeometry& theTracker(*tkGeom_);
1315       //    const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*> (theTracker.idToDet(hit_detId) );
1316       //    const RectangularPixelTopology * topol = dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
1317       //       }
1318 
1319       // get the enclosed persistent hit
1320       const TrackingRecHit *persistentHit = hit->hit();
1321       // check if it's not null, and if it's a valid pixel hit
1322       if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
1323         // tell the C++ compiler that the hit is a pixel hit
1324         const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
1325         // get the edm::Ref to the cluster
1326         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
1327         //  check if the ref is not null
1328         if (clust.isNonnull()) {
1329           numberOfClusters++;
1330           pixelHits++;
1331           float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
1332           int size = clust->size();
1333           int size_x = clust->sizeX();
1334           int size_y = clust->sizeY();
1335           float row = clust->x();
1336           float col = clust->y();
1337 
1338           //LocalPoint lp = topol->localPosition(MeasurementPoint(clust_.row,clust_.col));
1339           //float x = lp.x();
1340           //float y = lp.y();
1341 
1342           int maxPixelCol = clust->maxPixelCol();
1343           int maxPixelRow = clust->maxPixelRow();
1344           int minPixelCol = clust->minPixelCol();
1345           int minPixelRow = clust->minPixelRow();
1346 
1347           //int geoId = PixGeom->geographicalId().rawId();
1348 
1349           // Replace with the topology methods
1350           // edge method moved to topologi class
1351           //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
1352           //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
1353 
1354           // calculate alpha and beta from cluster position
1355           //LocalTrajectoryParameters ltp = tsos.localParameters();
1356           //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
1357 
1358           //float locx = localDir.x();
1359           //float locy = localDir.y();
1360           //float locz = localDir.z();
1361           //float loctheta = localDir.theta(); // currently unused
1362 
1363           //float alpha = atan2( locz, locx );
1364           //float beta = atan2( locz, locy );
1365 
1366           //clust_.normalized_charge = clust_.charge*sqrt(1.0/(1.0/pow(tan(clust_.clust_alpha),2)+1.0/pow(tan(clust_.clust_beta),2)+1.0));
1367         }  // valid cluster
1368       }    // valid peristant hit
1369 
1370     }  // loop over trajectory meas.
1371 
1372     if (PRINT)
1373       edm::LogPrint("TestWithTracks") << " Cluster for track " << trackNumber << " cluaters " << numberOfClusters << " "
1374                                       << pixelHits;
1375 
1376   }  // loop over tracks
1377 
1378 #endif  // USE_TRAJ
1379 
1380   edm::LogPrint("TestWithTracks") << " event with tracks = " << trackNumber << " " << countGoodTracks;
1381 
1382 }  // end
1383 
1384 //define this as a plug-in
1385 DEFINE_FWK_MODULE(TestWithTracks);
[ 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!