Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:21 

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 numOfPixelsPerLay4 = 0;
0504 
0505   int numOfClusPerTrk5 = 0;
0506   //int numOfPixelsPerLay5 = 0;
0507 
0508   int numOfClustersPerDisk1 = 0;
0509   int numOfClustersPerDisk2 = 0;
0510   int numOfClustersPerDisk3 = 0;
0511   int numOfClustersPerDisk4 = 0;
0512 
0513   RunNumber_t const run = e.id().run();
0514   EventNumber_t const event = e.id().event();
0515   LuminosityBlockNumber_t const lumiBlock = e.luminosityBlock();
0516 
0517   int bx = e.bunchCrossing();
0518   //int orbit     = e.orbitNumber(); // unused
0519 
0520   if (PRINT)
0521     edm::LogPrint("TestWithTracks") << "Run " << run << " Event " << event << " LS " << lumiBlock;
0522 
0523   hbx0->Fill(float(bx));
0524   hlumi0->Fill(float(lumiBlock));
0525 
0526   edm::LuminosityBlock const &iLumi = e.getLuminosityBlock();
0527   edm::Handle<LumiSummary> lumi;
0528   iLumi.getByToken(lumiToken_, lumi);
0529   edm::Handle<edm::ConditionsInLumiBlock> cond;
0530   float intlumi = 0, instlumi = 0;
0531   int beamint1 = 0, beamint2 = 0;
0532   iLumi.getByToken(condToken_, cond);
0533   // This will only work when running on RECO until (if) they fix it in the FW
0534   // When running on RAW and reconstructing, the LumiSummary will not appear
0535   // in the event before reaching endLuminosityBlock(). Therefore, it is not
0536   // possible to get this info in the event
0537   if (lumi.isValid()) {
0538     intlumi = (lumi->intgRecLumi()) / 1000.;  // 10^30 -> 10^33/cm2/sec  ->  1/nb/sec
0539     instlumi = (lumi->avgInsDelLumi()) / 1000.;
0540     beamint1 = (cond->totalIntensityBeam1) / 1000;
0541     beamint2 = (cond->totalIntensityBeam2) / 1000;
0542   } else {
0543     //std::edm::LogPrint("TestWithTracks") << "** ERROR: Event does not get lumi info\n";
0544   }
0545 
0546   hinstl->Fill(float(lumiBlock), float(instlumi));
0547   hintgl->Fill(float(lumiBlock), float(intlumi));
0548   hbeam1->Fill(float(lumiBlock), float(beamint1));
0549   hbeam2->Fill(float(lumiBlock), float(beamint2));
0550 
0551 #ifdef L1
0552   // Get L1
0553   Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
0554   e.getByToken(l1gtrrToken_, L1GTRR);
0555 
0556   if (L1GTRR.isValid()) {
0557     //bool l1a = L1GTRR->decision();  // global decission?
0558     //edm::LogPrint("TestWithTracks")<<" L1 status :"<<l1a<<" "<<hex;
0559     for (unsigned int i = 0; i < L1GTRR->decisionWord().size(); ++i) {
0560       int l1flag = L1GTRR->decisionWord()[i];
0561       int t1flag = L1GTRR->technicalTriggerWord()[i];
0562 
0563       if (l1flag > 0)
0564         hl1a->Fill(float(i));
0565       if (t1flag > 0 && i < 64)
0566         hl1t->Fill(float(i));
0567     }  // for loop
0568   }    // if l1a
0569 #endif
0570 
0571 #ifdef HLT
0572 
0573   bool hlt[256];
0574   for (int i = 0; i < 256; ++i)
0575     hlt[i] = false;
0576 
0577   edm::TriggerNames TrigNames;
0578   edm::Handle<edm::TriggerResults> HLTResults;
0579 
0580   // Extract the HLT results
0581   e.getByToken(hltToken_, HLTResults);
0582   if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) {
0583     //TrigNames.init(*HLTResults);
0584     const edm::TriggerNames &TrigNames = e.triggerNames(*HLTResults);
0585 
0586     //edm::LogPrint("TestWithTracks")<<TrigNames.triggerNames().size()<<endl;
0587 
0588     for (unsigned int i = 0; i < TrigNames.triggerNames().size(); i++) {  // loop over trigger
0589       //if(countAllEvents==1) edm::LogPrint("TestWithTracks")<<i<<" "<<TrigNames.triggerName(i)<<endl;
0590 
0591       if ((HLTResults->wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0592           (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0593           (HLTResults->error(TrigNames.triggerIndex(TrigNames.triggerName(i))) == false)) {
0594         hlt[i] = true;
0595         hlt1->Fill(float(i));
0596 
0597       }  // if hlt
0598 
0599     }  // loop
0600   }    // if valid
0601 #endif
0602 
0603   // Get event setup
0604   edm::ESHandle<TrackerGeometry> geom = es.getHandle(trackerGeomToken_);
0605   const TrackerGeometry &theTracker(*geom);
0606 
0607   // -- Primary vertices
0608   // ----------------------------------------------------------------------
0609   edm::Handle<reco::VertexCollection> vertices;
0610   e.getByToken(vtxToken_, vertices);
0611 
0612   if (PRINT)
0613     edm::LogPrint("TestWithTracks") << " PV list " << vertices->size();
0614   int pvNotFake = 0, pvsTrue = 0;
0615   vector<float> pvzVector;
0616   for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) {
0617     if ((iv->isFake()) == 1)
0618       continue;
0619     pvNotFake++;
0620     float pvx = iv->x();
0621     float pvy = iv->y();
0622     float pvz = iv->z();
0623     int numTracksPerPV = iv->tracksSize();
0624     //int numTracksPerPV = iv->nTracks();
0625 
0626     //float xe = iv->xError();
0627     //float ye = iv->yError();
0628     //float ze = iv->zError();
0629     //int chi2 = iv->chi2();
0630     //int dof = iv->ndof();
0631 
0632     if (PRINT)
0633       edm::LogPrint("TestWithTracks") << " PV " << pvNotFake << " pos = " << pvx << "/" << pvy << "/" << pvz
0634                                       << ", Num of tracks " << numTracksPerPV;
0635 
0636     hpvz->Fill(pvz);
0637     if (pvz > -22. && pvz < 22.) {
0638       float pvr = sqrt(pvx * pvx + pvy * pvy);
0639       hpvxy->Fill(pvx, pvy);
0640       hpvr->Fill(pvr);
0641       if (pvr < 0.3) {
0642         pvsTrue++;
0643         pvzVector.push_back(pvz);
0644         //if(PRINT) edm::LogPrint("TestWithTracks")<<"PV "<<pvsTrue<<" "<<pvz<<endl;
0645       }  //pvr
0646     }    // pvz
0647 
0648     //if(pvsTrue<1) continue; // skip events with no PV
0649 
0650   }  // loop pvs
0651   hNumPv->Fill(float(pvNotFake));
0652   hNumPvClean->Fill(float(pvsTrue));
0653 
0654   if (PRINT)
0655     edm::LogPrint("TestWithTracks") << " Not fake PVs = " << pvNotFake << " good position " << pvsTrue;
0656 
0657   // -- Tracks
0658   // ----------------------------------------------------------------------
0659   Handle<reco::TrackCollection> recTracks;
0660   e.getByToken(srcToken_, recTracks);
0661 
0662   if (PRINT)
0663     edm::LogPrint("TestWithTracks") << " Tracks " << recTracks->size();
0664   for (reco::TrackCollection::const_iterator t = recTracks->begin(); t != recTracks->end(); ++t) {
0665     trackNumber++;
0666     numOfClusPerTrk1 = 0;  // this is confusing, it is used as clus per track
0667     numOfClusPerTrk2 = 0;
0668     numOfClusPerTrk3 = 0;
0669     numOfClusPerTrk4 = 0;
0670     numOfClusPerTrk5 = 0;
0671     int pixelHits = 0;
0672 
0673     int size = t->recHitsSize();
0674     float pt = t->pt();
0675     float eta = t->eta();
0676     float phi = t->phi();
0677     //float trackCharge = t->charge(); // unused
0678     float d0 = t->d0();
0679     float dz = t->dz();
0680     //float tkvx = t->vx();  // unused
0681     //float tkvy = t->vy();
0682     //float tkvz = t->vz();
0683 
0684     if (PRINT)
0685       edm::LogPrint("TestWithTracks") << "Track " << trackNumber << " Pt " << pt << " Eta " << eta << " d0/dz " << d0
0686                                       << " " << dz << " Hits " << size;
0687 
0688     hEta->Fill(eta);
0689     hDz->Fill(dz);
0690     if (abs(eta) > 2.8 || abs(dz) > 25.)
0691       continue;  //  skip
0692 
0693     hD0->Fill(d0);
0694     if (d0 > 1.0)
0695       continue;  // skip
0696 
0697     bool goodTrack = false;
0698     for (vector<float>::iterator m = pvzVector.begin(); m != pvzVector.end(); ++m) {
0699       float z = *m;
0700       float tmp = abs(z - dz);
0701       hzdiff->Fill(tmp);
0702       if (tmp < 1.)
0703         goodTrack = true;
0704     }
0705 
0706     if (isData && !goodTrack)
0707       continue;
0708     countNiceTracks++;
0709     hPt->Fill(pt);
0710 
0711     // Loop over rechits
0712     for (trackingRecHit_iterator recHit = t->recHitsBegin(); recHit != t->recHitsEnd(); ++recHit) {
0713       if (!((*recHit)->isValid()))
0714         continue;
0715 
0716       if ((*recHit)->geographicalId().det() != DetId::Tracker)
0717         continue;
0718 
0719       const DetId &hit_detId = (*recHit)->geographicalId();
0720       uint IntSubDetID = (hit_detId.subdetId());
0721 
0722       // Select pixel rechits
0723       if (IntSubDetID == 0)
0724         continue;  // Select ??
0725 
0726       int layer = 0, ladder = 0, zindex = 0, ladderOn = 0, module = 0, shell = 0;
0727       unsigned int disk = 0;     //1,2,3
0728       unsigned int blade = 0;    //1-24
0729       unsigned int zindexF = 0;  //
0730       unsigned int side = 0;     //size=1 for -z, 2 for +z
0731       unsigned int panel = 0;    //panel=1
0732 
0733       if (IntSubDetID == PixelSubdetector::PixelBarrel) {  // bpix
0734 
0735         // Pixel detector
0736         PXBDetId pdetId = PXBDetId(hit_detId);
0737         //unsigned int detTypeP=pdetId.det();
0738         //unsigned int subidP=pdetId.subdetId();
0739         // Barell layer = 1,2,3
0740         layer = pdetId.layer();
0741         // Barrel ladder id 1-20,32,44.
0742         ladder = pdetId.ladder();
0743         // Barrel Z-index=1,8
0744         zindex = pdetId.module();
0745         if (zindex < 5)
0746           side = 1;
0747         else
0748           side = 2;
0749 
0750         // Convert to online
0751         PixelBarrelName pbn(pdetId);
0752         // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0753         PixelBarrelName::Shell sh = pbn.shell();  //enum
0754         //sector = pbn.sectorName();
0755         ladderOn = pbn.ladderName();
0756         //layerOn  = pbn.layerName();
0757         module = pbn.moduleName();  // 1 to 4
0758         //half  = pbn.isHalfModule();
0759         shell = int(sh);
0760         // change the module sign for z<0
0761         if (shell == 1 || shell == 2)
0762           module = -module;  // make -1 to -4 for -z
0763                              // change ladeer sign for Outer )x<0)
0764         if (shell == 1 || shell == 3)
0765           ladderOn = -ladderOn;
0766 
0767         if (PRINT)
0768           edm::LogPrint("TestWithTracks") << "barrel layer/ladder/module: " << layer << "/" << ladder << "/" << zindex;
0769 
0770       } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {  // fpix
0771 
0772         PXFDetId pdetId = PXFDetId(hit_detId);
0773         disk = pdetId.disk();       //1,2,3
0774         blade = pdetId.blade();     //1-24
0775         zindexF = pdetId.module();  //
0776         side = pdetId.side();       //size=1 for -z, 2 for +z
0777         panel = pdetId.panel();     //panel=1
0778 
0779         if (PRINT)
0780           edm::LogPrint("TestWithTracks") << " forward det, disk " << disk << ", blade " << blade << ", module "
0781                                           << zindexF << ", side " << side << ", panel " << panel;
0782 
0783       } else {     // nothings
0784         continue;  // skip
0785       }
0786 
0787       // Get the geom-detector
0788       const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
0789       //double detZ = theGeomDet->surface().position().z();  // unused
0790       //double detR = theGeomDet->surface().position().perp(); // unused
0791       const PixelTopology *topol = &(theGeomDet->specificTopology());  // pixel topology
0792 
0793       //std::vector<SiStripRecHit2D*> output = getRecHitComponents((*recHit).get());
0794       //std::vector<SiPixelRecHit*> TrkComparison::getRecHitComponents(const TrackingRecHit* rechit){
0795 
0796       const SiPixelRecHit *hit = dynamic_cast<const SiPixelRecHit *>((*recHit));
0797       //edm::Ref<edmNew::DetSetVector<SiStripCluster> ,SiStripCluster> cluster = hit->cluster();
0798       // get the edm::Ref to the cluster
0799 
0800       if (hit) {
0801         if (pt > 1.) {  // eliminate low pt tracks
0802           // RecHit (recthits are transient, so not available without track refit)
0803           double xloc = hit->localPosition().x();  // 1st meas coord
0804           double yloc = hit->localPosition().y();  // 2nd meas coord or zero
0805           //double zloc = hit->localPosition().z();// up, always zero
0806           LocalError lerr = hit->localPositionError();
0807           float lerr_x = sqrt(lerr.xx()) * 1E4;
0808           float lerr_y = sqrt(lerr.yy()) * 1E4;
0809 
0810           if (layer == 1) {
0811             recHitXError1->Fill(lerr_x);
0812             recHitYError1->Fill(lerr_y);
0813             hErrorXB->Fill(float(ladder + (110 * (side - 1))), lerr_x);
0814             hErrorYB->Fill(float(ladder + (110 * (side - 1))), lerr_y);
0815           } else if (layer == 2) {
0816             recHitXError2->Fill(lerr_x);
0817             recHitYError2->Fill(lerr_y);
0818             hErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), lerr_x);
0819             hErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), lerr_y);
0820 
0821           } else if (layer == 3) {
0822             recHitXError3->Fill(lerr_x);
0823             recHitYError3->Fill(lerr_y);
0824             hErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), lerr_x);
0825             hErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), lerr_y);
0826           } else if ((disk == 2) && (side == 1)) {
0827             recHitXError4->Fill(lerr_x);
0828             recHitYError4->Fill(lerr_y);
0829             hErrorXF->Fill(float(blade), lerr_x);
0830             hErrorYF->Fill(float(blade), lerr_y);
0831 
0832           } else if ((disk == 1) && (side == 1)) {
0833             recHitXError5->Fill(lerr_x);
0834             recHitYError5->Fill(lerr_y);
0835             hErrorXF->Fill(float(blade + 25), lerr_x);
0836             hErrorYF->Fill(float(blade + 25), lerr_y);
0837 
0838           } else if ((disk == 1) && (side == 2)) {
0839             recHitXError6->Fill(lerr_x);
0840             recHitYError6->Fill(lerr_y);
0841             hErrorXF->Fill(float(blade + 50), lerr_x);
0842             hErrorYF->Fill(float(blade + 50), lerr_y);
0843           } else if ((disk == 2) && (side == 2)) {
0844             recHitXError7->Fill(lerr_x);
0845             recHitYError7->Fill(lerr_y);
0846             hErrorXF->Fill(float(blade + 75), lerr_x);
0847             hErrorYF->Fill(float(blade + 75), lerr_y);
0848           }
0849 
0850           LocalError lape = theGeomDet->localAlignmentError();
0851           if (lape.valid()) {
0852             float tmp11 = 0.;
0853             if (lape.xx() > 0.)
0854               tmp11 = sqrt(lape.xx()) * 1E4;
0855             //float tmp12= sqrt(lape.xy())*1E4;
0856             float tmp13 = 0.;
0857             if (lape.yy() > 0.)
0858               tmp13 = sqrt(lape.yy()) * 1E4;
0859             //bool tmp14=tmp2<tmp1;
0860             if (layer == 1) {
0861               recHitXAlignError1->Fill(tmp11);
0862               recHitYAlignError1->Fill(tmp13);
0863               hAErrorXB->Fill(float(ladder + (110 * (side - 1))), tmp11);
0864               hAErrorYB->Fill(float(ladder + (110 * (side - 1))), tmp13);
0865             } else if (layer == 2) {
0866               recHitXAlignError2->Fill(tmp11);
0867               recHitYAlignError2->Fill(tmp13);
0868               hAErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), tmp11);
0869               hAErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), tmp13);
0870             } else if (layer == 3) {
0871               recHitXAlignError3->Fill(tmp11);
0872               recHitYAlignError3->Fill(tmp13);
0873               hAErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), tmp11);
0874               hAErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), tmp13);
0875 
0876             } else if ((disk == 2) && (side == 1)) {
0877               recHitXAlignError4->Fill(tmp11);
0878               recHitYAlignError4->Fill(tmp13);
0879               hAErrorXF->Fill(float(blade), tmp11);
0880               hAErrorYF->Fill(float(blade), tmp13);
0881 
0882             } else if ((disk == 1) && (side == 1)) {
0883               recHitXAlignError5->Fill(tmp11);
0884               recHitYAlignError5->Fill(tmp13);
0885               hAErrorXF->Fill(float(blade + 25), tmp11);
0886               hAErrorYF->Fill(float(blade + 25), tmp13);
0887             } else if ((disk == 1) && (side == 2)) {
0888               recHitXAlignError6->Fill(tmp11);
0889               recHitYAlignError6->Fill(tmp13);
0890               hAErrorXF->Fill(float(blade + 50), tmp11);
0891               hAErrorYF->Fill(float(blade + 50), tmp13);
0892             } else if ((disk == 2) && (side == 2)) {
0893               recHitXAlignError7->Fill(tmp11);
0894               recHitYAlignError7->Fill(tmp13);
0895               hAErrorXF->Fill(float(blade + 75), tmp11);
0896               hAErrorYF->Fill(float(blade + 75), tmp13);
0897             }
0898 
0899             //edm::LogPrint("TestWithTracks")<<tTopo->pxbLayer(detId)<<" "<<tTopo->pxbModule(detId)<<" "<<rows<<" "<<tmp14<<" "
0900             if (PRINT)
0901               edm::LogPrint("TestWithTracks") << " align error " << layer << tmp11 << " " << tmp13;
0902           } else {
0903             edm::LogPrint("TestWithTracks") << " lape = 0";
0904           }  // if lape
0905 
0906           if (PRINT)
0907             edm::LogPrint("TestWithTracks") << " rechit loc " << xloc << " " << yloc << " " << lerr_x << " " << lerr_y;
0908         }  // limit pt
0909 
0910         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = hit->cluster();
0911         //  check if the ref is not null
0912         if (!clust.isNonnull())
0913           continue;
0914 
0915         numberOfClusters++;
0916         pixelHits++;
0917         float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
0918         int size = clust->size();
0919         int sizeX = clust->sizeX();
0920         int sizeY = clust->sizeY();
0921         float row = clust->x();
0922         float col = clust->y();
0923         numberOfPixels += size;
0924 
0925         //edm::LogPrint("TestWithTracks")<<" clus loc "<<row<<" "<<col<<endl;
0926 
0927         if (PRINT)
0928           edm::LogPrint("TestWithTracks")
0929               << " cluster " << numberOfClusters << " charge = " << charge << " size = " << size;
0930 
0931         LocalPoint lp = topol->localPosition(MeasurementPoint(clust->x(), clust->y()));
0932         //float x = lp.x();
0933         //float y = lp.y();
0934         //edm::LogPrint("TestWithTracks")<<" clu loc "<<x<<" "<<y<<endl;
0935 
0936         GlobalPoint clustgp = theGeomDet->surface().toGlobal(lp);
0937         double gX = clustgp.x();
0938         double gY = clustgp.y();
0939         double gZ = clustgp.z();
0940 
0941         //edm::LogPrint("TestWithTracks")<<" CLU GLOBAL "<<gX<<" "<<gY<<" "<<gZ<<endl;
0942 
0943         TVector3 v(gX, gY, gZ);
0944         //float phi = v.Phi(); // unused
0945 
0946         //int maxPixelCol = clust->maxPixelCol();
0947         //int maxPixelRow = clust->maxPixelRow();
0948         //int minPixelCol = clust->minPixelCol();
0949         //int minPixelRow = clust->minPixelRow();
0950         //int geoId = PixGeom->geographicalId().rawId();
0951         // Replace with the topology methods
0952         // edge method moved to topologi class
0953         //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
0954         //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
0955 
0956         // calculate alpha and beta from cluster position
0957         //LocalTrajectoryParameters ltp = tsos.localParameters();
0958         //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0959         //float locx = localDir.x();
0960         //float locy = localDir.y();
0961         //float locz = localDir.z();
0962         //float loctheta = localDir.theta(); // currently unused
0963         //float alpha = atan2( locz, locx );
0964         //float beta = atan2( locz, locy );
0965 
0966         if (layer == 1) {
0967           hDetMap1->Fill(float(zindex), float(ladder));
0968           hcluDetMap1->Fill(col, row);
0969           hcharge1->Fill(charge);
0970           //hcols1->Fill(col);
0971           //hrows1->Fill(row);
0972 
0973           hclusMap1->Fill(gZ, phi);
0974           hmult1->Fill(zindex, float(size));
0975 
0976           if (pt > CLU_SIZE_PT_CUT) {
0977             hsize1->Fill(float(size));
0978             hsizex1->Fill(float(sizeX));
0979             hsizey1->Fill(float(sizeY));
0980 
0981             hclumult1->Fill(eta, float(size));
0982             hclumultx1->Fill(eta, float(sizeX));
0983             hclumulty1->Fill(eta, float(sizeY));
0984             hcluchar1->Fill(eta, float(charge));
0985 
0986             //edm::LogPrint("TestWithTracks")<<ladder<<" "<<ladderOn<<endl;
0987 
0988             hclumultld1->Fill(float(ladderOn), size);
0989             hclumultxld1->Fill(float(ladderOn), sizeX);
0990             hclumultyld1->Fill(float(ladderOn), sizeY);
0991             hclucharld1->Fill(float(ladderOn), charge);
0992           }
0993 
0994 #ifdef VDM_STUDIES
0995           hcharCluls->Fill(lumiBlock, charge);
0996           hsizeCluls->Fill(lumiBlock, size);
0997           hsizeXCluls->Fill(lumiBlock, sizeX);
0998           hcharCluls1->Fill(lumiBlock, charge);
0999           hsizeCluls1->Fill(lumiBlock, size);
1000           hsizeXCluls1->Fill(lumiBlock, sizeX);
1001 #endif
1002 
1003           numOfClusPerTrk1++;
1004           numOfClustersPerLay1++;
1005           //numOfPixelsPerLay1 += size;
1006 
1007         } else if (layer == 2) {
1008           hDetMap2->Fill(float(zindex), float(ladder));
1009           hcluDetMap2->Fill(col, row);
1010           hcharge2->Fill(charge);
1011           //hcols2->Fill(col);
1012           //hrows2->Fill(row);
1013 
1014           hclusMap2->Fill(gZ, phi);
1015           hmult2->Fill(zindex, float(size));
1016 
1017           if (pt > CLU_SIZE_PT_CUT) {
1018             hsize2->Fill(float(size));
1019             hsizex2->Fill(float(sizeX));
1020             hsizey2->Fill(float(sizeY));
1021 
1022             hclumult2->Fill(eta, float(size));
1023             hclumultx2->Fill(eta, float(sizeX));
1024             hclumulty2->Fill(eta, float(sizeY));
1025             hcluchar2->Fill(eta, float(charge));
1026 
1027             hclumultld2->Fill(float(ladderOn), size);
1028             hclumultxld2->Fill(float(ladderOn), sizeX);
1029             hclumultyld2->Fill(float(ladderOn), sizeY);
1030             hclucharld2->Fill(float(ladderOn), charge);
1031           }
1032 
1033 #ifdef VDM_STUDIES
1034           hcharCluls->Fill(lumiBlock, charge);
1035           hsizeCluls->Fill(lumiBlock, size);
1036           hsizeXCluls->Fill(lumiBlock, sizeX);
1037           hcharCluls2->Fill(lumiBlock, charge);
1038           hsizeCluls2->Fill(lumiBlock, size);
1039           hsizeXCluls2->Fill(lumiBlock, sizeX);
1040 #endif
1041 
1042           numOfClusPerTrk2++;
1043           numOfClustersPerLay2++;
1044           //numOfPixelsPerLay2 += size;
1045 
1046         } else if (layer == 3) {
1047           hDetMap3->Fill(float(zindex), float(ladder));
1048           hcluDetMap3->Fill(col, row);
1049           hcharge3->Fill(charge);
1050           //hcols3->Fill(col);
1051           //hrows3->Fill(row);
1052 
1053           hclusMap3->Fill(gZ, phi);
1054           hmult3->Fill(zindex, float(size));
1055 
1056           if (pt > CLU_SIZE_PT_CUT) {
1057             hsize3->Fill(float(size));
1058             hsizex3->Fill(float(sizeX));
1059             hsizey3->Fill(float(sizeY));
1060             hclumult3->Fill(eta, float(size));
1061             hclumultx3->Fill(eta, float(sizeX));
1062             hclumulty3->Fill(eta, float(sizeY));
1063             hcluchar3->Fill(eta, float(charge));
1064 
1065             hclumultld3->Fill(float(ladderOn), size);
1066             hclumultxld3->Fill(float(ladderOn), sizeX);
1067             hclumultyld3->Fill(float(ladderOn), sizeY);
1068             hclucharld3->Fill(float(ladderOn), charge);
1069           }
1070 
1071 #ifdef VDM_STUDIES
1072           hcharCluls->Fill(lumiBlock, charge);
1073           hsizeCluls->Fill(lumiBlock, size);
1074           hsizeXCluls->Fill(lumiBlock, sizeX);
1075           hcharCluls3->Fill(lumiBlock, charge);
1076           hsizeCluls3->Fill(lumiBlock, size);
1077           hsizeXCluls3->Fill(lumiBlock, sizeX);
1078 #endif
1079 
1080           numOfClusPerTrk3++;
1081           numOfClustersPerLay3++;
1082           //numOfPixelsPerLay3 += size;
1083 
1084         } else if (disk == 1) {
1085           numOfClusPerTrk4++;
1086           //numOfPixelsPerLay4 += size;
1087 
1088           hcharge4->Fill(charge);
1089           if (pt > CLU_SIZE_PT_CUT) {
1090             hsize4->Fill(float(size));
1091             hsizex4->Fill(float(sizeX));
1092             hsizey4->Fill(float(sizeY));
1093           }
1094 
1095           if (side == 1)
1096             numOfClustersPerDisk2++;  // -z
1097           else if (side == 2)
1098             numOfClustersPerDisk3++;  // +z
1099 
1100         } else if (disk == 2) {
1101           numOfClusPerTrk5++;
1102           //numOfPixelsPerLay5 += size;
1103 
1104           hcharge5->Fill(charge);
1105           if (pt > CLU_SIZE_PT_CUT) {
1106             hsize5->Fill(float(size));
1107             hsizex5->Fill(float(sizeX));
1108             hsizey5->Fill(float(sizeY));
1109           }
1110 
1111           if (side == 1)
1112             numOfClustersPerDisk1++;  // -z
1113           else if (side == 2)
1114             numOfClustersPerDisk4++;  // +z
1115 
1116         } else {
1117           edm::LogPrint("TestWithTracks") << " which layer is this? " << layer << " " << disk;
1118         }  // if layer
1119 
1120       }  // if valid
1121 
1122     }  // clusters
1123 
1124     if (pixelHits > 0)
1125       countPixTracks++;
1126 
1127     if (PRINT)
1128       edm::LogPrint("TestWithTracks") << " Clusters for track " << trackNumber << " num of clusters "
1129                                       << numberOfClusters << " num of pixels " << pixelHits;
1130 
1131 #ifdef HISTOS
1132     // per track histos
1133     if (numberOfClusters > 0) {
1134       hclusPerTrk1->Fill(float(numOfClusPerTrk1));
1135       if (PRINT)
1136         edm::LogPrint("TestWithTracks") << "Lay1: number of clusters per track = " << numOfClusPerTrk1;
1137       hclusPerTrk2->Fill(float(numOfClusPerTrk2));
1138       if (PRINT)
1139         edm::LogPrint("TestWithTracks") << "Lay2: number of clusters per track = " << numOfClusPerTrk1;
1140       hclusPerTrk3->Fill(float(numOfClusPerTrk3));
1141       if (PRINT)
1142         edm::LogPrint("TestWithTracks") << "Lay3: number of clusters per track = " << numOfClusPerTrk1;
1143       hclusPerTrk4->Fill(float(numOfClusPerTrk4));  // fpix  disk1
1144       hclusPerTrk5->Fill(float(numOfClusPerTrk5));  // fpix disk2
1145 
1146       float clusPerTrkB = numOfClusPerTrk1 + numOfClusPerTrk2 + numOfClusPerTrk3;
1147       float clusPerTrkF = numOfClusPerTrk4 + numOfClusPerTrk5;
1148       float clusPerTrk = clusPerTrkB + clusPerTrkF;
1149 
1150       hclusPerTrkB->Fill(clusPerTrkB);
1151       hclusPerTrkF->Fill(clusPerTrkF);
1152       hclusPerTrk->Fill(clusPerTrk);
1153 
1154       hclusPerTrkVsEta->Fill(eta, clusPerTrk);
1155       hclusPerTrkVsEtaB->Fill(eta, clusPerTrkB);
1156       hclusPerTrkVsEtaF->Fill(eta, clusPerTrkF);
1157       hclusPerTrkVsPt->Fill(pt, clusPerTrk);
1158       hclusPerTrkVsls->Fill(lumiBlock, clusPerTrk);
1159     }
1160 #endif  // HISTOS
1161 
1162   }  // tracks
1163 
1164 #ifdef HISTOS
1165   // total layer histos
1166   if (numberOfClusters > 0) {
1167     hclusPerLay1->Fill(float(numOfClustersPerLay1));
1168     hclusPerLay2->Fill(float(numOfClustersPerLay2));
1169     hclusPerLay3->Fill(float(numOfClustersPerLay3));
1170 
1171     hclusPerDisk1->Fill(float(numOfClustersPerDisk1));
1172     hclusPerDisk2->Fill(float(numOfClustersPerDisk2));
1173     hclusPerDisk3->Fill(float(numOfClustersPerDisk3));
1174     hclusPerDisk4->Fill(float(numOfClustersPerDisk4));
1175 
1176     //hdetsPerLay1->Fill(float(numberOfDetUnits1));
1177     //hdetsPerLay2->Fill(float(numberOfDetUnits2));
1178     //hdetsPerLay3->Fill(float(numberOfDetUnits3));
1179     //int tmp = numberOfDetUnits1+numberOfDetUnits2+numberOfDetUnits3;
1180     //hpixPerLay1->Fill(float(numOfPixelsPerLay1));
1181     //hpixPerLay2->Fill(float(numOfPixelsPerLay2));
1182     //hpixPerLay3->Fill(float(numOfPixelsPerLay3));
1183     //htest7->Fill(float(tmp));
1184     hclusBpix->Fill(float(numberOfClusters));
1185     hpixBpix->Fill(float(numberOfPixels));
1186   }
1187   htracksGood->Fill(float(countNiceTracks));
1188   htracksGoodInPix->Fill(float(countPixTracks));
1189   htracks->Fill(float(trackNumber));
1190 
1191   hbx->Fill(float(bx));
1192   hlumi->Fill(float(lumiBlock));
1193 
1194   htracksls->Fill(float(lumiBlock), float(countPixTracks));
1195   hpvsls->Fill(float(lumiBlock), float(pvsTrue));
1196   if (instlumi > 0.) {
1197     float tmp = float(countPixTracks) / instlumi;
1198     htrackslsn->Fill(float(lumiBlock), tmp);
1199     tmp = float(pvsTrue) / instlumi;
1200     hpvslsn->Fill(float(lumiBlock), tmp);
1201   }
1202 
1203 #ifdef VDM_STUDIES
1204 
1205   hclusls->Fill(float(lumiBlock), float(numberOfClusters));  // clusters fpix+bpix
1206   //hpixls->Fill(float(lumiBlock),float(numberOfPixels)); // pixels fpix+bpix
1207 
1208   hclubx->Fill(float(bx), float(numberOfClusters));  // clusters fpix+bpix
1209   //hpixbx->Fill(float(bx),float(numberOfPixels)); // pixels fpix+bpix
1210   hpvbx->Fill(float(bx), float(pvsTrue));            // pvs
1211   htrackbx->Fill(float(bx), float(countPixTracks));  // tracks
1212 
1213   hclusls1->Fill(float(lumiBlock), float(numOfClustersPerLay1));  // clusters bpix1
1214   //hpixls1->Fill( float(lumiBlock),float(numOfPixPerLay1)); // pixels bpix1
1215   hclusls2->Fill(float(lumiBlock), float(numOfClustersPerLay2));  // clusters bpix2
1216   //hpixls2->Fill( float(lumiBlock),float(numOfPixPerLay2)); // pixels bpix2
1217   hclusls3->Fill(float(lumiBlock), float(numOfClustersPerLay3));  // clusters bpix3
1218   //hpixls3->Fill( float(lumiBlock),float(numOfPixPerLay3)); // pixels bpix3
1219 #endif
1220 
1221 #endif  // HISTOS
1222 
1223   //
1224   countTracks += float(trackNumber);
1225   countGoodTracks += float(countNiceTracks);
1226   countTracksInPix += float(countPixTracks);
1227   countPVs += float(pvsTrue);
1228   countEvents++;
1229   if (lumiBlock != lumiBlockOld) {
1230     countLumi += intlumi;
1231     lumiBlockOld = lumiBlock;
1232   }
1233 
1234   if (PRINT)
1235     edm::LogPrint("TestWithTracks") << " event with tracks = " << trackNumber << " " << countNiceTracks;
1236 
1237   return;
1238 
1239 #ifdef USE_TRAJ
1240   // Not used
1241 
1242   //----------------------------------------------------------------------------
1243   // Use Trajectories
1244 
1245   edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
1246   e.getByToken(trackAssocToken_, trajTrackCollectionHandle);
1247 
1248   TrajectoryStateCombiner tsoscomb;
1249 
1250   int NbrTracks = trajTrackCollectionHandle->size();
1251   std::edm::LogPrint("TestWithTracks") << " track measurements " << trajTrackCollectionHandle->size() << std::endl;
1252 
1253   int trackNumber = 0;
1254   int numberOfClusters = 0;
1255 
1256   for (TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(),
1257                                                       itEnd = trajTrackCollectionHandle->end();
1258        it != itEnd;
1259        ++it) {
1260     int pixelHits = 0;
1261     int stripHits = 0;
1262     const Track &track = *it->val;
1263     const Trajectory &traj = *it->key;
1264 
1265     std::vector<TrajectoryMeasurement> checkColl = traj.measurements();
1266     for (std::vector<TrajectoryMeasurement>::const_iterator checkTraj = checkColl.begin(); checkTraj != checkColl.end();
1267          ++checkTraj) {
1268       if (!checkTraj->updatedState().isValid())
1269         continue;
1270       TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit();
1271       if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
1272         continue;
1273       uint testSubDetID = (testhit->geographicalId().subdetId());
1274       if (testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap)
1275         pixelHits++;
1276       else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB ||
1277                testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC)
1278         stripHits++;
1279     }
1280 
1281     if (pixelHits == 0)
1282       continue;
1283 
1284     trackNumber++;
1285     std::edm::LogPrint("TestWithTracks") << " track " << trackNumber << " has pixelhits " << pixelHits << std::endl;
1286     pixelHits = 0;
1287 
1288     //std::vector<TrajectoryMeasurement> tmColl = traj.measurements();
1289     for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = checkColl.begin(); itTraj != checkColl.end();
1290          ++itTraj) {
1291       if (!itTraj->updatedState().isValid())
1292         continue;
1293 
1294       TrajectoryStateOnSurface tsos = tsoscomb(itTraj->forwardPredictedState(), itTraj->backwardPredictedState());
1295       TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
1296       if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker)
1297         continue;
1298 
1299       const DetId &hit_detId = hit->geographicalId();
1300       uint IntSubDetID = (hit_detId.subdetId());
1301 
1302       if (IntSubDetID == 0)
1303         continue;  // Select ??
1304       if (IntSubDetID != PixelSubdetector::PixelBarrel)
1305         continue;  // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) {
1306 
1307       //       const GeomDetUnit* detUnit = hit->detUnit();
1308       //       if(detUnit) {
1309       //    const Surface& surface = hit->detUnit()->surface();
1310       //    const TrackerGeometry& theTracker(*tkGeom_);
1311       //    const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*> (theTracker.idToDet(hit_detId) );
1312       //    const RectangularPixelTopology * topol = dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
1313       //       }
1314 
1315       // get the enclosed persistent hit
1316       const TrackingRecHit *persistentHit = hit->hit();
1317       // check if it's not null, and if it's a valid pixel hit
1318       if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
1319         // tell the C++ compiler that the hit is a pixel hit
1320         const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
1321         // get the edm::Ref to the cluster
1322         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
1323         //  check if the ref is not null
1324         if (clust.isNonnull()) {
1325           numberOfClusters++;
1326           pixelHits++;
1327           float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
1328           int size = clust->size();
1329           int size_x = clust->sizeX();
1330           int size_y = clust->sizeY();
1331           float row = clust->x();
1332           float col = clust->y();
1333 
1334           //LocalPoint lp = topol->localPosition(MeasurementPoint(clust_.row,clust_.col));
1335           //float x = lp.x();
1336           //float y = lp.y();
1337 
1338           int maxPixelCol = clust->maxPixelCol();
1339           int maxPixelRow = clust->maxPixelRow();
1340           int minPixelCol = clust->minPixelCol();
1341           int minPixelRow = clust->minPixelRow();
1342 
1343           //int geoId = PixGeom->geographicalId().rawId();
1344 
1345           // Replace with the topology methods
1346           // edge method moved to topologi class
1347           //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
1348           //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
1349 
1350           // calculate alpha and beta from cluster position
1351           //LocalTrajectoryParameters ltp = tsos.localParameters();
1352           //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
1353 
1354           //float locx = localDir.x();
1355           //float locy = localDir.y();
1356           //float locz = localDir.z();
1357           //float loctheta = localDir.theta(); // currently unused
1358 
1359           //float alpha = atan2( locz, locx );
1360           //float beta = atan2( locz, locy );
1361 
1362           //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));
1363         }  // valid cluster
1364       }    // valid peristant hit
1365 
1366     }  // loop over trajectory meas.
1367 
1368     if (PRINT)
1369       edm::LogPrint("TestWithTracks") << " Cluster for track " << trackNumber << " cluaters " << numberOfClusters << " "
1370                                       << pixelHits;
1371 
1372   }  // loop over tracks
1373 
1374 #endif  // USE_TRAJ
1375 
1376   edm::LogPrint("TestWithTracks") << " event with tracks = " << trackNumber << " " << countGoodTracks;
1377 
1378 }  // end
1379 
1380 //define this as a plug-in
1381 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!