Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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/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 "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0038 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0042 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0043 //#include "Geometry/CommonTopologies/interface/PixelTopology.h"
0044 
0045 // For L1
0046 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0047 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0048 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0049 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0050 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0051 
0052 // For HLT
0053 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0054 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0055 #include "DataFormats/Common/interface/TriggerResults.h"
0056 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0057 #include "FWCore/Common/interface/TriggerNames.h"
0058 
0059 // For tracks
0060 #include "DataFormats/TrackReco/interface/Track.h"
0061 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0062 
0063 //#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0064 
0065 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0066 
0067 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0068 //#include "TrackingTools/PatternTools/interface/TrajectoryFitter.h"
0069 
0070 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0071 
0072 //#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0073 //#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0074 
0075 //#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
0076 
0077 //#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0078 //#include "TrackingTools/PatternTools/interface/Trajectory.h"
0079 //#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0080 //#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0081 //#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0082 
0083 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0084 #include <DataFormats/VertexReco/interface/VertexFwd.h>
0085 
0086 // For luminisoty
0087 #include "DataFormats/Luminosity/interface/LumiSummary.h"
0088 #include "DataFormats/Common/interface/ConditionsInEdm.h"
0089 
0090 // To use root histos
0091 #include "FWCore/ServiceRegistry/interface/Service.h"
0092 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0093 
0094 // For ROOT
0095 #include <TROOT.h>
0096 #include <TChain.h>
0097 #include <TFile.h>
0098 #include <TF1.h>
0099 #include <TH2F.h>
0100 #include <TH1F.h>
0101 #include <TProfile.h>
0102 #include <TVector3.h>
0103 
0104 #define HISTOS
0105 //#define L1
0106 //#define HLT
0107 
0108 using namespace std;
0109 
0110 class TestPixTracks : public edm::EDAnalyzer {
0111 public:
0112   explicit TestPixTracks(const edm::ParameterSet &conf);
0113   virtual ~TestPixTracks();
0114   virtual void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0115   virtual void beginRun(const edm::EventSetup &iSetup) override;
0116   virtual void beginJob() override;
0117   virtual void endJob() override;
0118 
0119 private:
0120   edm::EDGetTokenT<LumiSummary> lumiToken_;
0121   edm::EDGetTokenT<edm::ConditionsInLumiBlock> condToken_;
0122   edm::EDGetTokenT<L1GlobalTriggerReadoutRecord> l1gtrrToken_;
0123   edm::EDGetTokenT<edm::TriggerResults> hltToken_;
0124   edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0125   edm::EDGetTokenT<reco::TrackCollection> srcToken_;
0126   edm::EDGetTokenT<TrajTrackAssociationCollection> trackAssocToken_;
0127   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0128   //const static bool PRINT = false;
0129   bool PRINT;
0130   float countTracks, countGoodTracks, countTracksInPix, countPVs, countEvents, countLumi;
0131 
0132   //TFile* hFile;
0133   //TH1D *hdetunit;
0134   //TH1D *hpixid,*hpixsubid,
0135   //*hlayerid,
0136   //*hladder1id,*hladder2id,*hladder3id,
0137   //*hz1id,*hz2id,*hz3id;
0138 
0139   TH1D *hcharge1, *hcharge2, *hcharge3, *hcharge;
0140   TH1D *hpixcharge1, *hpixcharge2, *hpixcharge3, *hpixcharge;
0141   TH1D *hcols1, *hcols2, *hcols3, *hrows1, *hrows2, *hrows3;
0142   TH1D *hsize1, *hsize2, *hsize3, *hsizex1, *hsizex2, *hsizex3, *hsizey1, *hsizey2, *hsizey3;
0143 
0144   TH1D *hclusPerTrk1, *hclusPerTrk2, *hclusPerTrk3;
0145   TH1D *hclusPerLay1, *hclusPerLay2, *hclusPerLay3;
0146   TH1D *hpixPerLay1, *hpixPerLay2, *hpixPerLay3;
0147   //TH1D *hdetsPerLay1,*hdetsPerLay2,*hdetsPerLay3;
0148 
0149   //TH1D *hdetr, *hdetz;
0150   //   TH1D *hcolsB,  *hrowsB,  *hcolsF,  *hrowsF;
0151   TH2F *hDetMap1, *hDetMap2, *hDetMap3;  // clusters
0152   //TH2F *hpixDetMap1, *hpixDetMap2, *hpixDetMap3;
0153   TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3;
0154 
0155   TH2F *hpvxy, *hclusMap1, *hclusMap2, *hclusMap3;
0156 
0157   TH1D *hpvz, *hpvr, *hNumPv, *hNumPvClean;
0158   TH1D *hPt, *hEta, *hDz, *hD0, *hzdiff;
0159 
0160   //TH1D *hncharge1,*hncharge2, *hncharge3;
0161   //TH1D *hchargeMonoPix1,*hchargeMonoPix2, *hchargeMonoPix3;
0162   // TH1D *hnpixcharge1,*hnpixcharge2, *hnpixcharge3;
0163   //TH1D *htest1,*htest2,*htest3,*htest4,*htest5,*htest6,*htest7,*htest8,*htest9;
0164   TH1D *hl1a, *hl1t, *hlt1;
0165 
0166   TH1D *hclusBpix, *hpixBpix;
0167   TH1D *htracks, *htracksGood, *htracksGoodInPix;
0168 
0169   TProfile *hclumult1, *hclumult2, *hclumult3;
0170   TProfile *hclumultx1, *hclumultx2, *hclumultx3;
0171   TProfile *hclumulty1, *hclumulty2, *hclumulty3;
0172   TProfile *hcluchar1, *hcluchar2, *hcluchar3;
0173   TProfile *hpixchar1, *hpixchar2, *hpixchar3;
0174 
0175   TProfile *htracksls, *hpvsls, *htrackslsn, *hpvslsn, *hintgl, *hinstl, *hbeam1, *hbeam2;
0176 
0177   TH1D *hlumi, *hlumi0, *hbx, *hbx0;
0178 };
0179 /////////////////////////////////////////////////////////////////
0180 // Contructor,
0181 TestPixTracks::TestPixTracks(edm::ParameterSet const &conf) {
0182   PRINT = conf.getUntrackedParameter<bool>("Verbosity", false);
0183   lumiToken_ = consumes<LumiSummary>(edm::InputTag("lumiProducer"));
0184   condToken_ = consumes<edm::ConditionsInLumiBlock>(edm::InputTag("conditionsInEdm"));
0185   l1gtrrToken_ = consumes<L1GlobalTriggerReadoutRecord>(edm::InputTag("gtDigis"));
0186   hltToken_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", "HLT"));
0187   vtxToken_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0188   srcToken_ = consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("src"));
0189   trackAssocToken_ =
0190       consumes<TrajTrackAssociationCollection>(edm::InputTag(conf.getParameter<std::string>("trajectoryInput")));
0191   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0192   //if(PRINT) cout<<" Construct "<<endl;
0193 }
0194 
0195 // Virtual destructor needed.
0196 TestPixTracks::~TestPixTracks() {}
0197 
0198 // ------------ method called at the begining   ------------
0199 void TestPixTracks::beginRun(const edm::EventSetup &iSetup) { cout << "BeginRun, Verbosity =  " << PRINT << endl; }
0200 
0201 // ------------ method called at the begining   ------------
0202 void TestPixTracks::beginJob() {
0203   cout << "BeginJob, Verbosity " << PRINT << endl;
0204 
0205   countTracks = 0.;
0206   countGoodTracks = 0.;
0207   countTracksInPix = 0.;
0208   countPVs = 0.;
0209   countEvents = 0.;
0210   countLumi = 0.;
0211 
0212 #ifdef HISTOS
0213 
0214   // NEW way to use root (from 2.0.0?)
0215   edm::Service<TFileService> fs;
0216 
0217   // put here whatever you want to do at the beginning of the job
0218   //hFile = new TFile ( "histo.root", "RECREATE" );
0219 
0220   //hladder1id = fs->make<TH1D>( "hladder1id", "Ladder L1 id", 50, 0., 50.);
0221   //hladder2id = fs->make<TH1D>( "hladder2id", "Ladder L2 id", 50, 0., 50.);
0222   //hladder3id = fs->make<TH1D>( "hladder3id", "Ladder L3 id", 50, 0., 50.);
0223   //hz1id = fs->make<TH1D>( "hz1id", "Z-index id L1", 10, 0., 10.);
0224   //hz2id = fs->make<TH1D>( "hz2id", "Z-index id L2", 10, 0., 10.);
0225   //hz3id = fs->make<TH1D>( "hz3id", "Z-index id L3", 10, 0., 10.);
0226 
0227   int sizeH = 20;
0228   float lowH = -0.5;
0229   float highH = 19.5;
0230 
0231   hclusPerTrk1 = fs->make<TH1D>("hclusPerTrk1", "Clus per track l1", sizeH, lowH, highH);
0232   hclusPerTrk2 = fs->make<TH1D>("hclusPerTrk2", "Clus per track l2", sizeH, lowH, highH);
0233   hclusPerTrk3 = fs->make<TH1D>("hclusPerTrk3", "Clus per track l3", sizeH, lowH, highH);
0234 
0235   sizeH = 2000;
0236   highH = 1999.5;
0237   hclusPerLay1 = fs->make<TH1D>("hclusPerLay1", "Clus per layer l1", sizeH, lowH, highH);
0238   hclusPerLay2 = fs->make<TH1D>("hclusPerLay2", "Clus per layer l2", sizeH, lowH, highH);
0239   hclusPerLay3 = fs->make<TH1D>("hclusPerLay3", "Clus per layer l3", sizeH, lowH, highH);
0240 
0241   highH = 9999.5;
0242   hpixPerLay1 = fs->make<TH1D>("hpixPerLay1", "Pix per layer l1", sizeH, lowH, highH);
0243   hpixPerLay2 = fs->make<TH1D>("hpixPerLay2", "Pix per layer l2", sizeH, lowH, highH);
0244   hpixPerLay3 = fs->make<TH1D>("hpixPerLay3", "Pix per layer l3", sizeH, lowH, highH);
0245 
0246   //hdetsPerLay1 = fs->make<TH1D>( "hdetsPerLay1", "Full dets per layer l1",
0247   //             161, -0.5, 160.5);
0248   //hdetsPerLay3 = fs->make<TH1D>( "hdetsPerLay3", "Full dets per layer l3",
0249   //             353, -0.5, 352.5);
0250   //hdetsPerLay2 = fs->make<TH1D>( "hdetsPerLay2", "Full dets per layer l2",
0251   //             257, -0.5, 256.5);
0252 
0253   hcharge1 = fs->make<TH1D>("hcharge1", "Clu charge l1", 400, 0., 200.);  //in ke
0254   hcharge2 = fs->make<TH1D>("hcharge2", "Clu charge l2", 400, 0., 200.);
0255   hcharge3 = fs->make<TH1D>("hcharge3", "Clu charge l3", 400, 0., 200.);
0256 
0257   //hchargeMonoPix1 = fs->make<TH1D>( "hchargeMonoPix1", "Clu charge l1 MonPix", 200, 0.,100.); //in ke
0258   //hchargeMonoPix2 = fs->make<TH1D>( "hchargeMonoPix2", "Clu charge l2 MonPix", 200, 0.,100.);
0259   //hchargeMonoPix3 = fs->make<TH1D>( "hchargeMonoPix3", "Clu charge l3 MonPix", 200, 0.,100.);
0260 
0261   //hncharge1 = fs->make<TH1D>( "hncharge1", "Noise charge l1", 200, 0.,100.); //in ke
0262   //hncharge2 = fs->make<TH1D>( "hncharge2", "Noise charge l2", 200, 0.,100.);
0263   //hncharge3 = fs->make<TH1D>( "hncharge3", "Noise charge l3", 200, 0.,100.);
0264 
0265   //hpixcharge1 = fs->make<TH1D>( "hpixcharge1", "Pix charge l1", 100, 0.,50.); //in ke
0266   //hpixcharge2 = fs->make<TH1D>( "hpixcharge2", "Pix charge l2", 100, 0.,50.);
0267   //hpixcharge3 = fs->make<TH1D>( "hpixcharge3", "Pix charge l3", 100, 0.,50.);
0268 
0269   //hnpixcharge1 = fs->make<TH1D>( "hnpixcharge1", "Noise pix charge l1", 100, 0.,50.); //in ke
0270   //hnpixcharge2 = fs->make<TH1D>( "hnpixcharge2", "Noise pix charge l2", 100, 0.,50.);
0271   //hnpixcharge3 = fs->make<TH1D>( "hnpixcharge3", "Noise pix charge l3", 100, 0.,50.);
0272 
0273   //hpixcharge = fs->make<TH1D>( "hpixcharge", "Clu charge", 100, 0.,50.);
0274   //hcharge = fs->make<TH1D>( "hcharge", "Pix charge", 100, 0.,50.);
0275 
0276   hcols1 = fs->make<TH1D>("hcols1", "Layer 1 cols", 500, -0.5, 499.5);
0277   hcols2 = fs->make<TH1D>("hcols2", "Layer 2 cols", 500, -0.5, 499.5);
0278   hcols3 = fs->make<TH1D>("hcols3", "Layer 3 cols", 500, -0.5, 499.5);
0279 
0280   hrows1 = fs->make<TH1D>("hrows1", "Layer 1 rows", 200, -0.5, 199.5);
0281   hrows2 = fs->make<TH1D>("hrows2", "Layer 2 rows", 200, -0.5, 199.5);
0282   hrows3 = fs->make<TH1D>("hrows3", "layer 3 rows", 200, -0.5, 199.5);
0283 
0284   hsize1 = fs->make<TH1D>("hsize1", "layer 1 clu size", 100, -0.5, 99.5);
0285   hsize2 = fs->make<TH1D>("hsize2", "layer 2 clu size", 100, -0.5, 99.5);
0286   hsize3 = fs->make<TH1D>("hsize3", "layer 3 clu size", 100, -0.5, 99.5);
0287   hsizex1 = fs->make<TH1D>("hsizex1", "lay1 clu size in x", 20, -0.5, 19.5);
0288   hsizex2 = fs->make<TH1D>("hsizex2", "lay2 clu size in x", 20, -0.5, 19.5);
0289   hsizex3 = fs->make<TH1D>("hsizex3", "lay3 clu size in x", 20, -0.5, 19.5);
0290   hsizey1 = fs->make<TH1D>("hsizey1", "lay1 clu size in y", 30, -0.5, 29.5);
0291   hsizey2 = fs->make<TH1D>("hsizey2", "lay2 clu size in y", 30, -0.5, 29.5);
0292   hsizey3 = fs->make<TH1D>("hsizey3", "lay3 clu size in y", 30, -0.5, 29.5);
0293 
0294   hDetMap1 = fs->make<TH2F>("hDetMap1", "layer 1 clus map", 9, 0., 9., 23, 0., 23.);
0295   hDetMap2 = fs->make<TH2F>("hDetMap2", "layer 2 clus map", 9, 0., 9., 33, 0., 33.);
0296   hDetMap3 = fs->make<TH2F>("hDetMap3", "layer 3 clus map", 9, 0., 9., 45, 0., 45.);
0297   //hpixDetMap1 = fs->make<TH2F>( "hpixDetMap1", "pix det layer 1",
0298   //          416,0.,416.,160,0.,160.);
0299   //hpixDetMap2 = fs->make<TH2F>( "hpixDetMap2", "pix det layer 2",
0300   //          416,0.,416.,160,0.,160.);
0301   //hpixDetMap3 = fs->make<TH2F>( "hpixDetMap3", "pix det layer 3",
0302   //          416,0.,416.,160,0.,160.);
0303   hcluDetMap1 = fs->make<TH2F>("hcluDetMap1", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0304   hcluDetMap2 = fs->make<TH2F>("hcluDetMap2", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0305   hcluDetMap3 = fs->make<TH2F>("hcluDetMap3", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0306 
0307   htracksGoodInPix = fs->make<TH1D>("htracksGoodInPix", "count good tracks in pix", 2000, -0.5, 1999.5);
0308   htracksGood = fs->make<TH1D>("htracksGood", "count good tracks", 2000, -0.5, 1999.5);
0309   htracks = fs->make<TH1D>("htracks", "count tracks", 2000, -0.5, 1999.5);
0310   hclusBpix = fs->make<TH1D>("hclusBpix", "count clus in bpix", 200, -0.5, 1999.5);
0311   hpixBpix = fs->make<TH1D>("hpixBpix", "count pixels", 200, -0.5, 1999.5);
0312 
0313   hpvxy = fs->make<TH2F>("hpvxy", "pv xy", 100, -1., 1., 100, -1., 1.);
0314   hpvz = fs->make<TH1D>("hpvz", "pv z", 1000, -50., 50.);
0315   hpvr = fs->make<TH1D>("hpvr", "pv r", 100, 0., 1.);
0316   hNumPv = fs->make<TH1D>("hNumPv", "num of pv", 100, 0., 100.);
0317   hNumPvClean = fs->make<TH1D>("hNumPvClean", "num of pv clean", 100, 0., 100.);
0318 
0319   hPt = fs->make<TH1D>("hPt", "pt", 100, 0., 100.);
0320   hEta = fs->make<TH1D>("hEta", "eta", 50, -2.5, 2.5);
0321   hD0 = fs->make<TH1D>("hD0", "d0", 500, 0., 5.);
0322   hDz = fs->make<TH1D>("hDz", "pt", 250, -25., 25.);
0323   hzdiff = fs->make<TH1D>("hzdiff", "PVz-Trackz", 200, -10., 10.);
0324 
0325   hl1a = fs->make<TH1D>("hl1a", "l1a", 128, -0.5, 127.5);
0326   hl1t = fs->make<TH1D>("hl1t", "l1t", 128, -0.5, 127.5);
0327   hlt1 = fs->make<TH1D>("hlt1", "hlt1", 256, -0.5, 255.5);
0328 
0329   hclumult1 = fs->make<TProfile>("hclumult1", "cluster size layer 1", 60, -3., 3., 0.0, 100.);
0330   hclumult2 = fs->make<TProfile>("hclumult2", "cluster size layer 2", 60, -3., 3., 0.0, 100.);
0331   hclumult3 = fs->make<TProfile>("hclumult3", "cluster size layer 3", 60, -3., 3., 0.0, 100.);
0332 
0333   hclumultx1 = fs->make<TProfile>("hclumultx1", "cluster x-size layer 1", 60, -3., 3., 0.0, 100.);
0334   hclumultx2 = fs->make<TProfile>("hclumultx2", "cluster x-size layer 2", 60, -3., 3., 0.0, 100.);
0335   hclumultx3 = fs->make<TProfile>("hclumultx3", "cluster x-size layer 3", 60, -3., 3., 0.0, 100.);
0336 
0337   hclumulty1 = fs->make<TProfile>("hclumulty1", "cluster y-size layer 1", 60, -3., 3., 0.0, 100.);
0338   hclumulty2 = fs->make<TProfile>("hclumulty2", "cluster y-size layer 2", 60, -3., 3., 0.0, 100.);
0339   hclumulty3 = fs->make<TProfile>("hclumulty3", "cluster y-size layer 3", 60, -3., 3., 0.0, 100.);
0340 
0341   hcluchar1 = fs->make<TProfile>("hcluchar1", "cluster char layer 1", 60, -3., 3., 0.0, 1000.);
0342   hcluchar2 = fs->make<TProfile>("hcluchar2", "cluster char layer 2", 60, -3., 3., 0.0, 1000.);
0343   hcluchar3 = fs->make<TProfile>("hcluchar3", "cluster char layer 3", 60, -3., 3., 0.0, 1000.);
0344 
0345   hpixchar1 = fs->make<TProfile>("hpixchar1", "pix char layer 1", 60, -3., 3., 0.0, 1000.);
0346   hpixchar2 = fs->make<TProfile>("hpixchar2", "pix char layer 2", 60, -3., 3., 0.0, 1000.);
0347   hpixchar3 = fs->make<TProfile>("hpixchar3", "pix char layer 3", 60, -3., 3., 0.0, 1000.);
0348 
0349   hintgl = fs->make<TProfile>("hintgl", "inst lumi vs ls ", 1000, 0., 3000., 0.0, 10000.);
0350   hinstl = fs->make<TProfile>("hinstl", "intg lumi vs ls ", 1000, 0., 3000., 0.0, 100.);
0351   hbeam1 = fs->make<TProfile>("hbeam1", "beam1 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0352   hbeam2 = fs->make<TProfile>("hbeam2", "beam2 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0353 
0354   htracksls = fs->make<TProfile>("htracksls", "tracks with pix hits  vs ls", 1000, 0., 3000., 0.0, 10000.);
0355   hpvsls = fs->make<TProfile>("hpvsls", "pvs  vs ls", 1000, 0., 3000., 0.0, 1000.);
0356   htrackslsn = fs->make<TProfile>("htrackslsn", "tracks with pix hits/lumi  vs ls", 1000, 0., 3000., 0.0, 10000.);
0357   hpvslsn = fs->make<TProfile>("hpvslsn", "pvs/lumi  vs ls", 1000, 0., 3000., 0.0, 1000.);
0358 
0359   hlumi0 = fs->make<TH1D>("hlumi0", "lumi", 2000, 0, 2000.);
0360   hlumi = fs->make<TH1D>("hlumi", "lumi", 2000, 0, 2000.);
0361   hbx0 = fs->make<TH1D>("hbx0", "bx", 4000, 0, 4000.);
0362   hbx = fs->make<TH1D>("hbx", "bx", 4000, 0, 4000.);
0363 
0364   hclusMap1 = fs->make<TH2F>("hclusMap1", "clus - lay1", 260, -26., 26., 350, -3.5, 3.5);
0365   hclusMap2 = fs->make<TH2F>("hclusMap2", "clus - lay2", 260, -26., 26., 350, -3.5, 3.5);
0366   hclusMap3 = fs->make<TH2F>("hclusMap3", "clus - lay3", 260, -26., 26., 350, -3.5, 3.5);
0367 
0368 #endif
0369 }
0370 // ------------ method called to at the end of the job  ------------
0371 void TestPixTracks::endJob() {
0372   cout << " End PixelTracksTest " << endl;
0373 
0374   if (countEvents > 0.) {
0375     countTracks /= countEvents;
0376     countGoodTracks /= countEvents;
0377     countTracksInPix /= countEvents;
0378     countPVs /= countEvents;
0379     countLumi /= 1000.;
0380     cout << " Average tracks/event " << countTracks << " good " << countGoodTracks << " in pix " << countTracksInPix
0381          << " PVs " << countPVs << " events " << countEvents << " lumi pb-1 " << countLumi << "/10, bug!" << endl;
0382   }
0383 }
0384 //////////////////////////////////////////////////////////////////
0385 // Functions that gets called by framework every event
0386 void TestPixTracks::analyze(const edm::Event &e, const edm::EventSetup &es) {
0387   using namespace edm;
0388   using namespace reco;
0389   static int lumiBlockOld = -9999;
0390 
0391   const float CLU_SIZE_PT_CUT = 1.;
0392 
0393   int trackNumber = 0;
0394   int countNiceTracks = 0;
0395   int countPixTracks = 0;
0396 
0397   int numberOfDetUnits = 0;
0398   int numberOfClusters = 0;
0399   int numberOfPixels = 0;
0400 
0401   int numberOfDetUnits1 = 0;
0402   int numOfClustersPerDet1 = 0;
0403   int numOfClustersPerLay1 = 0;
0404   int numOfPixelsPerLay1 = 0;
0405 
0406   int numberOfDetUnits2 = 0;
0407   int numOfClustersPerDet2 = 0;
0408   int numOfClustersPerLay2 = 0;
0409   int numOfPixelsPerLay2 = 0;
0410 
0411   int numberOfDetUnits3 = 0;
0412   int numOfClustersPerDet3 = 0;
0413   int numOfClustersPerLay3 = 0;
0414   int numOfPixelsPerLay3 = 0;
0415 
0416   int run = e.id().run();
0417   int event = e.id().event();
0418   int lumiBlock = e.luminosityBlock();
0419   int bx = e.bunchCrossing();
0420   int orbit = e.orbitNumber();
0421 
0422   if (PRINT)
0423     cout << "Run " << run << " Event " << event << " LS " << lumiBlock << endl;
0424 
0425   hbx0->Fill(float(bx));
0426   hlumi0->Fill(float(lumiBlock));
0427 
0428   edm::LuminosityBlock const &iLumi = e.getLuminosityBlock();
0429   edm::Handle<LumiSummary> lumi;
0430   iLumi.getByToken(lumiToken_, lumi);
0431   edm::Handle<edm::ConditionsInLumiBlock> cond;
0432   float intlumi = 0, instlumi = 0;
0433   int beamint1 = 0, beamint2 = 0;
0434   iLumi.getByToken(condToken_, cond);
0435   // This will only work when running on RECO until (if) they fix it in the FW
0436   // When running on RAW and reconstructing, the LumiSummary will not appear
0437   // in the event before reaching endLuminosityBlock(). Therefore, it is not
0438   // possible to get this info in the event
0439   if (lumi.isValid()) {
0440     intlumi = (lumi->intgRecLumi()) / 1000.;  // 10^30 -> 10^33/cm2/sec  ->  1/nb/sec
0441     instlumi = (lumi->avgInsDelLumi()) / 1000.;
0442     beamint1 = (cond->totalIntensityBeam1) / 1000;
0443     beamint2 = (cond->totalIntensityBeam2) / 1000;
0444   } else {
0445     //std::cout << "** ERROR: Event does not get lumi info\n";
0446   }
0447 
0448   hinstl->Fill(float(lumiBlock), float(instlumi));
0449   hintgl->Fill(float(lumiBlock), float(intlumi));
0450   hbeam1->Fill(float(lumiBlock), float(beamint1));
0451   hbeam2->Fill(float(lumiBlock), float(beamint2));
0452 
0453 #ifdef L1
0454   // Get L1
0455   Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
0456   e.getByToken(l1gtrrToken_, L1GTRR);
0457 
0458   if (L1GTRR.isValid()) {
0459     //bool l1a = L1GTRR->decision();  // global decission?
0460     //cout<<" L1 status :"<<l1a<<" "<<hex;
0461     for (unsigned int i = 0; i < L1GTRR->decisionWord().size(); ++i) {
0462       int l1flag = L1GTRR->decisionWord()[i];
0463       int t1flag = L1GTRR->technicalTriggerWord()[i];
0464 
0465       if (l1flag > 0)
0466         hl1a->Fill(float(i));
0467       if (t1flag > 0 && i < 64)
0468         hl1t->Fill(float(i));
0469     }  // for loop
0470   }    // if l1a
0471 #endif
0472 
0473 #ifdef HLT
0474 
0475   bool hlt[256];
0476   for (int i = 0; i < 256; ++i)
0477     hlt[i] = false;
0478 
0479   edm::TriggerNames TrigNames;
0480   edm::Handle<edm::TriggerResults> HLTResults;
0481 
0482   // Extract the HLT results
0483   e.getByToken(hltToken_, HLTResults);
0484   if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) {
0485     //TrigNames.init(*HLTResults);
0486     const edm::TriggerNames &TrigNames = e.triggerNames(*HLTResults);
0487 
0488     //cout<<TrigNames.triggerNames().size()<<endl;
0489 
0490     for (unsigned int i = 0; i < TrigNames.triggerNames().size(); i++) {  // loop over trigger
0491       //if(countAllEvents==1) cout<<i<<" "<<TrigNames.triggerName(i)<<endl;
0492 
0493       if ((HLTResults->wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0494           (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0495           (HLTResults->error(TrigNames.triggerIndex(TrigNames.triggerName(i))) == false)) {
0496         hlt[i] = true;
0497         hlt1->Fill(float(i));
0498 
0499       }  // if hlt
0500 
0501     }  // loop
0502   }    // if valid
0503 #endif
0504 
0505   // Get event setup
0506   edm::ESHandle<TrackerGeometry> geom = es.getHandle(trackerGeomToken_);
0507   const TrackerGeometry &theTracker(*geom);
0508 
0509   // -- Primary vertices
0510   // ----------------------------------------------------------------------
0511   edm::Handle<reco::VertexCollection> vertices;
0512   e.getByToken(vtxToken_, vertices);
0513 
0514   if (PRINT)
0515     cout << " PV list " << vertices->size() << endl;
0516   int pvNotFake = 0, pvsTrue = 0;
0517   vector<float> pvzVector;
0518   for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) {
0519     if ((iv->isFake()) == 1)
0520       continue;
0521     pvNotFake++;
0522     float pvx = iv->x();
0523     float pvy = iv->y();
0524     float pvz = iv->z();
0525     int numTracksPerPV = iv->tracksSize();
0526     //int numTracksPerPV = iv->nTracks();
0527 
0528     //float xe = iv->xError();
0529     //float ye = iv->yError();
0530     //float ze = iv->zError();
0531     //int chi2 = iv->chi2();
0532     //int dof = iv->ndof();
0533 
0534     if (PRINT)
0535       cout << " PV " << pvNotFake << " pos = " << pvx << "/" << pvy << "/" << pvz << ", Num of tracks "
0536            << numTracksPerPV << endl;
0537 
0538     hpvz->Fill(pvz);
0539     if (pvz > -22. && pvz < 22.) {
0540       float pvr = sqrt(pvx * pvx + pvy * pvy);
0541       hpvxy->Fill(pvx, pvy);
0542       hpvr->Fill(pvr);
0543       if (pvr < 0.3) {
0544         pvsTrue++;
0545         pvzVector.push_back(pvz);
0546         //if(PRINT) cout<<"PV "<<pvsTrue<<" "<<pvz<<endl;
0547       }  //pvr
0548     }    // pvz
0549 
0550     //if(pvsTrue<1) continue; // skip events with no PV
0551 
0552   }  // loop pvs
0553   hNumPv->Fill(float(pvNotFake));
0554   hNumPvClean->Fill(float(pvsTrue));
0555 
0556   if (PRINT)
0557     cout << " Not fake PVs = " << pvNotFake << " good position " << pvsTrue << endl;
0558 
0559   Handle<reco::TrackCollection> recTracks;
0560   e.getByToken(srcToken_, recTracks);
0561 
0562   if (PRINT)
0563     cout << " Tracks " << recTracks->size() << endl;
0564   for (reco::TrackCollection::const_iterator t = recTracks->begin(); t != recTracks->end(); ++t) {
0565     trackNumber++;
0566     numOfClustersPerDet1 = 0;
0567     numOfClustersPerDet2 = 0;
0568     numOfClustersPerDet3 = 0;
0569     int pixelHits = 0;
0570 
0571     int size = t->recHitsSize();
0572     float pt = t->pt();
0573     float eta = t->eta();
0574     float phi = t->phi();
0575     float trackCharge = t->charge();
0576     float d0 = t->d0();
0577     float dz = t->dz();
0578     float tkvx = t->vx();
0579     float tkvy = t->vy();
0580     float tkvz = t->vz();
0581 
0582     if (PRINT)
0583       cout << "Track " << trackNumber << " Pt " << pt << " Eta " << eta << " d0/dz " << d0 << " " << dz << " Hits "
0584            << size << endl;
0585 
0586     hEta->Fill(eta);
0587     hDz->Fill(dz);
0588     if (abs(eta) > 2.8 || abs(dz) > 25.)
0589       continue;  //  skip
0590 
0591     hD0->Fill(d0);
0592     if (d0 > 1.0)
0593       continue;  // skip
0594 
0595     bool goodTrack = false;
0596     for (vector<float>::iterator m = pvzVector.begin(); m != pvzVector.end(); ++m) {
0597       float z = *m;
0598       float tmp = abs(z - dz);
0599       hzdiff->Fill(tmp);
0600       if (tmp < 1.)
0601         goodTrack = true;
0602     }
0603 
0604     if (!goodTrack)
0605       continue;
0606     countNiceTracks++;
0607     hPt->Fill(pt);
0608 
0609     // Loop over rechits
0610     for (trackingRecHit_iterator recHit = t->recHitsBegin(); recHit != t->recHitsEnd(); ++recHit) {
0611       if (!((*recHit)->isValid()))
0612         continue;
0613 
0614       if ((*recHit)->geographicalId().det() != DetId::Tracker)
0615         continue;
0616 
0617       const DetId &hit_detId = (*recHit)->geographicalId();
0618       uint IntSubDetID = (hit_detId.subdetId());
0619 
0620       // Select pixel rechits
0621       if (IntSubDetID == 0)
0622         continue;  // Select ??
0623       if (IntSubDetID != PixelSubdetector::PixelBarrel)
0624         continue;  // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) {
0625 
0626       // Pixel detector
0627       PXBDetId pdetId = PXBDetId(hit_detId);
0628       //unsigned int detTypeP=pdetId.det();
0629       //unsigned int subidP=pdetId.subdetId();
0630       // Barell layer = 1,2,3
0631       int layer = pdetId.layer();
0632       // Barrel ladder id 1-20,32,44.
0633       int ladder = pdetId.ladder();
0634       // Barrel Z-index=1,8
0635       int zindex = pdetId.module();
0636 
0637       if (PRINT)
0638         cout << "barrel layer/ladder/module: " << layer << "/" << ladder << "/" << zindex << endl;
0639 
0640       // Get the geom-detector
0641       const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
0642       double detZ = theGeomDet->surface().position().z();
0643       double detR = theGeomDet->surface().position().perp();
0644       const PixelTopology *topol = &(theGeomDet->specificTopology());  // pixel topology
0645 
0646       //std::vector<SiStripRecHit2D*> output = getRecHitComponents((*recHit).get());
0647       //std::vector<SiPixelRecHit*> TrkComparison::getRecHitComponents(const TrackingRecHit* rechit){
0648 
0649       const SiPixelRecHit *hit = dynamic_cast<const SiPixelRecHit *>((*recHit).get());
0650       //edm::Ref<edmNew::DetSetVector<SiStripCluster> ,SiStripCluster> cluster = hit->cluster();
0651       // get the edm::Ref to the cluster
0652 
0653       if (hit) {
0654         // RecHit (recthits are transient, so not available without ttrack refit)
0655         //double xloc = hit->localPosition().x();// 1st meas coord
0656         //double yloc = hit->localPosition().y();// 2nd meas coord or zero
0657         //double zloc = hit->localPosition().z();// up, always zero
0658         //cout<<" rechit loc "<<xloc<<" "<<yloc<<endl;
0659 
0660         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = hit->cluster();
0661         //  check if the ref is not null
0662         if (!clust.isNonnull())
0663           continue;
0664 
0665         numberOfClusters++;
0666         pixelHits++;
0667         float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
0668         int size = clust->size();
0669         int sizeX = clust->sizeX();
0670         int sizeY = clust->sizeY();
0671         float row = clust->x();
0672         float col = clust->y();
0673         numberOfPixels += size;
0674 
0675         //cout<<" clus loc "<<row<<" "<<col<<endl;
0676 
0677         if (PRINT)
0678           cout << " cluster " << numberOfClusters << " charge = " << charge << " size = " << size << endl;
0679 
0680         LocalPoint lp = topol->localPosition(MeasurementPoint(clust->x(), clust->y()));
0681         //float x = lp.x();
0682         //float y = lp.y();
0683         //cout<<" clu loc "<<x<<" "<<y<<endl;
0684 
0685         GlobalPoint clustgp = theGeomDet->surface().toGlobal(lp);
0686         double gX = clustgp.x();
0687         double gY = clustgp.y();
0688         double gZ = clustgp.z();
0689 
0690         //cout<<" CLU GLOBAL "<<gX<<" "<<gY<<" "<<gZ<<endl;
0691 
0692         TVector3 v(gX, gY, gZ);
0693         float phi = v.Phi();
0694 
0695         //int maxPixelCol = clust->maxPixelCol();
0696         //int maxPixelRow = clust->maxPixelRow();
0697         //int minPixelCol = clust->minPixelCol();
0698         //int minPixelRow = clust->minPixelRow();
0699         //int geoId = PixGeom->geographicalId().rawId();
0700         // Replace with the topology methods
0701         // edge method moved to topologi class
0702         //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
0703         //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
0704 
0705         // calculate alpha and beta from cluster position
0706         //LocalTrajectoryParameters ltp = tsos.localParameters();
0707         //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0708         //float locx = localDir.x();
0709         //float locy = localDir.y();
0710         //float locz = localDir.z();
0711         //float loctheta = localDir.theta(); // currently unused
0712         //float alpha = atan2( locz, locx );
0713         //float beta = atan2( locz, locy );
0714 
0715         if (layer == 1) {
0716           hDetMap1->Fill(float(zindex), float(ladder));
0717           hcluDetMap1->Fill(col, row);
0718           hcharge1->Fill(charge);
0719           hcols1->Fill(col);
0720           hrows1->Fill(row);
0721 
0722           hclusMap1->Fill(gZ, phi);
0723 
0724           if (pt > CLU_SIZE_PT_CUT) {
0725             hsize1->Fill(float(size));
0726             hsizex1->Fill(float(sizeX));
0727             hsizey1->Fill(float(sizeY));
0728 
0729             hclumult1->Fill(eta, float(size));
0730             hclumultx1->Fill(eta, float(sizeX));
0731             hclumulty1->Fill(eta, float(sizeY));
0732             hcluchar1->Fill(eta, float(charge));
0733           }
0734 
0735           numOfClustersPerDet1++;
0736           numOfClustersPerLay1++;
0737           numOfPixelsPerLay1 += size;
0738 
0739         } else if (layer == 2) {
0740           hDetMap2->Fill(float(zindex), float(ladder));
0741           hcluDetMap2->Fill(col, row);
0742           hcharge2->Fill(charge);
0743           hcols2->Fill(col);
0744           hrows2->Fill(row);
0745 
0746           hclusMap2->Fill(gZ, phi);
0747 
0748           if (pt > CLU_SIZE_PT_CUT) {
0749             hsize2->Fill(float(size));
0750             hsizex2->Fill(float(sizeX));
0751             hsizey2->Fill(float(sizeY));
0752 
0753             hclumult2->Fill(eta, float(size));
0754             hclumultx2->Fill(eta, float(sizeX));
0755             hclumulty2->Fill(eta, float(sizeY));
0756             hcluchar2->Fill(eta, float(charge));
0757           }
0758 
0759           numOfClustersPerDet2++;
0760           numOfClustersPerLay2++;
0761           numOfPixelsPerLay2 += size;
0762 
0763         } else if (layer == 3) {
0764           hDetMap3->Fill(float(zindex), float(ladder));
0765           hcluDetMap3->Fill(col, row);
0766           hcharge3->Fill(charge);
0767           hcols3->Fill(col);
0768           hrows3->Fill(row);
0769 
0770           hclusMap3->Fill(gZ, phi);
0771 
0772           if (pt > CLU_SIZE_PT_CUT) {
0773             hsize3->Fill(float(size));
0774             hsizex3->Fill(float(sizeX));
0775             hsizey3->Fill(float(sizeY));
0776             hclumult3->Fill(eta, float(size));
0777             hclumultx3->Fill(eta, float(sizeX));
0778             hclumulty3->Fill(eta, float(sizeY));
0779             hcluchar3->Fill(eta, float(charge));
0780           }
0781 
0782           numOfClustersPerDet3++;
0783           numOfClustersPerLay3++;
0784           numOfPixelsPerLay3 += size;
0785         }  // if layer
0786 
0787       }  // if valid
0788 
0789     }  // clusters
0790 
0791     if (pixelHits > 0)
0792       countPixTracks++;
0793 
0794     if (PRINT)
0795       cout << " Clusters for track " << trackNumber << " num of clusters " << numberOfClusters << " num of pixels "
0796            << pixelHits << endl;
0797 
0798 #ifdef HISTOS
0799     if (numberOfClusters > 0) {
0800       hclusPerTrk1->Fill(float(numOfClustersPerDet1));
0801       if (PRINT)
0802         cout << "Lay1: number of clusters per track = " << numOfClustersPerDet1 << endl;
0803       hclusPerTrk2->Fill(float(numOfClustersPerDet2));
0804       if (PRINT)
0805         cout << "Lay2: number of clusters per track = " << numOfClustersPerDet1 << endl;
0806       hclusPerTrk3->Fill(float(numOfClustersPerDet3));
0807       if (PRINT)
0808         cout << "Lay3: number of clusters per track = " << numOfClustersPerDet1 << endl;
0809     }
0810 #endif  // HISTOS
0811 
0812   }  // tracks
0813 
0814 #ifdef HISTOS
0815   if (numberOfClusters > 0) {
0816     hclusPerLay1->Fill(float(numOfClustersPerLay1));
0817     hclusPerLay2->Fill(float(numOfClustersPerLay2));
0818     hclusPerLay3->Fill(float(numOfClustersPerLay3));
0819     //hdetsPerLay1->Fill(float(numberOfDetUnits1));
0820     //hdetsPerLay2->Fill(float(numberOfDetUnits2));
0821     //hdetsPerLay3->Fill(float(numberOfDetUnits3));
0822     //int tmp = numberOfDetUnits1+numberOfDetUnits2+numberOfDetUnits3;
0823     hpixPerLay1->Fill(float(numOfPixelsPerLay1));
0824     hpixPerLay2->Fill(float(numOfPixelsPerLay2));
0825     hpixPerLay3->Fill(float(numOfPixelsPerLay3));
0826     //htest7->Fill(float(tmp));
0827     hclusBpix->Fill(float(numberOfClusters));
0828     hpixBpix->Fill(float(numberOfPixels));
0829   }
0830   htracksGood->Fill(float(countNiceTracks));
0831   htracksGoodInPix->Fill(float(countPixTracks));
0832   htracks->Fill(float(trackNumber));
0833 
0834   hbx->Fill(float(bx));
0835   hlumi->Fill(float(lumiBlock));
0836 
0837   htracksls->Fill(float(lumiBlock), float(countPixTracks));
0838   hpvsls->Fill(float(lumiBlock), float(pvsTrue));
0839   if (instlumi > 0.) {
0840     float tmp = float(countPixTracks) / instlumi;
0841     htrackslsn->Fill(float(lumiBlock), tmp);
0842     tmp = float(pvsTrue) / instlumi;
0843     hpvslsn->Fill(float(lumiBlock), tmp);
0844   }
0845 
0846 #endif  // HISTOS
0847 
0848   //
0849   countTracks += float(trackNumber);
0850   countGoodTracks += float(countNiceTracks);
0851   countTracksInPix += float(countPixTracks);
0852   countPVs += float(pvsTrue);
0853   countEvents++;
0854   if (lumiBlock != lumiBlockOld) {
0855     countLumi += intlumi;
0856     lumiBlockOld = lumiBlock;
0857   }
0858 
0859   if (PRINT)
0860     cout << " event with tracks = " << trackNumber << " " << countNiceTracks << endl;
0861 
0862   return;
0863 
0864 #ifdef USE_TRAJ
0865 
0866   //------------------------------------------------------------------------------------
0867   // Use Trajectories
0868 
0869   edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
0870   e.getByToken(trackAssocToken_, trajTrackCollectionHandle);
0871 
0872   TrajectoryStateCombiner tsoscomb;
0873 
0874   int NbrTracks = trajTrackCollectionHandle->size();
0875   std::cout << " track measurements " << trajTrackCollectionHandle->size() << std::endl;
0876 
0877   int trackNumber = 0;
0878   int numberOfClusters = 0;
0879 
0880   for (TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(),
0881                                                       itEnd = trajTrackCollectionHandle->end();
0882        it != itEnd;
0883        ++it) {
0884     int pixelHits = 0;
0885     int stripHits = 0;
0886     const Track &track = *it->val;
0887     const Trajectory &traj = *it->key;
0888 
0889     std::vector<TrajectoryMeasurement> checkColl = traj.measurements();
0890     for (std::vector<TrajectoryMeasurement>::const_iterator checkTraj = checkColl.begin(); checkTraj != checkColl.end();
0891          ++checkTraj) {
0892       if (!checkTraj->updatedState().isValid())
0893         continue;
0894       TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit();
0895       if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
0896         continue;
0897       uint testSubDetID = (testhit->geographicalId().subdetId());
0898       if (testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap)
0899         pixelHits++;
0900       else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB ||
0901                testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC)
0902         stripHits++;
0903     }
0904 
0905     if (pixelHits == 0)
0906       continue;
0907 
0908     trackNumber++;
0909     std::cout << " track " << trackNumber << " has pixelhits " << pixelHits << std::endl;
0910     pixelHits = 0;
0911 
0912     //std::vector<TrajectoryMeasurement> tmColl = traj.measurements();
0913     for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = checkColl.begin(); itTraj != checkColl.end();
0914          ++itTraj) {
0915       if (!itTraj->updatedState().isValid())
0916         continue;
0917 
0918       TrajectoryStateOnSurface tsos = tsoscomb(itTraj->forwardPredictedState(), itTraj->backwardPredictedState());
0919       TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
0920       if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker)
0921         continue;
0922 
0923       const DetId &hit_detId = hit->geographicalId();
0924       uint IntSubDetID = (hit_detId.subdetId());
0925 
0926       if (IntSubDetID == 0)
0927         continue;  // Select ??
0928       if (IntSubDetID != PixelSubdetector::PixelBarrel)
0929         continue;  // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) {
0930 
0931       //       const GeomDetUnit* detUnit = hit->detUnit();
0932       //       if(detUnit) {
0933       //    const Surface& surface = hit->detUnit()->surface();
0934       //    const TrackerGeometry& theTracker(*tkGeom_);
0935       //    const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*> (theTracker.idToDet(hit_detId) );
0936       //    const RectangularPixelTopology * topol = dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
0937       //       }
0938 
0939       // get the enclosed persistent hit
0940       const TrackingRecHit *persistentHit = hit->hit();
0941       // check if it's not null, and if it's a valid pixel hit
0942       if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
0943         // tell the C++ compiler that the hit is a pixel hit
0944         const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
0945         // get the edm::Ref to the cluster
0946         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
0947         //  check if the ref is not null
0948         if (clust.isNonnull()) {
0949           numberOfClusters++;
0950           pixelHits++;
0951           float charge = (clust->charge()) / 1000.0;  // convert electrons to kilo-electrons
0952           int size = clust->size();
0953           int size_x = clust->sizeX();
0954           int size_y = clust->sizeY();
0955           float row = clust->x();
0956           float col = clust->y();
0957 
0958           //LocalPoint lp = topol->localPosition(MeasurementPoint(clust_.row,clust_.col));
0959           //float x = lp.x();
0960           //float y = lp.y();
0961 
0962           int maxPixelCol = clust->maxPixelCol();
0963           int maxPixelRow = clust->maxPixelRow();
0964           int minPixelCol = clust->minPixelCol();
0965           int minPixelRow = clust->minPixelRow();
0966 
0967           //int geoId = PixGeom->geographicalId().rawId();
0968 
0969           // Replace with the topology methods
0970           // edge method moved to topologi class
0971           //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) );
0972           //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) );
0973 
0974           // calculate alpha and beta from cluster position
0975           //LocalTrajectoryParameters ltp = tsos.localParameters();
0976           //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
0977 
0978           //float locx = localDir.x();
0979           //float locy = localDir.y();
0980           //float locz = localDir.z();
0981           //float loctheta = localDir.theta(); // currently unused
0982 
0983           //float alpha = atan2( locz, locx );
0984           //float beta = atan2( locz, locy );
0985 
0986           //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));
0987         }  // valid cluster
0988       }    // valid peristant hit
0989 
0990     }  // loop over trajectory meas.
0991 
0992     if (PRINT)
0993       cout << " Cluster for track " << trackNumber << " cluaters " << numberOfClusters << " " << pixelHits << endl;
0994 
0995   }  // loop over tracks
0996 
0997 #endif  // USE_TRAJ
0998 
0999   cout << " event with tracks = " << trackNumber << " " << countGoodTracks << endl;
1000 
1001 }  // end
1002 
1003 //define this as a plug-in
1004 DEFINE_FWK_MODULE(TestPixTracks);
[ 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!