diff --git a/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader.cc b/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader.cc index 0d062e0587dc5..3d19608288a21 100644 --- a/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader.cc +++ b/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader.cc @@ -14,11 +14,13 @@ #include "DataFormats/SiPixelDetId/interface/PixelEndcapName.h" #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" + #include #include #include - using namespace cms; SiPixelLorentzAngleReader::SiPixelLorentzAngleReader( const edm::ParameterSet& iConfig ): @@ -31,25 +33,79 @@ SiPixelLorentzAngleReader::~SiPixelLorentzAngleReader(){} void SiPixelLorentzAngleReader::analyze( const edm::Event& e, const edm::EventSetup& iSetup){ edm::ESHandle SiPixelLorentzAngle_; - if(useSimRcd_ == true) - iSetup.get().get(SiPixelLorentzAngle_); - else - iSetup.get().get(SiPixelLorentzAngle_); - edm::LogInfo("SiPixelLorentzAngleReader") << "[SiPixelLorentzAngleReader::analyze] End Reading SiPixelLorentzAngle" << std::endl; - edm::Service fs; - LorentzAngleBarrel_ = fs->make("LorentzAngleBarrelPixel","LorentzAngleBarrelPixel",150,0,0.15); - LorentzAngleForward_= fs->make("LorentzAngleForwardPixel","LorentzAngleForwardPixel",150,0,0.15); - std::map detid_la= SiPixelLorentzAngle_->getLorentzAngles(); - std::map::const_iterator it; - for (it=detid_la.begin();it!=detid_la.end();it++) - { - // std::cout << "detid " << it->first << " \t" << " Lorentz angle " << it->second << std::endl; - //edm::LogInfo("SiPixelLorentzAngleReader") << "detid " << it->first << " \t" << " Lorentz angle " << it->second; - unsigned int subdet = DetId(it->first).subdetId(); - if(subdet == static_cast(PixelSubdetector::PixelBarrel)){ - LorentzAngleBarrel_->Fill(it->second); - }else if(subdet == static_cast(PixelSubdetector::PixelEndcap)){ - LorentzAngleForward_->Fill(it->second); - } - } + + if(useSimRcd_ == true) { + iSetup.get().get(SiPixelLorentzAngle_); + edm::LogInfo("SiPixelLorentzAngleReader") <<" Show LA for reconstruction "<().get(SiPixelLorentzAngle_); + //iSetup.get().get("fromAlignment",SiPixelLorentzAngle_); + //iSetup.get().get("forWidth",SiPixelLorentzAngle_); + edm::LogInfo("SiPixelLorentzAngleReader") <<" Show LA for reconstruction "< fs; + LorentzAngleBarrel_ = fs->make("LorentzAngleBarrelPixel","LorentzAngleBarrelPixel",150,0,0.15); + LorentzAngleForward_= fs->make("LorentzAngleForwardPixel","LorentzAngleForwardPixel",150,0,0.15); + std::map detid_la= SiPixelLorentzAngle_->getLorentzAngles(); + std::map::const_iterator it; + double la_old=-1.; + for (it=detid_la.begin();it!=detid_la.end();it++) { + //if(printdebug_) std::cout << "detid " << it->first << " \t" << " Lorentz angle " << it->second << std::endl; + //if(printdebug_) edm::LogInfo("SiPixelLorentzAngleReader") << "detid " << it->first << " \t" << " Lorentz angle " << it->second; + + unsigned int subdet = DetId(it->first).subdetId(); + int detid = it->first; + + if(subdet == static_cast(PixelSubdetector::PixelBarrel)){ + LorentzAngleBarrel_->Fill(it->second); + //std::cout << " bpix detid " << it->first << " \t" << " Lorentz angle " << it->second << std::endl; + //edm::LogInfo("SiPixelLorentzAngleReader") << " bpix detid " << it->first << " \t" << " Lorentz angle " << it->second; + + PXBDetId pdetId = PXBDetId(detid); + //unsigned int detTypeP=pdetId.det(); + //unsigned int subidP=pdetId.subdetId(); + // Barell layer = 1,2,3 + int layerC=pdetId.layer(); + // Barrel ladder id 1-20,32,44. + int ladderC=pdetId.ladder(); + // Barrel Z-index=1,8 + int zindex=pdetId.module(); + + if(printdebug_) { + //std::cout<<"BPix - layer "<second << std::endl; + edm::LogInfo("SiPixelLorentzAngleReader") <<"BPix - layer "<second; + } else { + + if(ladderC==1) { // print once per ring + + if(it->second != la_old) { + //std::cout<<"BPix - layer "<second << std::endl; + edm::LogInfo("SiPixelLorentzAngleReader") <<"BPix - layer "<second; + } // else {std::cout<<"same"<second; + } + } + + }else if(subdet == static_cast(PixelSubdetector::PixelEndcap)){ + LorentzAngleForward_->Fill(it->second); + + PXFDetId pdetId = PXFDetId(detid); + int disk=pdetId.disk(); //1,2,3 + int blade=pdetId.blade(); //1-24 + int moduleF=pdetId.module(); // + int side=pdetId.side(); //size=1 for -z, 2 for +z + int panel=pdetId.panel(); //panel=1 + + if(blade==1 && moduleF==1 && side==1 && panel==1) { // print once per disk + //std::cout<<"FPix - disk "<second << std::endl; + edm::LogInfo("SiPixelLorentzAngleReader") <<"FPix - disk "<second; + } + + } + } } diff --git a/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader_cfg.py b/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader_cfg.py index bd8fde01c6e2a..f239e4e0e853a 100644 --- a/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader_cfg.py +++ b/CalibTracker/SiPixelLorentzAngle/test/SiPixelLorentzAngleReader_cfg.py @@ -1,14 +1,21 @@ import FWCore.ParameterSet.Config as cms process = cms.Process("Test") + process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1) ) +# firsRun defines the IOV used process.source = cms.Source("EmptySource", - lastRun = cms.untracked.uint32(1), - timetype = cms.string('runnumber'), +# firstRun = cms.untracked.uint32(180000), +# firstRun = cms.untracked.uint32(190000), +# firstRun = cms.untracked.uint32(196000), +# firstRun = cms.untracked.uint32(198000), +# firstRun = cms.untracked.uint32(202000), +# firstRun = cms.untracked.uint32(204000), +# firstRun = cms.untracked.uint32(205000), +# firstRun = cms.untracked.uint32(208686), firstRun = cms.untracked.uint32(1), - interval = cms.uint32(1) ) @@ -16,7 +23,6 @@ fileName = cms.string("histo.root") ) - process.MessageLogger = cms.Service("MessageLogger", cout = cms.untracked.PSet( threshold = cms.untracked.string('WARNING') @@ -24,31 +30,36 @@ destinations = cms.untracked.vstring('cout') ) -process.Timing = cms.Service("Timing") +#process.Timing = cms.Service("Timing") process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") -process.load("Configuration.StandardSequences.GeometryIdeal_cff") +process.GlobalTag.globaltag = 'MC_70_V1::All' +#process.GlobalTag.globaltag = 'FT_R_53_V18::All' +#process.GlobalTag.globaltag = 'GR_P_V43::All' + +process.load("Configuration.Geometry.GeometryIdeal_cff") -process.QualityReader = cms.ESSource("PoolDBESSource", -# BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), - DBParameters = cms.PSet( - messageLevel = cms.untracked.int32(0), - authenticationPath = cms.untracked.string('') - ), - toGet = cms.VPSet( - cms.PSet( - record = cms.string("SiPixelLorentzAngleRcd"), - tag = cms.string("trivial_LorentzAngle") - ), - cms.PSet( - record = cms.string("SiPixelLorentzAngleSimRcd"), - tag = cms.string("trivial_LorentzAngle_Sim") - ) - ), - connect = cms.string('sqlite_file:SiPixelLorentzAngle.db') -) -process.es_prefer_QualityReader = cms.ESPrefer("PoolDBESSource","QualityReader") +# read from sqlite file +## process.QualityReader = cms.ESSource("PoolDBESSource", +## # BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), +## DBParameters = cms.PSet( +## messageLevel = cms.untracked.int32(0), +## authenticationPath = cms.untracked.string('') +## ), +## toGet = cms.VPSet( +## cms.PSet( +## record = cms.string("SiPixelLorentzAngleRcd"), +## tag = cms.string("trivial_LorentzAngle") +## ), +## cms.PSet( +## record = cms.string("SiPixelLorentzAngleSimRcd"), +## tag = cms.string("trivial_LorentzAngle_Sim") +## ) +## ), +## connect = cms.string('sqlite_file:SiPixelLorentzAngle.db') +## ) +## process.es_prefer_QualityReader = cms.ESPrefer("PoolDBESSource","QualityReader") process.LorentzAngleReader = cms.EDAnalyzer("SiPixelLorentzAngleReader", printDebug = cms.untracked.bool(False), @@ -63,4 +74,5 @@ process.p = cms.Path(process.LorentzAngleReader*process.LorentzAngleSimReader) +# process.p = cms.Path(process.LorentzAngleReader) diff --git a/EventFilter/SiPixelRawToDigi/test/runRawToClus_cfg.py b/EventFilter/SiPixelRawToDigi/test/runRawToClus_cfg.py new file mode 100644 index 0000000000000..7aa3183bd64ad --- /dev/null +++ b/EventFilter/SiPixelRawToDigi/test/runRawToClus_cfg.py @@ -0,0 +1,108 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("MyRawToClus") + +process.load("FWCore.MessageLogger.MessageLogger_cfi") +#process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +#process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") + +# for strips +#process.load("CalibTracker.SiStripESProducers.SiStripGainSimESProducer_cfi") + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000)) + +process.source = cms.Source("PoolSource", +# fileNames = cms.untracked.vstring('file:rawdata.root') +fileNames = cms.untracked.vstring( +# "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/205/217/2EF61B7D-F216-E211-98C3-001D09F28D54.root", + "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/208/686/A88F66A0-393F-E211-9287-002481E0D524.root", + ) +) + +#process.load("Geometry.TrackerSimData.trackerSimGeometryXML_cfi") +#process.load("Geometry.TrackerGeometryBuilder.trackerGeometry_cfi") +#process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi") +#process.load("Configuration.StandardSequences.MagneticField_cff") + +# Cabling +# include "CalibTracker/Configuration/data/Tracker_FakeConditions.cff" +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#process.GlobalTag.connect = "frontier://FrontierProd/CMS_COND_21X_GLOBALTAG" +#process.GlobalTag.globaltag = "CRAFT_V3P::All" +#process.es_prefer_GlobalTag = cms.ESPrefer('PoolDBESSource','GlobalTag') + +#process.load("CalibTracker.Configuration.SiPixel_FakeConditions_cff") +#process.load("CalibTracker.Configuration.SiPixelCabling.SiPixelCabling_SQLite_cff") +#process.siPixelCabling.connect = 'sqlite_file:cabling.db' +#process.siPixelCabling.toGet = cms.VPSet(cms.PSet( +# record = cms.string('SiPixelFedCablingMapRcd'), +# tag = cms.string('SiPixelFedCablingMap_v14') +#)) + + +# Choose the global tag here: +#process.GlobalTag.globaltag = "GR_P_V40::All" +process.GlobalTag.globaltag = "GR_R_62_V1::All" + +# process.load("EventFilter.SiPixelRawToDigi.SiPixelRawToDigi_cfi") +process.load('Configuration.StandardSequences.RawToDigi_cff') + +# needed for pixel RecHits (TkPixelCPERecord) +process.load("Configuration.StandardSequences.Reconstruction_cff") + +# clusterizer +process.load("RecoLocalTracker.Configuration.RecoLocalTracker_cff") + +# for Raw2digi for data +process.siPixelDigis.InputLabel = 'rawDataCollector' +process.siStripDigis.ProductLabel = 'rawDataCollector' + + +# for digi to clu +#process.siPixelClusters.src = 'siPixelDigis' + +process.MessageLogger = cms.Service("MessageLogger", + debugModules = cms.untracked.vstring('SiPixelClusterizer'), + destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), + cout = cms.untracked.PSet( +# threshold = cms.untracked.string('INFO') +# threshold = cms.untracked.string('ERROR') + threshold = cms.untracked.string('WARNING') + ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +) + +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +process.out = cms.OutputModule("PoolOutputModule", +# fileName = cms.untracked.string('file:digis.root'), + fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/data/clus/clus_1k.root'), + + #outputCommands = cms.untracked.vstring("drop *","keep *_*_*_MyRawToClus") # 13.1MB per 10 events + splitLevel = cms.untracked.int32(0), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + outputCommands = process.RECOEventContent.outputCommands, # 4.9MB per 10 events + dataset = cms.untracked.PSet( + filterName = cms.untracked.string(''), + dataTier = cms.untracked.string('RECO') + ) + +) + +#process.p = cms.Path(process.siPixelDigis) +#process.p = cms.Path(process.siPixelDigis*process.siPixelClusters) +#process.p = cms.Path(process.siPixelDigis*process.pixeltrackerlocalreco) + +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis) +# crash on strip clusters, calibration records missing? works with the 620 tag +process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco) + +process.ep = cms.EndPath(process.out) diff --git a/EventFilter/SiPixelRawToDigi/test/runRawToDigi_cfg.py b/EventFilter/SiPixelRawToDigi/test/runRawToDigi_cfg.py index 4c8d23045f10d..699c04346b05e 100644 --- a/EventFilter/SiPixelRawToDigi/test/runRawToDigi_cfg.py +++ b/EventFilter/SiPixelRawToDigi/test/runRawToDigi_cfg.py @@ -2,21 +2,25 @@ process = cms.Process("MyRawToDigi") +process.load("FWCore.MessageLogger.MessageLogger_cfi") +#process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +#process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") + process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000)) process.source = cms.Source("PoolSource", # fileNames = cms.untracked.vstring('file:rawdata.root') -# fileNames = cms.untracked.vstring('file:/scratch/dkotlins/COSMIC/RAW/002C2E77-8CD5-DE11-9533-000423D944FC.root') fileNames = cms.untracked.vstring( - "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/205/217/2EF61B7D-F216-E211-98C3-001D09F28D54.root", +# "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/205/217/2EF61B7D-F216-E211-98C3-001D09F28D54.root", + "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/208/686/A88F66A0-393F-E211-9287-002481E0D524.root", ) ) -process.load("Geometry.TrackerSimData.trackerSimGeometryXML_cfi") -process.load("Geometry.TrackerGeometryBuilder.trackerGeometry_cfi") -process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi") -process.load("Configuration.StandardSequences.MagneticField_cff") - # Cabling # include "CalibTracker/Configuration/data/Tracker_FakeConditions.cff" #process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") @@ -33,7 +37,7 @@ #)) -process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: +# Choose the global tag here: process.GlobalTag.globaltag = "GR_P_V40::All" @@ -54,7 +58,8 @@ ) process.out = cms.OutputModule("PoolOutputModule", - fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/digis.root'), +# fileName = cms.untracked.string('file:digis.root'), + fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/data/digis/digis_1k.root'), outputCommands = cms.untracked.vstring("drop *","keep *_siPixelDigis_*_*") ) diff --git a/EventFilter/SiPixelRawToDigi/test/runRawToTracks_cfg.py b/EventFilter/SiPixelRawToDigi/test/runRawToTracks_cfg.py new file mode 100644 index 0000000000000..537a0b494c85d --- /dev/null +++ b/EventFilter/SiPixelRawToDigi/test/runRawToTracks_cfg.py @@ -0,0 +1,133 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("MyRawToTracks") + +process.load("FWCore.MessageLogger.MessageLogger_cfi") +#process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +#process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") + +# for strips +#process.load("CalibTracker.SiStripESProducers.SiStripGainSimESProducer_cfi") + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10)) + +process.source = cms.Source("PoolSource", +# fileNames = cms.untracked.vstring('file:rawdata.root') +fileNames = cms.untracked.vstring( +# "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/205/217/2EF61B7D-F216-E211-98C3-001D09F28D54.root", + "rfio:/castor/cern.ch/cms/store/data/Run2012D/MinimumBias/RAW/v1/000/208/686/A88F66A0-393F-E211-9287-002481E0D524.root", + ) +# skipEvents = cms.untracked.uint32(5000) +) + +#process.load("Geometry.TrackerSimData.trackerSimGeometryXML_cfi") +#process.load("Geometry.TrackerGeometryBuilder.trackerGeometry_cfi") +#process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi") +#process.load("Configuration.StandardSequences.MagneticField_cff") + +# Cabling +# include "CalibTracker/Configuration/data/Tracker_FakeConditions.cff" +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#process.GlobalTag.connect = "frontier://FrontierProd/CMS_COND_21X_GLOBALTAG" +#process.GlobalTag.globaltag = "CRAFT_V3P::All" +#process.es_prefer_GlobalTag = cms.ESPrefer('PoolDBESSource','GlobalTag') + +#process.load("CalibTracker.Configuration.SiPixel_FakeConditions_cff") +#process.load("CalibTracker.Configuration.SiPixelCabling.SiPixelCabling_SQLite_cff") +#process.siPixelCabling.connect = 'sqlite_file:cabling.db' +#process.siPixelCabling.toGet = cms.VPSet(cms.PSet( +# record = cms.string('SiPixelFedCablingMapRcd'), +# tag = cms.string('SiPixelFedCablingMap_v14') +#)) + + +# Choose the global tag here: +#process.GlobalTag.globaltag = "GR_P_V40::All" +process.GlobalTag.globaltag = "GR_R_62_V1::All" + +# process.load("EventFilter.SiPixelRawToDigi.SiPixelRawToDigi_cfi") +process.load('Configuration.StandardSequences.RawToDigi_cff') + +# needed for pixel RecHits (TkPixelCPERecord) +process.load("Configuration.StandardSequences.Reconstruction_cff") + +# clusterizer +process.load("RecoLocalTracker.Configuration.RecoLocalTracker_cff") + +# for Raw2digi for data +process.siPixelDigis.InputLabel = 'rawDataCollector' +process.siStripDigis.ProductLabel = 'rawDataCollector' + + +# for digi to clu +#process.siPixelClusters.src = 'siPixelDigis' + +process.MessageLogger = cms.Service("MessageLogger", + debugModules = cms.untracked.vstring('SiPixelClusterizer'), + destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), + cout = cms.untracked.PSet( +# threshold = cms.untracked.string('INFO') +# threshold = cms.untracked.string('ERROR') + threshold = cms.untracked.string('WARNING') + ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +) + +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +process.out = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('file:tracks.root'), +# fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/data/tracks/tracks_1_0.root'), + + #outputCommands = cms.untracked.vstring("drop *","keep *_*_*_MyRawToClus") # 13.1MB per 10 events + splitLevel = cms.untracked.int32(0), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + outputCommands = process.RECOEventContent.outputCommands, # 4.9MB per 10 events + dataset = cms.untracked.PSet( + filterName = cms.untracked.string(''), + dataTier = cms.untracked.string('RECO') + ) + +) + +# copy the sequence below from +# RecoTracker/IterativeTracking/python/iterativeTk_cff.py +process.myTracking = cms.Sequence(process.InitialStep* + process.LowPtTripletStep* + process.PixelPairStep* + process.DetachedTripletStep* + process.MixedTripletStep* + process.PixelLessStep* + process.TobTecStep* + process.earlyGeneralTracks* + #process.muonSeededStep* + process.preDuplicateMergingGeneralTracks* + process.generalTracksSequence* + process.ConvStep* + process.conversionStepTracks + ) + +#process.p = cms.Path(process.siPixelDigis) +#process.p = cms.Path(process.siPixelDigis*process.siPixelClusters) +#process.p = cms.Path(process.siPixelDigis*process.pixeltrackerlocalreco) + +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis) +# crash on strip clusters, calibration records missing? works with the 620 tag +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco) + +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot) +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.recopixelvertexing) +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.recopixelvertexing*process.MeasurementTrackerEvent) +# trackingGlobalReco, ckftracks, iterTracking - +#process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.recopixelvertexing*process.MeasurementTrackerEvent*process.myTracking) +process.p1 = cms.Path(process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.recopixelvertexing*process.MeasurementTrackerEvent*process.myTracking*process.vertexreco) + +process.ep = cms.EndPath(process.out) diff --git a/RecoLocalTracker/SiPixelClusterizer/test/BuildFile.xml b/RecoLocalTracker/SiPixelClusterizer/test/BuildFile.xml index f073e6abba4af..b8a4daa03783f 100644 --- a/RecoLocalTracker/SiPixelClusterizer/test/BuildFile.xml +++ b/RecoLocalTracker/SiPixelClusterizer/test/BuildFile.xml @@ -28,9 +28,9 @@ -# -# -# -# -# -# + + + + + + diff --git a/RecoLocalTracker/SiPixelClusterizer/test/ReadPixClusters.cc b/RecoLocalTracker/SiPixelClusterizer/test/ReadPixClusters.cc index 2b05c27dac3ec..f4ebaea38fa06 100644 --- a/RecoLocalTracker/SiPixelClusterizer/test/ReadPixClusters.cc +++ b/RecoLocalTracker/SiPixelClusterizer/test/ReadPixClusters.cc @@ -113,9 +113,9 @@ class ReadPixClusters : public edm::EDAnalyzer { TH1F *hpixcharge1,*hpixcharge2, *hpixcharge3, *hpixcharge4, *hpixcharge5; TH1F *hcols1,*hcols2,*hcols3,*hrows1,*hrows2,*hrows3; TH1F *hpcols1,*hpcols2,*hpcols3,*hprows1,*hprows2,*hprows3; - TH1F *hsize1,*hsize2,*hsize3, - *hsizex1,*hsizex2,*hsizex3, - *hsizey1,*hsizey2,*hsizey3; + TH1F *hsize1,*hsize2,*hsize3,*hsize4,*hsize5, + *hsizex1,*hsizex2,*hsizex3,*hsizex4,*hsizex5, + *hsizey1,*hsizey2,*hsizey3,*hsizey4,*hsizey5; TH1F *hclusPerDet1,*hclusPerDet2,*hclusPerDet3; TH1F *hpixPerDet1,*hpixPerDet2,*hpixPerDet3; @@ -261,17 +261,17 @@ void ReadPixClusters::beginJob() { hdetsPerLay2 = fs->make( "hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5); - sizeH=1000; + sizeH=120; lowH = 0.; - highH = 100.0; // charge limit in kelec + highH = 121.0; // charge limit in kelec hcharge1 = fs->make( "hcharge1", "Clu charge l1", sizeH, 0.,highH); //in ke hcharge2 = fs->make( "hcharge2", "Clu charge l2", sizeH, 0.,highH); hcharge3 = fs->make( "hcharge3", "Clu charge l3", sizeH, 0.,highH); hcharge4 = fs->make( "hcharge4", "Clu charge d1", sizeH, 0.,highH); hcharge5 = fs->make( "hcharge5", "Clu charge d2", sizeH, 0.,highH); - sizeH=600; - highH = 60.0; // charge limit in kelec + sizeH=90; + highH = 61.0; // charge limit in kelec hpixcharge1 = fs->make( "hpixcharge1", "Pix charge l1",sizeH, 0.,highH);//in ke hpixcharge2 = fs->make( "hpixcharge2", "Pix charge l2",sizeH, 0.,highH); hpixcharge3 = fs->make( "hpixcharge3", "Pix charge l3",sizeH, 0.,highH); @@ -300,6 +300,8 @@ void ReadPixClusters::beginJob() { hsize1 = fs->make( "hsize1", "layer 1 clu size",sizeH,-0.5,highH); hsize2 = fs->make( "hsize2", "layer 2 clu size",sizeH,-0.5,highH); hsize3 = fs->make( "hsize3", "layer 3 clu size",sizeH,-0.5,highH); + hsize4 = fs->make( "hsize4", "disk 1 clu size",sizeH,-0.5,highH); + hsize5 = fs->make( "hsize5", "disk 2 clu size",sizeH,-0.5,highH); hsizex1 = fs->make( "hsizex1", "lay1 clu size in x", 10,-0.5,9.5); @@ -307,12 +309,21 @@ void ReadPixClusters::beginJob() { 10,-0.5,9.5); hsizex3 = fs->make( "hsizex3", "lay3 clu size in x", 10,-0.5,9.5); + hsizex4 = fs->make( "hsizex4", "d1 clu size in x", + 10,-0.5,9.5); + hsizex5 = fs->make( "hsizex5", "d2 clu size in x", + 10,-0.5,9.5); + hsizey1 = fs->make( "hsizey1", "lay1 clu size in y", 20,-0.5,19.5); hsizey2 = fs->make( "hsizey2", "lay2 clu size in y", 20,-0.5,19.5); hsizey3 = fs->make( "hsizey3", "lay3 clu size in y", 20,-0.5,19.5); + hsizey4 = fs->make( "hsizey4", "d1 clu size in y", + 20,-0.5,19.5); + hsizey5 = fs->make( "hsizey5", "d2 clu size in y", + 20,-0.5,19.5); hevent = fs->make("hevent","event",1000,0,10000000.); horbit = fs->make("horbit","orbit",100, 0,100000000.); @@ -408,7 +419,7 @@ void ReadPixClusters::analyze(const edm::Event& e, hdets->Fill(float(numOf)); // number of modules with pix // Select events with pixels - if(numOf<1) return; // skip events with pixel dets + //if(numOf<1) return; // skip events with pixel dets //if(numOf<4) return; // skip events with few pixel dets hevent->Fill(float(event)); @@ -846,6 +857,9 @@ void ReadPixClusters::analyze(const edm::Event& e, hcharge4->Fill(ch); aveCharge4 += ch; + hsize4->Fill(float(size)); + hsizex4->Fill(float(sizeX)); + hsizey4->Fill(float(sizeY)); } else if(disk==2) { // disk2 -+z @@ -855,6 +869,9 @@ void ReadPixClusters::analyze(const edm::Event& e, hcharge5->Fill(ch); aveCharge5 += ch; + hsize5->Fill(float(size)); + hsizex5->Fill(float(sizeX)); + hsizey5->Fill(float(sizeY)); } else cout<<" unknown disk "< +#include +#include + +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +//#include "DataFormats/Common/interface/Handle.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "DataFormats/Common/interface/EDProduct.h" + +//#include "DataFormats/SiPixelCluster/interface/SiPixelClusterCollection.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/DetId/interface/DetId.h" + +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" + +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" +//#include "Geometry/CommonTopologies/interface/PixelTopology.h" + +// For L1 +#include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h" + +// For HLT +#include "DataFormats/HLTReco/interface/TriggerEvent.h" +#include "DataFormats/HLTReco/interface/TriggerTypeDefs.h" +#include "DataFormats/Common/interface/TriggerResults.h" +#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h" +#include "FWCore/Common/interface/TriggerNames.h" + + +// For tracks +#include "DataFormats/TrackReco/interface/Track.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" + +//#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" + +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" +//#include "TrackingTools/PatternTools/interface/TrajectoryFitter.h" + +#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" + +//#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h" +//#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h" + +//#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h" + + +//#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h" +//#include "TrackingTools/PatternTools/interface/Trajectory.h" +//#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +//#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +//#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include + +// For luminisoty +#include "DataFormats/Luminosity/interface/LumiSummary.h" +#include "DataFormats/Common/interface/ConditionsInEdm.h" + +// To use root histos +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +// For ROOT +#include +#include +#include +#include +#include +#include +#include +#include + + +#define HISTOS +//#define L1 +//#define HLT + +using namespace std; + +class TestPixTracks : public edm::EDAnalyzer { + public: + + explicit TestPixTracks(const edm::ParameterSet& conf); + virtual ~TestPixTracks(); + virtual void analyze(const edm::Event& e, const edm::EventSetup& c); + virtual void beginRun(const edm::EventSetup& iSetup); + virtual void beginJob(); + virtual void endJob(); + + private: + edm::ParameterSet conf_; + edm::InputTag src_; + //const static bool PRINT = false; + bool PRINT; + float countTracks, countGoodTracks, countTracksInPix, countPVs, countEvents, countLumi; + + //TFile* hFile; + //TH1D *hdetunit; + //TH1D *hpixid,*hpixsubid, + //*hlayerid, + //*hladder1id,*hladder2id,*hladder3id, + //*hz1id,*hz2id,*hz3id; + + TH1D *hcharge1,*hcharge2, *hcharge3, *hcharge; + TH1D *hpixcharge1,*hpixcharge2, *hpixcharge3, *hpixcharge; + TH1D *hcols1,*hcols2,*hcols3,*hrows1,*hrows2,*hrows3; + TH1D *hsize1,*hsize2,*hsize3, + *hsizex1,*hsizex2,*hsizex3, + *hsizey1,*hsizey2,*hsizey3; + + TH1D *hclusPerTrk1,*hclusPerTrk2,*hclusPerTrk3; + TH1D *hclusPerLay1,*hclusPerLay2,*hclusPerLay3; + TH1D *hpixPerLay1,*hpixPerLay2,*hpixPerLay3; + //TH1D *hdetsPerLay1,*hdetsPerLay2,*hdetsPerLay3; + + //TH1D *hdetr, *hdetz; + // TH1D *hcolsB, *hrowsB, *hcolsF, *hrowsF; + TH2F *hDetMap1, *hDetMap2, *hDetMap3; // clusters + //TH2F *hpixDetMap1, *hpixDetMap2, *hpixDetMap3; + TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3; + + TH2F *hpvxy, *hclusMap1, *hclusMap2, *hclusMap3; + + TH1D *hpvz, *hpvr, *hNumPv, *hNumPvClean; + TH1D *hPt, *hEta, *hDz, *hD0,*hzdiff; + + //TH1D *hncharge1,*hncharge2, *hncharge3; + //TH1D *hchargeMonoPix1,*hchargeMonoPix2, *hchargeMonoPix3; + // TH1D *hnpixcharge1,*hnpixcharge2, *hnpixcharge3; + //TH1D *htest1,*htest2,*htest3,*htest4,*htest5,*htest6,*htest7,*htest8,*htest9; + TH1D *hl1a, *hl1t, *hlt1; + + TH1D *hclusBpix, *hpixBpix; + TH1D *htracks, *htracksGood, *htracksGoodInPix; + + TProfile *hclumult1, *hclumult2, *hclumult3; + TProfile *hclumultx1, *hclumultx2, *hclumultx3; + TProfile *hclumulty1, *hclumulty2, *hclumulty3; + TProfile *hcluchar1, *hcluchar2, *hcluchar3; + TProfile *hpixchar1, *hpixchar2, *hpixchar3; + + TProfile *htracksls, *hpvsls, *htrackslsn, *hpvslsn, *hintgl, *hinstl, *hbeam1, *hbeam2; + + TH1D *hlumi, *hlumi0, *hbx, *hbx0; + +}; +///////////////////////////////////////////////////////////////// +// Contructor, +TestPixTracks::TestPixTracks(edm::ParameterSet const& conf) +// : conf_(conf), src_(conf.getParameter( "src" )) { } + : conf_(conf) { + + PRINT = conf.getUntrackedParameter("Verbosity",false); + src_ = conf.getParameter( "src" ); + //if(PRINT) cout<<" Construct "< fs; + + // put here whatever you want to do at the beginning of the job + //hFile = new TFile ( "histo.root", "RECREATE" ); + + //hladder1id = fs->make( "hladder1id", "Ladder L1 id", 50, 0., 50.); + //hladder2id = fs->make( "hladder2id", "Ladder L2 id", 50, 0., 50.); + //hladder3id = fs->make( "hladder3id", "Ladder L3 id", 50, 0., 50.); + //hz1id = fs->make( "hz1id", "Z-index id L1", 10, 0., 10.); + //hz2id = fs->make( "hz2id", "Z-index id L2", 10, 0., 10.); + //hz3id = fs->make( "hz3id", "Z-index id L3", 10, 0., 10.); + + int sizeH=20; + float lowH = -0.5; + float highH = 19.5; + + hclusPerTrk1 = fs->make( "hclusPerTrk1", "Clus per track l1", + sizeH, lowH, highH); + hclusPerTrk2 = fs->make( "hclusPerTrk2", "Clus per track l2", + sizeH, lowH, highH); + hclusPerTrk3 = fs->make( "hclusPerTrk3", "Clus per track l3", + sizeH, lowH, highH); + + sizeH=2000; + highH = 1999.5; + hclusPerLay1 = fs->make( "hclusPerLay1", "Clus per layer l1", + sizeH, lowH, highH); + hclusPerLay2 = fs->make( "hclusPerLay2", "Clus per layer l2", + sizeH, lowH, highH); + hclusPerLay3 = fs->make( "hclusPerLay3", "Clus per layer l3", + sizeH, lowH, highH); + + highH = 9999.5; + hpixPerLay1 = fs->make( "hpixPerLay1", "Pix per layer l1", + sizeH, lowH, highH); + hpixPerLay2 = fs->make( "hpixPerLay2", "Pix per layer l2", + sizeH, lowH, highH); + hpixPerLay3 = fs->make( "hpixPerLay3", "Pix per layer l3", + sizeH, lowH, highH); + + //hdetsPerLay1 = fs->make( "hdetsPerLay1", "Full dets per layer l1", + // 161, -0.5, 160.5); + //hdetsPerLay3 = fs->make( "hdetsPerLay3", "Full dets per layer l3", + // 353, -0.5, 352.5); + //hdetsPerLay2 = fs->make( "hdetsPerLay2", "Full dets per layer l2", + // 257, -0.5, 256.5); + + hcharge1 = fs->make( "hcharge1", "Clu charge l1", 400, 0.,200.); //in ke + hcharge2 = fs->make( "hcharge2", "Clu charge l2", 400, 0.,200.); + hcharge3 = fs->make( "hcharge3", "Clu charge l3", 400, 0.,200.); + + //hchargeMonoPix1 = fs->make( "hchargeMonoPix1", "Clu charge l1 MonPix", 200, 0.,100.); //in ke + //hchargeMonoPix2 = fs->make( "hchargeMonoPix2", "Clu charge l2 MonPix", 200, 0.,100.); + //hchargeMonoPix3 = fs->make( "hchargeMonoPix3", "Clu charge l3 MonPix", 200, 0.,100.); + + //hncharge1 = fs->make( "hncharge1", "Noise charge l1", 200, 0.,100.); //in ke + //hncharge2 = fs->make( "hncharge2", "Noise charge l2", 200, 0.,100.); + //hncharge3 = fs->make( "hncharge3", "Noise charge l3", 200, 0.,100.); + + //hpixcharge1 = fs->make( "hpixcharge1", "Pix charge l1", 100, 0.,50.); //in ke + //hpixcharge2 = fs->make( "hpixcharge2", "Pix charge l2", 100, 0.,50.); + //hpixcharge3 = fs->make( "hpixcharge3", "Pix charge l3", 100, 0.,50.); + + //hnpixcharge1 = fs->make( "hnpixcharge1", "Noise pix charge l1", 100, 0.,50.); //in ke + //hnpixcharge2 = fs->make( "hnpixcharge2", "Noise pix charge l2", 100, 0.,50.); + //hnpixcharge3 = fs->make( "hnpixcharge3", "Noise pix charge l3", 100, 0.,50.); + + //hpixcharge = fs->make( "hpixcharge", "Clu charge", 100, 0.,50.); + //hcharge = fs->make( "hcharge", "Pix charge", 100, 0.,50.); + + hcols1 = fs->make( "hcols1", "Layer 1 cols", 500,-0.5,499.5); + hcols2 = fs->make( "hcols2", "Layer 2 cols", 500,-0.5,499.5); + hcols3 = fs->make( "hcols3", "Layer 3 cols", 500,-0.5,499.5); + + hrows1 = fs->make( "hrows1", "Layer 1 rows", 200,-0.5,199.5); + hrows2 = fs->make( "hrows2", "Layer 2 rows", 200,-0.5,199.5); + hrows3 = fs->make( "hrows3", "layer 3 rows", 200,-0.5,199.5); + + hsize1 = fs->make( "hsize1", "layer 1 clu size",100,-0.5,99.5); + hsize2 = fs->make( "hsize2", "layer 2 clu size",100,-0.5,99.5); + hsize3 = fs->make( "hsize3", "layer 3 clu size",100,-0.5,99.5); + hsizex1 = fs->make( "hsizex1", "lay1 clu size in x", + 20,-0.5,19.5); + hsizex2 = fs->make( "hsizex2", "lay2 clu size in x", + 20,-0.5,19.5); + hsizex3 = fs->make( "hsizex3", "lay3 clu size in x", + 20,-0.5,19.5); + hsizey1 = fs->make( "hsizey1", "lay1 clu size in y", + 30,-0.5,29.5); + hsizey2 = fs->make( "hsizey2", "lay2 clu size in y", + 30,-0.5,29.5); + hsizey3 = fs->make( "hsizey3", "lay3 clu size in y", + 30,-0.5,29.5); + + + hDetMap1 = fs->make( "hDetMap1", "layer 1 clus map", + 9,0.,9.,23,0.,23.); + hDetMap2 = fs->make( "hDetMap2", "layer 2 clus map", + 9,0.,9.,33,0.,33.); + hDetMap3 = fs->make( "hDetMap3", "layer 3 clus map", + 9,0.,9.,45,0.,45.); + //hpixDetMap1 = fs->make( "hpixDetMap1", "pix det layer 1", + // 416,0.,416.,160,0.,160.); + //hpixDetMap2 = fs->make( "hpixDetMap2", "pix det layer 2", + // 416,0.,416.,160,0.,160.); + //hpixDetMap3 = fs->make( "hpixDetMap3", "pix det layer 3", + // 416,0.,416.,160,0.,160.); + hcluDetMap1 = fs->make( "hcluDetMap1", "clu det layer 1", + 416,0.,416.,160,0.,160.); + hcluDetMap2 = fs->make( "hcluDetMap2", "clu det layer 1", + 416,0.,416.,160,0.,160.); + hcluDetMap3 = fs->make( "hcluDetMap3", "clu det layer 1", + 416,0.,416.,160,0.,160.); + + htracksGoodInPix = fs->make( "htracksGoodInPix", "count good tracks in pix",2000,-0.5,1999.5); + htracksGood = fs->make( "htracksGood", "count good tracks",2000,-0.5,1999.5); + htracks = fs->make( "htracks", "count tracks",2000,-0.5,1999.5); + hclusBpix = fs->make( "hclusBpix", "count clus in bpix",200,-0.5,1999.5); + hpixBpix = fs->make( "hpixBpix", "count pixels",200,-0.5,1999.5); + + hpvxy = fs->make( "hpvxy", "pv xy",100,-1.,1.,100,-1.,1.); + hpvz = fs->make( "hpvz", "pv z",1000,-50.,50.); + hpvr = fs->make( "hpvr", "pv r",100,0.,1.); + hNumPv = fs->make( "hNumPv", "num of pv",100,0.,100.); + hNumPvClean = fs->make( "hNumPvClean", "num of pv clean",100,0.,100.); + + hPt = fs->make( "hPt", "pt",100,0.,100.); + hEta = fs->make( "hEta", "eta",50,-2.5,2.5); + hD0 = fs->make( "hD0", "d0",500,0.,5.); + hDz = fs->make( "hDz", "pt",250,-25.,25.); + hzdiff = fs->make( "hzdiff", "PVz-Trackz",200,-10.,10.); + + hl1a = fs->make("hl1a", "l1a", 128,-0.5,127.5); + hl1t = fs->make("hl1t", "l1t", 128,-0.5,127.5); + hlt1 = fs->make("hlt1","hlt1",256,-0.5,255.5); + + hclumult1 = fs->make("hclumult1","cluster size layer 1",60,-3.,3.,0.0,100.); + hclumult2 = fs->make("hclumult2","cluster size layer 2",60,-3.,3.,0.0,100.); + hclumult3 = fs->make("hclumult3","cluster size layer 3",60,-3.,3.,0.0,100.); + + hclumultx1 = fs->make("hclumultx1","cluster x-size layer 1",60,-3.,3.,0.0,100.); + hclumultx2 = fs->make("hclumultx2","cluster x-size layer 2",60,-3.,3.,0.0,100.); + hclumultx3 = fs->make("hclumultx3","cluster x-size layer 3",60,-3.,3.,0.0,100.); + + hclumulty1 = fs->make("hclumulty1","cluster y-size layer 1",60,-3.,3.,0.0,100.); + hclumulty2 = fs->make("hclumulty2","cluster y-size layer 2",60,-3.,3.,0.0,100.); + hclumulty3 = fs->make("hclumulty3","cluster y-size layer 3",60,-3.,3.,0.0,100.); + + hcluchar1 = fs->make("hcluchar1","cluster char layer 1",60,-3.,3.,0.0,1000.); + hcluchar2 = fs->make("hcluchar2","cluster char layer 2",60,-3.,3.,0.0,1000.); + hcluchar3 = fs->make("hcluchar3","cluster char layer 3",60,-3.,3.,0.0,1000.); + + hpixchar1 = fs->make("hpixchar1","pix char layer 1",60,-3.,3.,0.0,1000.); + hpixchar2 = fs->make("hpixchar2","pix char layer 2",60,-3.,3.,0.0,1000.); + hpixchar3 = fs->make("hpixchar3","pix char layer 3",60,-3.,3.,0.0,1000.); + + hintgl = fs->make("hintgl", "inst lumi vs ls ",1000,0.,3000.,0.0,10000.); + hinstl = fs->make("hinstl", "intg lumi vs ls ",1000,0.,3000.,0.0,100.); + hbeam1 = fs->make("hbeam1", "beam1 vs ls ",1000,0.,3000.,0.0,1000.); + hbeam2 = fs->make("hbeam2", "beam2 vs ls ",1000,0.,3000.,0.0,1000.); + + htracksls = fs->make("htracksls","tracks with pix hits vs ls",1000,0.,3000.,0.0,10000.); + hpvsls = fs->make("hpvsls","pvs vs ls",1000,0.,3000.,0.0,1000.); + htrackslsn = fs->make("htrackslsn","tracks with pix hits/lumi vs ls",1000,0.,3000.,0.0,10000.); + hpvslsn = fs->make("hpvslsn","pvs/lumi vs ls",1000,0.,3000.,0.0,1000.); + + hlumi0 = fs->make("hlumi0", "lumi", 2000,0,2000.); + hlumi = fs->make("hlumi", "lumi", 2000,0,2000.); + hbx0 = fs->make("hbx0", "bx", 4000,0,4000.); + hbx = fs->make("hbx", "bx", 4000,0,4000.); + + hclusMap1 = fs->make("hclusMap1","clus - lay1",260,-26.,26.,350,-3.5,3.5); + hclusMap2 = fs->make("hclusMap2","clus - lay2",260,-26.,26.,350,-3.5,3.5); + hclusMap3 = fs->make("hclusMap3","clus - lay3",260,-26.,26.,350,-3.5,3.5); + +#endif + +} +// ------------ method called to at the end of the job ------------ +void TestPixTracks::endJob(){ + cout << " End PixelTracksTest " << endl; + + if(countEvents>0.) { + countTracks /= countEvents; + countGoodTracks /= countEvents; + countTracksInPix /= countEvents; + countPVs /= countEvents; + countLumi /= 1000.; + cout<<" Average tracks/event "<< countTracks<<" good "<< countGoodTracks<<" in pix "<< countTracksInPix + <<" PVs "<< countPVs<<" events " << countEvents<<" lumi pb-1 "<< countLumi<<"/10, bug!"<Fill(float(bx)); + hlumi0->Fill(float(lumiBlock)); + + edm::LuminosityBlock const& iLumi = e.getLuminosityBlock(); + edm::Handle lumi; + iLumi.getByLabel("lumiProducer", lumi); + edm::Handle cond; + float intlumi = 0, instlumi=0; + int beamint1=0, beamint2=0; + iLumi.getByLabel("conditionsInEdm", cond); + // This will only work when running on RECO until (if) they fix it in the FW + // When running on RAW and reconstructing, the LumiSummary will not appear + // in the event before reaching endLuminosityBlock(). Therefore, it is not + // possible to get this info in the event + if (lumi.isValid()) { + intlumi =(lumi->intgRecLumi())/1000.; // 10^30 -> 10^33/cm2/sec -> 1/nb/sec + instlumi=(lumi->avgInsDelLumi())/1000.; + beamint1=(cond->totalIntensityBeam1)/1000; + beamint2=(cond->totalIntensityBeam2)/1000; + } else { + //std::cout << "** ERROR: Event does not get lumi info\n"; + } + + hinstl->Fill(float(lumiBlock),float(instlumi)); + hintgl->Fill(float(lumiBlock),float(intlumi)); + hbeam1->Fill(float(lumiBlock),float(beamint1)); + hbeam2->Fill(float(lumiBlock),float(beamint2)); + + + + +#ifdef L1 + // Get L1 + Handle L1GTRR; + e.getByLabel("gtDigis",L1GTRR); + + if (L1GTRR.isValid()) { + //bool l1a = L1GTRR->decision(); // global decission? + //cout<<" L1 status :"<decisionWord().size(); ++i) { + int l1flag = L1GTRR->decisionWord()[i]; + int t1flag = L1GTRR->technicalTriggerWord()[i]; + + if( l1flag>0 ) hl1a->Fill(float(i)); + if( t1flag>0 && i<64) hl1t->Fill(float(i)); + } // for loop + } // if l1a +#endif + +#ifdef HLT + + bool hlt[256]; + for(int i=0;i<256;++i) hlt[i]=false; + + edm::TriggerNames TrigNames; + edm::Handle HLTResults; + + // Extract the HLT results + e.getByLabel(edm::InputTag("TriggerResults","","HLT"),HLTResults); + if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) { + + //TrigNames.init(*HLTResults); + const edm::TriggerNames & TrigNames = e.triggerNames(*HLTResults); + + //cout<wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) && + (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) && + (HLTResults->error( TrigNames.triggerIndex(TrigNames.triggerName(i))) == false) ) { + + hlt[i]=true; + hlt1->Fill(float(i)); + + } // if hlt + + } // loop + } // if valid +#endif + + + // -- Does this belong into beginJob()? + //ESHandle TG; + //iSetup.get().get(TG); + //const TrackerGeometry* theTrackerGeometry = TG.product(); + //const TrackerGeometry& theTracker(*theTrackerGeometry); + // Get event setup + edm::ESHandle geom; + es.get().get( geom ); + const TrackerGeometry& theTracker(*geom); + + + // -- Primary vertices + // ---------------------------------------------------------------------- + edm::Handle vertices; + e.getByLabel("offlinePrimaryVertices", vertices); + + if(PRINT) cout<<" PV list "<size()< pvzVector; + for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) { + + if( (iv->isFake()) == 1 ) continue; + pvNotFake++; + float pvx = iv->x(); + float pvy = iv->y(); + float pvz = iv->z(); + int numTracksPerPV = iv->tracksSize(); + //int numTracksPerPV = iv->nTracks(); + + //float xe = iv->xError(); + //float ye = iv->yError(); + //float ze = iv->zError(); + //int chi2 = iv->chi2(); + //int dof = iv->ndof(); + + if(PRINT) cout<<" PV "<Fill(pvz); + if(pvz>-22. && pvz<22.) { + float pvr = sqrt(pvx*pvx + pvy*pvy); + hpvxy->Fill(pvx,pvy); + hpvr->Fill(pvr); + if(pvr<0.3) { + pvsTrue++; + pvzVector.push_back(pvz); + //if(PRINT) cout<<"PV "<Fill(float(pvNotFake)); + hNumPvClean->Fill(float(pvsTrue)); + + if(PRINT) cout<<" Not fake PVs = "< recTracks; + // e.getByLabel("generalTracks", recTracks); + // e.getByLabel("ctfWithMaterialTracksP5", recTracks); + // e.getByLabel("splittedTracksP5", recTracks); + //e.getByLabel("cosmictrackfinderP5", recTracks); + e.getByLabel(src_ , recTracks); + + + if(PRINT) cout<<" Tracks "<size()<begin(); + t!=recTracks->end(); ++t){ + + trackNumber++; + numOfClustersPerDet1=0; + numOfClustersPerDet2=0; + numOfClustersPerDet3=0; + int pixelHits=0; + + int size = t->recHitsSize(); + float pt = t->pt(); + float eta = t->eta(); + float phi = t->phi(); + float trackCharge = t->charge(); + float d0 = t->d0(); + float dz = t->dz(); + float tkvx = t->vx(); + float tkvy = t->vy(); + float tkvz = t->vz(); + + if(PRINT) cout<<"Track "<Fill(eta); + hDz->Fill(dz); + if(abs(eta)>2.8 || abs(dz)>25.) continue; // skip + + hD0->Fill(d0); + if(d0>1.0) continue; // skip + + bool goodTrack = false; + for(vector::iterator m=pvzVector.begin(); m!=pvzVector.end();++m) { + float z = *m; + float tmp = abs(z-dz); + hzdiff->Fill(tmp); + if( tmp < 1.) goodTrack=true; + } + + if(!goodTrack) continue; + countNiceTracks++; + hPt->Fill(pt); + + // Loop over rechits + for ( trackingRecHit_iterator recHit = t->recHitsBegin(); + recHit != t->recHitsEnd(); ++recHit ) { + + if ( !((*recHit)->isValid()) ) continue; + + if( (*recHit)->geographicalId().det() != DetId::Tracker ) continue; + + const DetId & hit_detId = (*recHit)->geographicalId(); + uint IntSubDetID = (hit_detId.subdetId()); + + // Select pixel rechits + if(IntSubDetID == 0 ) continue; // Select ?? + if(IntSubDetID != PixelSubdetector::PixelBarrel) continue; // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) { + + // Pixel detector + PXBDetId pdetId = PXBDetId(hit_detId); + //unsigned int detTypeP=pdetId.det(); + //unsigned int subidP=pdetId.subdetId(); + // Barell layer = 1,2,3 + int layer=pdetId.layer(); + // Barrel ladder id 1-20,32,44. + int ladder=pdetId.ladder(); + // Barrel Z-index=1,8 + int zindex=pdetId.module(); + + if(PRINT) cout<<"barrel layer/ladder/module: "< (theTracker.idToDet(hit_detId) ); + double detZ = theGeomDet->surface().position().z(); + double detR = theGeomDet->surface().position().perp(); + const PixelTopology * topol = &(theGeomDet->specificTopology()); // pixel topology + + //std::vector output = getRecHitComponents((*recHit).get()); + //std::vector TrkComparison::getRecHitComponents(const TrackingRecHit* rechit){ + + const SiPixelRecHit* hit = dynamic_cast((*recHit).get()); + //edm::Ref ,SiStripCluster> cluster = hit->cluster(); + // get the edm::Ref to the cluster + + if(hit) { + + // RecHit (recthits are transient, so not available without ttrack refit) + //double xloc = hit->localPosition().x();// 1st meas coord + //double yloc = hit->localPosition().y();// 2nd meas coord or zero + //double zloc = hit->localPosition().z();// up, always zero + //cout<<" rechit loc "<, SiPixelCluster> const& clust = hit->cluster(); + // check if the ref is not null + if (!clust.isNonnull()) continue; + + numberOfClusters++; + pixelHits++; + float charge = (clust->charge())/1000.0; // convert electrons to kilo-electrons + int size = clust->size(); + int sizeX = clust->sizeX(); + int sizeY = clust->sizeY(); + float row = clust->x(); + float col = clust->y(); + numberOfPixels += size; + + //cout<<" clus loc "<surface().toGlobal( lp ); + double gX = clustgp.x(); + double gY = clustgp.y(); + double gZ = clustgp.z(); + + //cout<<" CLU GLOBAL "<maxPixelCol(); + //int maxPixelRow = clust->maxPixelRow(); + //int minPixelCol = clust->minPixelCol(); + //int minPixelRow = clust->minPixelRow(); + //int geoId = PixGeom->geographicalId().rawId(); + // Replace with the topology methods + // edge method moved to topologi class + //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) ); + //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) ); + + // calculate alpha and beta from cluster position + //LocalTrajectoryParameters ltp = tsos.localParameters(); + //LocalVector localDir = ltp.momentum()/ltp.momentum().mag(); + //float locx = localDir.x(); + //float locy = localDir.y(); + //float locz = localDir.z(); + //float loctheta = localDir.theta(); // currently unused + //float alpha = atan2( locz, locx ); + //float beta = atan2( locz, locy ); + + if(layer==1) { + + hDetMap1->Fill(float(zindex),float(ladder)); + hcluDetMap1->Fill(col,row); + hcharge1->Fill(charge); + hcols1->Fill(col); + hrows1->Fill(row); + + hclusMap1->Fill(gZ,phi); + + if(pt>CLU_SIZE_PT_CUT) { + hsize1->Fill(float(size)); + hsizex1->Fill(float(sizeX)); + hsizey1->Fill(float(sizeY)); + + hclumult1->Fill(eta,float(size)); + hclumultx1->Fill(eta,float(sizeX)); + hclumulty1->Fill(eta,float(sizeY)); + hcluchar1->Fill(eta,float(charge)); + } + + numOfClustersPerDet1++; + numOfClustersPerLay1++; + numOfPixelsPerLay1 += size; + + } else if(layer==2) { + + hDetMap2->Fill(float(zindex),float(ladder)); + hcluDetMap2->Fill(col,row); + hcharge2->Fill(charge); + hcols2->Fill(col); + hrows2->Fill(row); + + hclusMap2->Fill(gZ,phi); + + if(pt>CLU_SIZE_PT_CUT) { + hsize2->Fill(float(size)); + hsizex2->Fill(float(sizeX)); + hsizey2->Fill(float(sizeY)); + + hclumult2->Fill(eta,float(size)); + hclumultx2->Fill(eta,float(sizeX)); + hclumulty2->Fill(eta,float(sizeY)); + hcluchar2->Fill(eta,float(charge)); + } + + numOfClustersPerDet2++; + numOfClustersPerLay2++; + numOfPixelsPerLay2 += size; + + } else if(layer==3) { + + hDetMap3->Fill(float(zindex),float(ladder)); + hcluDetMap3->Fill(col,row); + hcharge3->Fill(charge); + hcols3->Fill(col); + hrows3->Fill(row); + + hclusMap3->Fill(gZ,phi); + + if(pt>CLU_SIZE_PT_CUT) { + hsize3->Fill(float(size)); + hsizex3->Fill(float(sizeX)); + hsizey3->Fill(float(sizeY)); + hclumult3->Fill(eta,float(size)); + hclumultx3->Fill(eta,float(sizeX)); + hclumulty3->Fill(eta,float(sizeY)); + hcluchar3->Fill(eta,float(charge)); + } + + numOfClustersPerDet3++; + numOfClustersPerLay3++; + numOfPixelsPerLay3 += size; + } // if layer + + } // if valid + + } // clusters + + if(pixelHits>0) countPixTracks++; + + if(PRINT) cout<<" Clusters for track "<0) { + hclusPerTrk1->Fill(float(numOfClustersPerDet1)); + if(PRINT) cout<<"Lay1: number of clusters per track = "<Fill(float(numOfClustersPerDet2)); + if(PRINT) cout<<"Lay2: number of clusters per track = "<Fill(float(numOfClustersPerDet3)); + if(PRINT) cout<<"Lay3: number of clusters per track = "<0) { + hclusPerLay1->Fill(float(numOfClustersPerLay1)); + hclusPerLay2->Fill(float(numOfClustersPerLay2)); + hclusPerLay3->Fill(float(numOfClustersPerLay3)); + //hdetsPerLay1->Fill(float(numberOfDetUnits1)); + //hdetsPerLay2->Fill(float(numberOfDetUnits2)); + //hdetsPerLay3->Fill(float(numberOfDetUnits3)); + //int tmp = numberOfDetUnits1+numberOfDetUnits2+numberOfDetUnits3; + hpixPerLay1->Fill(float(numOfPixelsPerLay1)); + hpixPerLay2->Fill(float(numOfPixelsPerLay2)); + hpixPerLay3->Fill(float(numOfPixelsPerLay3)); + //htest7->Fill(float(tmp)); + hclusBpix->Fill(float(numberOfClusters)); + hpixBpix->Fill(float(numberOfPixels)); + + } + htracksGood->Fill(float(countNiceTracks)); + htracksGoodInPix->Fill(float(countPixTracks)); + htracks->Fill(float(trackNumber)); + + hbx->Fill(float(bx)); + hlumi->Fill(float(lumiBlock)); + + htracksls->Fill(float(lumiBlock),float(countPixTracks)); + hpvsls->Fill(float(lumiBlock),float(pvsTrue)); + if(instlumi>0.) { + float tmp = float(countPixTracks)/instlumi; + htrackslsn->Fill(float(lumiBlock),tmp); + tmp = float(pvsTrue)/instlumi; + hpvslsn->Fill(float(lumiBlock),tmp); + } + +#endif // HISTOS + + // + countTracks += float(trackNumber); + countGoodTracks += float(countNiceTracks); + countTracksInPix += float(countPixTracks); + countPVs += float(pvsTrue); + countEvents++; + if(lumiBlock != lumiBlockOld) { + countLumi += intlumi; + lumiBlockOld = lumiBlock; + } + + + if(PRINT) cout<<" event with tracks = "< trajTrackCollectionHandle; + e.getByLabel(conf_.getParameter("trajectoryInput"),trajTrackCollectionHandle); + + TrajectoryStateCombiner tsoscomb; + + int NbrTracks = trajTrackCollectionHandle->size(); + std::cout << " track measurements " << trajTrackCollectionHandle->size() << std::endl; + + int trackNumber = 0; + int numberOfClusters = 0; + + for(TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(), itEnd = trajTrackCollectionHandle->end(); it!=itEnd;++it){ + + int pixelHits = 0; + int stripHits = 0; + const Track& track = *it->val; + const Trajectory& traj = *it->key; + + std::vector checkColl = traj.measurements(); + for(std::vector::const_iterator checkTraj = checkColl.begin(); + checkTraj != checkColl.end(); ++checkTraj) { + + if(! checkTraj->updatedState().isValid()) continue; + TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit(); + if(! testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker ) continue; + uint testSubDetID = (testhit->geographicalId().subdetId()); + if(testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap) pixelHits++; + else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB || + testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC) stripHits++; + + } + + + if (pixelHits == 0) continue; + + trackNumber++; + std::cout << " track " << trackNumber <<" has pixelhits "< tmColl = traj.measurements(); + for(std::vector::const_iterator itTraj = checkColl.begin(); + itTraj != checkColl.end(); ++itTraj) { + if(! itTraj->updatedState().isValid()) continue; + + TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() ); + TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit(); + if(! hit->isValid() || hit->geographicalId().det() != DetId::Tracker ) continue; + + const DetId & hit_detId = hit->geographicalId(); + uint IntSubDetID = (hit_detId.subdetId()); + + if(IntSubDetID == 0 ) continue; // Select ?? + if(IntSubDetID != PixelSubdetector::PixelBarrel) continue; // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) { + + +// const GeomDetUnit* detUnit = hit->detUnit(); +// if(detUnit) { +// const Surface& surface = hit->detUnit()->surface(); +// const TrackerGeometry& theTracker(*tkGeom_); +// const PixelGeomDetUnit* theGeomDet = dynamic_cast (theTracker.idToDet(hit_detId) ); +// const RectangularPixelTopology * topol = dynamic_cast(&(theGeomDet->specificTopology())); +// } + + + // get the enclosed persistent hit + const TrackingRecHit *persistentHit = hit->hit(); + // check if it's not null, and if it's a valid pixel hit + if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) { + // tell the C++ compiler that the hit is a pixel hit + const SiPixelRecHit* pixhit = dynamic_cast( hit->hit() ); + // get the edm::Ref to the cluster + edm::Ref, SiPixelCluster> const& clust = (*pixhit).cluster(); + // check if the ref is not null + if (clust.isNonnull()) { + + numberOfClusters++; + pixelHits++; + float charge = (clust->charge())/1000.0; // convert electrons to kilo-electrons + int size = clust->size(); + int size_x = clust->sizeX(); + int size_y = clust->sizeY(); + float row = clust->x(); + float col = clust->y(); + + //LocalPoint lp = topol->localPosition(MeasurementPoint(clust_.row,clust_.col)); + //float x = lp.x(); + //float y = lp.y(); + + int maxPixelCol = clust->maxPixelCol(); + int maxPixelRow = clust->maxPixelRow(); + int minPixelCol = clust->minPixelCol(); + int minPixelRow = clust->minPixelRow(); + + //int geoId = PixGeom->geographicalId().rawId(); + + // Replace with the topology methods + // edge method moved to topologi class + //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) ); + //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) ); + + // calculate alpha and beta from cluster position + //LocalTrajectoryParameters ltp = tsos.localParameters(); + //LocalVector localDir = ltp.momentum()/ltp.momentum().mag(); + + //float locx = localDir.x(); + //float locy = localDir.y(); + //float locz = localDir.z(); + //float loctheta = localDir.theta(); // currently unused + + //float alpha = atan2( locz, locx ); + //float beta = atan2( locz, locy ); + + //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)); + } // valid cluster + } // valid peristant hit + + } // loop over trajectory meas. + + if(PRINT) cout<<" Cluster for track "< +#include +#include + +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +//#include "DataFormats/Common/interface/Handle.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "DataFormats/Common/interface/EDProduct.h" + +//#include "DataFormats/SiPixelCluster/interface/SiPixelClusterCollection.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/DetId/interface/DetId.h" + +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" + +#include "DataFormats/SiPixelDetId/interface/PixelBarrelName.h" + + +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" +//#include "Geometry/CommonTopologies/interface/PixelTopology.h" + +// For L1 +#include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h" +#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h" + +// For HLT +#include "DataFormats/HLTReco/interface/TriggerEvent.h" +#include "DataFormats/HLTReco/interface/TriggerTypeDefs.h" +#include "DataFormats/Common/interface/TriggerResults.h" +#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h" +#include "FWCore/Common/interface/TriggerNames.h" + + +// For tracks +#include "DataFormats/TrackReco/interface/Track.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" + +//#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" + +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" +//#include "TrackingTools/PatternTools/interface/TrajectoryFitter.h" + +#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" + +//#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h" +//#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h" + +//#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h" + + +//#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h" +//#include "TrackingTools/PatternTools/interface/Trajectory.h" +//#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +//#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +//#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include + +// For luminisoty +#include "DataFormats/Luminosity/interface/LumiSummary.h" +#include "DataFormats/Common/interface/ConditionsInEdm.h" + +// To use root histos +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +// For ROOT +#include +#include +#include +#include +#include +#include +#include +#include + + +#define HISTOS +//#define L1 +//#define HLT + +#define VDM_STUDIES + +const bool isData = false; // set false for MC + +using namespace std; + +class TestWithTracks : public edm::EDAnalyzer { + public: + + explicit TestWithTracks(const edm::ParameterSet& conf); + virtual ~TestWithTracks(); + virtual void analyze(const edm::Event& e, const edm::EventSetup& c); + virtual void beginRun(const edm::EventSetup& iSetup); + virtual void beginJob(); + virtual void endJob(); + + private: + edm::ParameterSet conf_; + edm::InputTag src_; + //const static bool PRINT = false; + bool PRINT; + float countTracks, countGoodTracks, countTracksInPix, countPVs, countEvents, countLumi; + + TH1D *hcharge1,*hcharge2, *hcharge3, *hcharge4, *hcharge5; // ADD FPIX + TH1D *hsize1,*hsize2,*hsize3, + *hsizex1,*hsizex2,*hsizex3, + *hsizey1,*hsizey2,*hsizey3; + TH1D *hsize4,*hsize5, // ADD FPIX + *hsizex4,*hsizex5, + *hsizey4,*hsizey5; + + TH1D *hclusPerTrk1,*hclusPerTrk2,*hclusPerTrk3,*hclusPerTrk4,*hclusPerTrk5; + TH1D *hclusPerLay1,*hclusPerLay2,*hclusPerLay3; + TH1D *hclusPerDisk1,*hclusPerDisk2,*hclusPerDisk3,*hclusPerDisk4; + TH1D *hclusPerTrk,*hclusPerTrkB,*hclusPerTrkF; + + TH2F *hDetMap1, *hDetMap2, *hDetMap3; // clusters + TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3; // MODULE PROJECTION + + TH2F *hpvxy, *hclusMap1, *hclusMap2, *hclusMap3; // Z vs PHI + + TH1D *hpvz, *hpvr, *hNumPv, *hNumPvClean; + TH1D *hPt, *hEta, *hDz, *hD0,*hzdiff; + + TH1D *hl1a, *hl1t, *hlt1; + + TH1D *hclusBpix, *hpixBpix; + TH1D *htracks, *htracksGood, *htracksGoodInPix; + + TProfile *hclumult1, *hclumult2, *hclumult3; + TProfile *hclumultx1, *hclumultx2, *hclumultx3; + TProfile *hclumulty1, *hclumulty2, *hclumulty3; + TProfile *hcluchar1, *hcluchar2, *hcluchar3; + + TProfile *hclumultld1, *hclumultld2, *hclumultld3; + TProfile *hclumultxld1, *hclumultxld2, *hclumultxld3; + TProfile *hclumultyld1, *hclumultyld2, *hclumultyld3; + TProfile *hclucharld1, *hclucharld2, *hclucharld3; + + TProfile *htracksls, *hpvsls, *htrackslsn, *hpvslsn, *hintgl, *hinstl, *hbeam1, *hbeam2; + TProfile *hmult1, *hmult2, *hmult3; + TProfile *hclusPerTrkVsEta,*hclusPerTrkVsPt,*hclusPerTrkVsls,*hclusPerTrkVsEtaB,*hclusPerTrkVsEtaF; + + TH1D *hlumi, *hlumi0, *hbx, *hbx0; + +#ifdef VDM_STUDIES + TProfile *hcharCluls, *hcharPixls, *hsizeCluls, *hsizeXCluls; + TProfile *hcharCluls1, *hcharPixls1, *hsizeCluls1, *hsizeXCluls1; + TProfile *hcharCluls2, *hcharPixls2, *hsizeCluls2, *hsizeXCluls2; + TProfile *hcharCluls3, *hcharPixls3, *hsizeCluls3, *hsizeXCluls3; + TProfile *hclusls; // *hpixls; + TProfile *hclusls1, *hclusls2,*hclusls3; + TProfile *hclubx, *hpvbx, *htrackbx; // *hcharClubx, *hcharPixbx,*hsizeClubx, *hsizeYClubx; +#endif + +}; +///////////////////////////////////////////////////////////////// +// Contructor, +TestWithTracks::TestWithTracks(edm::ParameterSet const& conf) +// : conf_(conf), src_(conf.getParameter( "src" )) { } + : conf_(conf) { + + PRINT = conf.getUntrackedParameter("Verbosity",false); + src_ = conf.getParameter( "src" ); + //if(PRINT) cout<<" Construct "< fs; + + int sizeH=20; + float lowH = -0.5; + float highH = 19.5; + hclusPerTrk1 = fs->make( "hclusPerTrk1", "Clus per track l1", + sizeH, lowH, highH); + hclusPerTrk2 = fs->make( "hclusPerTrk2", "Clus per track l2", + sizeH, lowH, highH); + hclusPerTrk3 = fs->make( "hclusPerTrk3", "Clus per track l3", + sizeH, lowH, highH); + hclusPerTrk4 = fs->make( "hclusPerTrk4", "Clus per track d1", + sizeH, lowH, highH); + hclusPerTrk5 = fs->make( "hclusPerTrk5", "Clus per track d2", + sizeH, lowH, highH); + hclusPerTrk = fs->make( "hclusPerTrk", "Clus per track", + sizeH, lowH, highH); + hclusPerTrkB = fs->make( "hclusPerTrkB", "B Clus per track", + sizeH, lowH, highH); + hclusPerTrkF = fs->make( "hclusPerTrkF", "F Clus per track", + sizeH, lowH, highH); + + sizeH=2000; + highH = 1999.5; + hclusPerLay1 = fs->make( "hclusPerLay1", "Clus per layer l1", + sizeH, lowH, highH); + hclusPerLay2 = fs->make( "hclusPerLay2", "Clus per layer l2", + sizeH, lowH, highH); + hclusPerLay3 = fs->make( "hclusPerLay3", "Clus per layer l3", + sizeH, lowH, highH); + hclusPerDisk1 = fs->make( "hclusPerDisk1", "Clus per disk 1", + sizeH, lowH, highH); + hclusPerDisk2 = fs->make( "hclusPerDisk2", "Clus per disk 2", + sizeH, lowH, highH); + hclusPerDisk3 = fs->make( "hclusPerDisk3", "Clus per disk 3", + sizeH, lowH, highH); + hclusPerDisk4 = fs->make( "hclusPerDisk4", "Clus per disk 4", + sizeH, lowH, highH); + + + hcharge1 = fs->make( "hcharge1", "Clu charge l1", 400, 0.,200.); //in ke + hcharge2 = fs->make( "hcharge2", "Clu charge l2", 400, 0.,200.); + hcharge3 = fs->make( "hcharge3", "Clu charge l3", 400, 0.,200.); + hcharge4 = fs->make( "hcharge4", "Clu charge d1", 400, 0.,200.); + hcharge5 = fs->make( "hcharge5", "Clu charge d2", 400, 0.,200.); + + hsize1 = fs->make( "hsize1", "layer 1 clu size",100,-0.5,99.5); + hsize2 = fs->make( "hsize2", "layer 2 clu size",100,-0.5,99.5); + hsize3 = fs->make( "hsize3", "layer 3 clu size",100,-0.5,99.5); + hsizex1 = fs->make( "hsizex1", "lay1 clu size in x", + 20,-0.5,19.5); + hsizex2 = fs->make( "hsizex2", "lay2 clu size in x", + 20,-0.5,19.5); + hsizex3 = fs->make( "hsizex3", "lay3 clu size in x", + 20,-0.5,19.5); + hsizey1 = fs->make( "hsizey1", "lay1 clu size in y", + 30,-0.5,29.5); + hsizey2 = fs->make( "hsizey2", "lay2 clu size in y", + 30,-0.5,29.5); + hsizey3 = fs->make( "hsizey3", "lay3 clu size in y", + 30,-0.5,29.5); + + hsize4 = fs->make( "hsize4", "disk 1 clu size",100,-0.5,99.5); + hsize5 = fs->make( "hsize5", "disk 2 clu size",100,-0.5,99.5); + hsizex4 = fs->make( "hsizex4", "d1 clu size in x", + 20,-0.5,19.5); + hsizex5 = fs->make( "hsizex5", "d2 clu size in x", + 20,-0.5,19.5); + hsizey4 = fs->make( "hsizey4", "d1 clu size in y", + 30,-0.5,29.5); + hsizey5 = fs->make( "hsizey5", "d2 clu size in y", + 30,-0.5,29.5); + + + hDetMap1 = fs->make( "hDetMap1", "layer 1 clus map", + 9,0.,9.,23,0.,23.); + hDetMap2 = fs->make( "hDetMap2", "layer 2 clus map", + 9,0.,9.,33,0.,33.); + hDetMap3 = fs->make( "hDetMap3", "layer 3 clus map", + 9,0.,9.,45,0.,45.); + //hpixDetMap1 = fs->make( "hpixDetMap1", "pix det layer 1", + // 416,0.,416.,160,0.,160.); + //hpixDetMap2 = fs->make( "hpixDetMap2", "pix det layer 2", + // 416,0.,416.,160,0.,160.); + //hpixDetMap3 = fs->make( "hpixDetMap3", "pix det layer 3", + // 416,0.,416.,160,0.,160.); + hcluDetMap1 = fs->make( "hcluDetMap1", "clu det layer 1", + 416,0.,416.,160,0.,160.); + hcluDetMap2 = fs->make( "hcluDetMap2", "clu det layer 2", + 416,0.,416.,160,0.,160.); + hcluDetMap3 = fs->make( "hcluDetMap3", "clu det layer 3", + 416,0.,416.,160,0.,160.); + + htracksGoodInPix = fs->make( "htracksGoodInPix", "count good tracks in pix",2000,-0.5,1999.5); + htracksGood = fs->make( "htracksGood", "count good tracks",2000,-0.5,1999.5); + htracks = fs->make( "htracks", "count tracks",2000,-0.5,1999.5); + hclusBpix = fs->make( "hclusBpix", "count clus in bpix",200,-0.5,1999.5); + hpixBpix = fs->make( "hpixBpix", "count pixels",200,-0.5,1999.5); + + hpvxy = fs->make( "hpvxy", "pv xy",100,-1.,1.,100,-1.,1.); + hpvz = fs->make( "hpvz", "pv z",1000,-50.,50.); + hpvr = fs->make( "hpvr", "pv r",100,0.,1.); + hNumPv = fs->make( "hNumPv", "num of pv",100,0.,100.); + hNumPvClean = fs->make( "hNumPvClean", "num of pv clean",100,0.,100.); + + hPt = fs->make( "hPt", "pt",120,0.,120.); + hEta = fs->make( "hEta", "eta",50,-2.5,2.5); + hD0 = fs->make( "hD0", "d0",500,0.,5.); + hDz = fs->make( "hDz", "pt",250,-25.,25.); + hzdiff = fs->make( "hzdiff", "PVz-Trackz",200,-10.,10.); + + hl1a = fs->make("hl1a", "l1a", 128,-0.5,127.5); + hl1t = fs->make("hl1t", "l1t", 128,-0.5,127.5); + hlt1 = fs->make("hlt1","hlt1",256,-0.5,255.5); + + hclumult1 = fs->make("hclumult1","cluster size layer 1",60,-3.,3.,0.0,100.); + hclumult2 = fs->make("hclumult2","cluster size layer 2",60,-3.,3.,0.0,100.); + hclumult3 = fs->make("hclumult3","cluster size layer 3",60,-3.,3.,0.0,100.); + + hclumultx1 = fs->make("hclumultx1","cluster x-size layer 1",60,-3.,3.,0.0,100.); + hclumultx2 = fs->make("hclumultx2","cluster x-size layer 2",60,-3.,3.,0.0,100.); + hclumultx3 = fs->make("hclumultx3","cluster x-size layer 3",60,-3.,3.,0.0,100.); + + hclumulty1 = fs->make("hclumulty1","cluster y-size layer 1",60,-3.,3.,0.0,100.); + hclumulty2 = fs->make("hclumulty2","cluster y-size layer 2",60,-3.,3.,0.0,100.); + hclumulty3 = fs->make("hclumulty3","cluster y-size layer 3",60,-3.,3.,0.0,100.); + + hcluchar1 = fs->make("hcluchar1","cluster char layer 1",60,-3.,3.,0.0,1000.); + hcluchar2 = fs->make("hcluchar2","cluster char layer 2",60,-3.,3.,0.0,1000.); + hcluchar3 = fs->make("hcluchar3","cluster char layer 3",60,-3.,3.,0.0,1000.); + + // profiles versus ladder + hclumultld1 = fs->make("hclumultld1","cluster size layer 1",23,-11.5,11.5,0.0,100.); + hclumultld2 = fs->make("hclumultld2","cluster size layer 2",35,-17.5,17.5,0.0,100.); + hclumultld3 = fs->make("hclumultld3","cluster size layer 3",47,-23.5,23.5,0.0,100.); + + hclumultxld1 = fs->make("hclumultxld1","cluster x-size layer 1",23,-11.5,11.5,0.0,100.); + hclumultxld2 = fs->make("hclumultxld2","cluster x-size layer 2",35,-17.5,17.5,0.0,100.); + hclumultxld3 = fs->make("hclumultxld3","cluster x-size layer 3",47,-23.5,23.5,0.0,100.); + + hclumultyld1 = fs->make("hclumultyld1","cluster y-size layer 1",23,-11.5,11.5,0.0,100.); + hclumultyld2 = fs->make("hclumultyld2","cluster y-size layer 2",35,-17.5,17.5,0.0,100.); + hclumultyld3 = fs->make("hclumultyld3","cluster y-size layer 3",47,-23.5,23.5,0.0,100.); + + hclucharld1 = fs->make("hclucharld1","cluster char layer 1",23,-11.5,11.5,0.0,1000.); + hclucharld2 = fs->make("hclucharld2","cluster char layer 2",35,-17.5,17.5,0.0,1000.); + hclucharld3 = fs->make("hclucharld3","cluster char layer 3",47,-23.5,23.5,0.0,1000.); + + hintgl = fs->make("hintgl", "inst lumi vs ls ",1000,0.,3000.,0.0,10000.); + hinstl = fs->make("hinstl", "intg lumi vs ls ",1000,0.,3000.,0.0,100.); + hbeam1 = fs->make("hbeam1", "beam1 vs ls ",1000,0.,3000.,0.0,1000.); + hbeam2 = fs->make("hbeam2", "beam2 vs ls ",1000,0.,3000.,0.0,1000.); + + htracksls = fs->make("htracksls","tracks with pix hits vs ls",1000,0.,3000.,0.0,10000.); + hpvsls = fs->make("hpvsls","pvs vs ls",1000,0.,3000.,0.0,1000.); + htrackslsn = fs->make("htrackslsn","tracks with pix hits/lumi vs ls",1000,0.,3000.,0.0,10000.); + hpvslsn = fs->make("hpvslsn","pvs/lumi vs ls",1000,0.,3000.,0.0,1000.); + + hmult1 = fs->make("hmult1","clu mult layer 1",10,0.,10.,0.0,1000.); + hmult2 = fs->make("hmult2","clu mult layer 2",10,0.,10.,0.0,1000.); + hmult3 = fs->make("hmult3","clu mult layer 3",10,0.,10.,0.0,1000.); + + hclusPerTrkVsEta = fs->make("hclusPerTrkVsEta","clus per trk vs.eta",60,-3.,3.,0.0,100.); + hclusPerTrkVsPt = fs->make("hclusPerTrkVsPt","clus per trk vs.pt",120,0.,120.,0.0,100.); + hclusPerTrkVsls = fs->make("hclusPerTrkVsls","clus per trk vs.ls",300,0.,3000.,0.0,100.); + hclusPerTrkVsEtaF = fs->make("hclusPerTrkVsEtaF","F clus per trk vs.eta",60,-3.,3.,0.0,100.); + hclusPerTrkVsEtaB = fs->make("hclusPerTrkVsEtaB","B clus per trk vs.eta",60,-3.,3.,0.0,100.); + + hlumi0 = fs->make("hlumi0", "lumi", 2000,0,2000.); + hlumi = fs->make("hlumi", "lumi", 2000,0,2000.); + hbx0 = fs->make("hbx0", "bx", 4000,0,4000.); + hbx = fs->make("hbx", "bx", 4000,0,4000.); + + hclusMap1 = fs->make("hclusMap1","clus - lay1",260,-26.,26.,350,-3.5,3.5); + hclusMap2 = fs->make("hclusMap2","clus - lay2",260,-26.,26.,350,-3.5,3.5); + hclusMap3 = fs->make("hclusMap3","clus - lay3",260,-26.,26.,350,-3.5,3.5); + + +#ifdef VDM_STUDIES + highH = 3000.; + sizeH = 1000; + + hclusls = fs->make("hclusls","clus vs ls",sizeH,0.,highH,0.0,30000.); + //hpixls = fs->make("hpixls", "pix vs ls ",sizeH,0.,highH,0.0,100000.); + + hcharCluls = fs->make("hcharCluls","clu char vs ls",sizeH,0.,highH,0.0,100.); + //hcharPixls = fs->make("hcharPixls","pix char vs ls",sizeH,0.,highH,0.0,100.); + hsizeCluls = fs->make("hsizeCluls","clu size vs ls",sizeH,0.,highH,0.0,1000.); + hsizeXCluls= fs->make("hsizeXCluls","clu size-x vs ls",sizeH,0.,highH,0.0,100.); + + hcharCluls1 = fs->make("hcharCluls1","clu char vs ls",sizeH,0.,highH,0.0,100.); + //hcharPixls1 = fs->make("hcharPixls1","pix char vs ls",sizeH,0.,highH,0.0,100.); + hsizeCluls1 = fs->make("hsizeCluls1","clu size vs ls",sizeH,0.,highH,0.0,1000.); + hsizeXCluls1= fs->make("hsizeXCluls1","clu size-x vs ls",sizeH,0.,highH,0.0,100.); + hcharCluls2 = fs->make("hcharCluls2","clu char vs ls",sizeH,0.,highH,0.0,100.); + //hcharPixls2 = fs->make("hcharPixls2","pix char vs ls",sizeH,0.,highH,0.0,100.); + hsizeCluls2 = fs->make("hsizeCluls2","clu size vs ls",sizeH,0.,highH,0.0,1000.); + hsizeXCluls2= fs->make("hsizeXCluls2","clu size-x vs ls",sizeH,0.,highH,0.0,100.); + hcharCluls3 = fs->make("hcharCluls3","clu char vs ls",sizeH,0.,highH,0.0,100.); + //hcharPixls3 = fs->make("hcharPixls3","pix char vs ls",sizeH,0.,highH,0.0,100.); + hsizeCluls3 = fs->make("hsizeCluls3","clu size vs ls",sizeH,0.,highH,0.0,1000.); + hsizeXCluls3= fs->make("hsizeXCluls3","clu size-x vs ls",sizeH,0.,highH,0.0,100.); + hclusls1 = fs->make("hclusls1","clus vs ls",sizeH,0.,highH,0.0,30000.); + //hpixls1 = fs->make("hpixls1", "pix vs ls ",sizeH,0.,highH,0.0,100000.); + hclusls2 = fs->make("hclusls2","clus vs ls",sizeH,0.,highH,0.0,30000.); + //hpixls2 = fs->make("hpixls2", "pix vs ls ",sizeH,0.,highH,0.0,100000.); + hclusls3 = fs->make("hclusls3","clus vs ls",sizeH,0.,highH,0.0,30000.); + //hpixls3 = fs->make("hpixls3", "pix vs ls ",sizeH,0.,highH,0.0,100000.); + + // Profiles versus bx + //hpixbx = fs->make("hpixbx", "pixs vs bx ",4000,-0.5,3999.5,0.0,1000000.); + hclubx = fs->make("hclubx", "clus vs bx ",4000,-0.5,3999.5,0.0,1000000.); + hpvbx = fs->make("hpvbx", "pvs vs bx ", 4000,-0.5,3999.5,0.0,1000000.); + htrackbx = fs->make("htrackbx", "tracks vs bx ", 4000,-0.5,3999.5,0.0,1000000.); + +#endif + +#endif + +} +// ------------ method called to at the end of the job ------------ +void TestWithTracks::endJob(){ + cout << " End PixelTracksTest, events = " << countEvents << endl; + + if(countEvents>0.) { + countTracks /= countEvents; + countGoodTracks /= countEvents; + countTracksInPix /= countEvents; + countPVs /= countEvents; + countLumi /= 1000.; + cout<<" Average tracks/event "<< countTracks<<" good "<< countGoodTracks<<" in pix "<< countTracksInPix + <<" PVs "<< countPVs<<" events " << countEvents<<" lumi pb-1 "<< countLumi<<"/10, bug!"<Fill(float(bx)); + hlumi0->Fill(float(lumiBlock)); + + edm::LuminosityBlock const& iLumi = e.getLuminosityBlock(); + edm::Handle lumi; + iLumi.getByLabel("lumiProducer", lumi); + edm::Handle cond; + float intlumi = 0, instlumi=0; + int beamint1=0, beamint2=0; + iLumi.getByLabel("conditionsInEdm", cond); + // This will only work when running on RECO until (if) they fix it in the FW + // When running on RAW and reconstructing, the LumiSummary will not appear + // in the event before reaching endLuminosityBlock(). Therefore, it is not + // possible to get this info in the event + if (lumi.isValid()) { + intlumi =(lumi->intgRecLumi())/1000.; // 10^30 -> 10^33/cm2/sec -> 1/nb/sec + instlumi=(lumi->avgInsDelLumi())/1000.; + beamint1=(cond->totalIntensityBeam1)/1000; + beamint2=(cond->totalIntensityBeam2)/1000; + } else { + //std::cout << "** ERROR: Event does not get lumi info\n"; + } + + hinstl->Fill(float(lumiBlock),float(instlumi)); + hintgl->Fill(float(lumiBlock),float(intlumi)); + hbeam1->Fill(float(lumiBlock),float(beamint1)); + hbeam2->Fill(float(lumiBlock),float(beamint2)); + + +#ifdef L1 + // Get L1 + Handle L1GTRR; + e.getByLabel("gtDigis",L1GTRR); + + if (L1GTRR.isValid()) { + //bool l1a = L1GTRR->decision(); // global decission? + //cout<<" L1 status :"<decisionWord().size(); ++i) { + int l1flag = L1GTRR->decisionWord()[i]; + int t1flag = L1GTRR->technicalTriggerWord()[i]; + + if( l1flag>0 ) hl1a->Fill(float(i)); + if( t1flag>0 && i<64) hl1t->Fill(float(i)); + } // for loop + } // if l1a +#endif + +#ifdef HLT + + bool hlt[256]; + for(int i=0;i<256;++i) hlt[i]=false; + + edm::TriggerNames TrigNames; + edm::Handle HLTResults; + + // Extract the HLT results + e.getByLabel(edm::InputTag("TriggerResults","","HLT"),HLTResults); + if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) { + + //TrigNames.init(*HLTResults); + const edm::TriggerNames & TrigNames = e.triggerNames(*HLTResults); + + //cout<wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) && + (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) && + (HLTResults->error( TrigNames.triggerIndex(TrigNames.triggerName(i))) == false) ) { + + hlt[i]=true; + hlt1->Fill(float(i)); + + } // if hlt + + } // loop + } // if valid +#endif + + + // -- Does this belong into beginJob()? + //ESHandle TG; + //iSetup.get().get(TG); + //const TrackerGeometry* theTrackerGeometry = TG.product(); + //const TrackerGeometry& theTracker(*theTrackerGeometry); + // Get event setup + edm::ESHandle geom; + es.get().get( geom ); + const TrackerGeometry& theTracker(*geom); + + + // -- Primary vertices + // ---------------------------------------------------------------------- + edm::Handle vertices; + e.getByLabel("offlinePrimaryVertices", vertices); + + if(PRINT) cout<<" PV list "<size()< pvzVector; + for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) { + + if( (iv->isFake()) == 1 ) continue; + pvNotFake++; + float pvx = iv->x(); + float pvy = iv->y(); + float pvz = iv->z(); + int numTracksPerPV = iv->tracksSize(); + //int numTracksPerPV = iv->nTracks(); + + //float xe = iv->xError(); + //float ye = iv->yError(); + //float ze = iv->zError(); + //int chi2 = iv->chi2(); + //int dof = iv->ndof(); + + if(PRINT) cout<<" PV "<Fill(pvz); + if(pvz>-22. && pvz<22.) { + float pvr = sqrt(pvx*pvx + pvy*pvy); + hpvxy->Fill(pvx,pvy); + hpvr->Fill(pvr); + if(pvr<0.3) { + pvsTrue++; + pvzVector.push_back(pvz); + //if(PRINT) cout<<"PV "<Fill(float(pvNotFake)); + hNumPvClean->Fill(float(pvsTrue)); + + if(PRINT) cout<<" Not fake PVs = "< recTracks; + // e.getByLabel("generalTracks", recTracks); + // e.getByLabel("ctfWithMaterialTracksP5", recTracks); + // e.getByLabel("splittedTracksP5", recTracks); + //e.getByLabel("cosmictrackfinderP5", recTracks); + e.getByLabel(src_ , recTracks); + + + if(PRINT) cout<<" Tracks "<size()<begin(); + t!=recTracks->end(); ++t){ + + trackNumber++; + numOfClusPerTrk1=0; // this is confusing, it is used as clus per track + numOfClusPerTrk2=0; + numOfClusPerTrk3=0; + numOfClusPerTrk4=0; + numOfClusPerTrk5=0; + int pixelHits=0; + + int size = t->recHitsSize(); + float pt = t->pt(); + float eta = t->eta(); + float phi = t->phi(); + //float trackCharge = t->charge(); // unused + float d0 = t->d0(); + float dz = t->dz(); + //float tkvx = t->vx(); // unused + //float tkvy = t->vy(); + //float tkvz = t->vz(); + + if(PRINT) cout<<"Track "<Fill(eta); + hDz->Fill(dz); + if(abs(eta)>2.8 || abs(dz)>25.) continue; // skip + + hD0->Fill(d0); + if(d0>1.0) continue; // skip + + bool goodTrack = false; + for(vector::iterator m=pvzVector.begin(); m!=pvzVector.end();++m) { + float z = *m; + float tmp = abs(z-dz); + hzdiff->Fill(tmp); + if( tmp < 1.) goodTrack=true; + } + + if(isData && !goodTrack) continue; + countNiceTracks++; + hPt->Fill(pt); + + // Loop over rechits + for ( trackingRecHit_iterator recHit = t->recHitsBegin(); + recHit != t->recHitsEnd(); ++recHit ) { + + if ( !((*recHit)->isValid()) ) continue; + + if( (*recHit)->geographicalId().det() != DetId::Tracker ) continue; + + const DetId & hit_detId = (*recHit)->geographicalId(); + uint IntSubDetID = (hit_detId.subdetId()); + + // Select pixel rechits + if(IntSubDetID == 0 ) continue; // Select ?? + + + int layer=0, ladder=0, zindex=0, ladderOn=0, module=0, shell=0; + unsigned int disk=0; //1,2,3 + unsigned int blade=0; //1-24 + unsigned int zindexF=0; // + unsigned int side=0; //size=1 for -z, 2 for +z + unsigned int panel=0; //panel=1 + + if(IntSubDetID == PixelSubdetector::PixelBarrel) { // bpix + + // Pixel detector + PXBDetId pdetId = PXBDetId(hit_detId); + //unsigned int detTypeP=pdetId.det(); + //unsigned int subidP=pdetId.subdetId(); + // Barell layer = 1,2,3 + layer=pdetId.layer(); + // Barrel ladder id 1-20,32,44. + ladder=pdetId.ladder(); + // Barrel Z-index=1,8 + zindex=pdetId.module(); + + // Convert to online + PixelBarrelName pbn(pdetId); + // Shell { mO = 1, mI = 2 , pO =3 , pI =4 }; + PixelBarrelName::Shell sh = pbn.shell(); //enum + //sector = pbn.sectorName(); + ladderOn = pbn.ladderName(); + //layerOn = pbn.layerName(); + module = pbn.moduleName(); // 1 to 4 + //half = pbn.isHalfModule(); + shell = int(sh); + // change the module sign for z<0 + if(shell==1 || shell==2) module = -module; // make -1 to -4 for -z + // change ladeer sign for Outer )x<0) + if(shell==1 || shell==3) ladderOn = -ladderOn; + + if(PRINT) cout<<"barrel layer/ladder/module: "< (theTracker.idToDet(hit_detId) ); + //double detZ = theGeomDet->surface().position().z(); // unused + //double detR = theGeomDet->surface().position().perp(); // unused + const PixelTopology * topol = &(theGeomDet->specificTopology()); // pixel topology + + //std::vector output = getRecHitComponents((*recHit).get()); + //std::vector TrkComparison::getRecHitComponents(const TrackingRecHit* rechit){ + + const SiPixelRecHit* hit = dynamic_cast((*recHit).get()); + //edm::Ref ,SiStripCluster> cluster = hit->cluster(); + // get the edm::Ref to the cluster + + if(hit) { + + // RecHit (recthits are transient, so not available without ttrack refit) + //double xloc = hit->localPosition().x();// 1st meas coord + //double yloc = hit->localPosition().y();// 2nd meas coord or zero + //double zloc = hit->localPosition().z();// up, always zero + //cout<<" rechit loc "<, SiPixelCluster> const& clust = hit->cluster(); + // check if the ref is not null + if (!clust.isNonnull()) continue; + + numberOfClusters++; + pixelHits++; + float charge = (clust->charge())/1000.0; // convert electrons to kilo-electrons + int size = clust->size(); + int sizeX = clust->sizeX(); + int sizeY = clust->sizeY(); + float row = clust->x(); + float col = clust->y(); + numberOfPixels += size; + + //cout<<" clus loc "<surface().toGlobal( lp ); + double gX = clustgp.x(); + double gY = clustgp.y(); + double gZ = clustgp.z(); + + //cout<<" CLU GLOBAL "<maxPixelCol(); + //int maxPixelRow = clust->maxPixelRow(); + //int minPixelCol = clust->minPixelCol(); + //int minPixelRow = clust->minPixelRow(); + //int geoId = PixGeom->geographicalId().rawId(); + // Replace with the topology methods + // edge method moved to topologi class + //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) ); + //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) ); + + // calculate alpha and beta from cluster position + //LocalTrajectoryParameters ltp = tsos.localParameters(); + //LocalVector localDir = ltp.momentum()/ltp.momentum().mag(); + //float locx = localDir.x(); + //float locy = localDir.y(); + //float locz = localDir.z(); + //float loctheta = localDir.theta(); // currently unused + //float alpha = atan2( locz, locx ); + //float beta = atan2( locz, locy ); + + if(layer==1) { + + hDetMap1->Fill(float(zindex),float(ladder)); + hcluDetMap1->Fill(col,row); + hcharge1->Fill(charge); + //hcols1->Fill(col); + //hrows1->Fill(row); + + hclusMap1->Fill(gZ,phi); + hmult1->Fill(zindex,float(size)); + + if(pt>CLU_SIZE_PT_CUT) { + hsize1->Fill(float(size)); + hsizex1->Fill(float(sizeX)); + hsizey1->Fill(float(sizeY)); + + hclumult1->Fill(eta,float(size)); + hclumultx1->Fill(eta,float(sizeX)); + hclumulty1->Fill(eta,float(sizeY)); + hcluchar1->Fill(eta,float(charge)); + + //cout<Fill(float(ladderOn),size); + hclumultxld1->Fill(float(ladderOn),sizeX); + hclumultyld1->Fill(float(ladderOn),sizeY); + hclucharld1->Fill(float(ladderOn),charge); + + } + +#ifdef VDM_STUDIES + hcharCluls->Fill(lumiBlock,charge); + hsizeCluls->Fill(lumiBlock,size); + hsizeXCluls->Fill(lumiBlock,sizeX); + hcharCluls1->Fill(lumiBlock,charge); + hsizeCluls1->Fill(lumiBlock,size); + hsizeXCluls1->Fill(lumiBlock,sizeX); +#endif + + numOfClusPerTrk1++; + numOfClustersPerLay1++; + numOfPixelsPerLay1 += size; + + } else if(layer==2) { + + hDetMap2->Fill(float(zindex),float(ladder)); + hcluDetMap2->Fill(col,row); + hcharge2->Fill(charge); + //hcols2->Fill(col); + //hrows2->Fill(row); + + hclusMap2->Fill(gZ,phi); + hmult2->Fill(zindex,float(size)); + + if(pt>CLU_SIZE_PT_CUT) { + hsize2->Fill(float(size)); + hsizex2->Fill(float(sizeX)); + hsizey2->Fill(float(sizeY)); + + hclumult2->Fill(eta,float(size)); + hclumultx2->Fill(eta,float(sizeX)); + hclumulty2->Fill(eta,float(sizeY)); + hcluchar2->Fill(eta,float(charge)); + + hclumultld2->Fill(float(ladderOn),size); + hclumultxld2->Fill(float(ladderOn),sizeX); + hclumultyld2->Fill(float(ladderOn),sizeY); + hclucharld2->Fill(float(ladderOn),charge); + + } + +#ifdef VDM_STUDIES + hcharCluls->Fill(lumiBlock,charge); + hsizeCluls->Fill(lumiBlock,size); + hsizeXCluls->Fill(lumiBlock,sizeX); + hcharCluls2->Fill(lumiBlock,charge); + hsizeCluls2->Fill(lumiBlock,size); + hsizeXCluls2->Fill(lumiBlock,sizeX); +#endif + + numOfClusPerTrk2++; + numOfClustersPerLay2++; + numOfPixelsPerLay2 += size; + + } else if(layer==3) { + + hDetMap3->Fill(float(zindex),float(ladder)); + hcluDetMap3->Fill(col,row); + hcharge3->Fill(charge); + //hcols3->Fill(col); + //hrows3->Fill(row); + + hclusMap3->Fill(gZ,phi); + hmult3->Fill(zindex,float(size)); + + if(pt>CLU_SIZE_PT_CUT) { + hsize3->Fill(float(size)); + hsizex3->Fill(float(sizeX)); + hsizey3->Fill(float(sizeY)); + hclumult3->Fill(eta,float(size)); + hclumultx3->Fill(eta,float(sizeX)); + hclumulty3->Fill(eta,float(sizeY)); + hcluchar3->Fill(eta,float(charge)); + + hclumultld3->Fill(float(ladderOn),size); + hclumultxld3->Fill(float(ladderOn),sizeX); + hclumultyld3->Fill(float(ladderOn),sizeY); + hclucharld3->Fill(float(ladderOn),charge); + } + +#ifdef VDM_STUDIES + hcharCluls->Fill(lumiBlock,charge); + hsizeCluls->Fill(lumiBlock,size); + hsizeXCluls->Fill(lumiBlock,sizeX); + hcharCluls3->Fill(lumiBlock,charge); + hsizeCluls3->Fill(lumiBlock,size); + hsizeXCluls3->Fill(lumiBlock,sizeX); +#endif + + numOfClusPerTrk3++; + numOfClustersPerLay3++; + numOfPixelsPerLay3 += size; + + } else if (disk==1) { + + numOfClusPerTrk4++; + numOfClustersPerLay4++; + numOfPixelsPerLay4 += size; + + hcharge4->Fill(charge); + if(pt>CLU_SIZE_PT_CUT) { + hsize4->Fill(float(size)); + hsizex4->Fill(float(sizeX)); + hsizey4->Fill(float(sizeY)); + } + + if(side==1) numOfClustersPerDisk2++; // -z + else if(side==2) numOfClustersPerDisk3++; // +z + + + } else if (disk==2) { + + numOfClusPerTrk5++; + numOfClustersPerLay5++; + numOfPixelsPerLay5 += size; + + hcharge5->Fill(charge); + if(pt>CLU_SIZE_PT_CUT) { + hsize5->Fill(float(size)); + hsizex5->Fill(float(sizeX)); + hsizey5->Fill(float(sizeY)); + } + + if(side==1) numOfClustersPerDisk1++; // -z + else if(side==2) numOfClustersPerDisk4++; // +z + + } else { + + cout<<" which layer is this? "<0) countPixTracks++; + + if(PRINT) cout<<" Clusters for track "<0) { + hclusPerTrk1->Fill(float(numOfClusPerTrk1)); + if(PRINT) cout<<"Lay1: number of clusters per track = "<Fill(float(numOfClusPerTrk2)); + if(PRINT) cout<<"Lay2: number of clusters per track = "<Fill(float(numOfClusPerTrk3)); + if(PRINT) cout<<"Lay3: number of clusters per track = "<Fill(float(numOfClusPerTrk4)); // fpix disk1 + hclusPerTrk5->Fill(float(numOfClusPerTrk5)); // fpix disk2 + + float clusPerTrkB = numOfClusPerTrk1 + numOfClusPerTrk2 + numOfClusPerTrk3; + float clusPerTrkF = numOfClusPerTrk4 + numOfClusPerTrk5; + float clusPerTrk = clusPerTrkB + clusPerTrkF; + + hclusPerTrkB->Fill(clusPerTrkB); + hclusPerTrkF->Fill(clusPerTrkF); + hclusPerTrk->Fill(clusPerTrk); + + hclusPerTrkVsEta->Fill(eta,clusPerTrk); + hclusPerTrkVsEtaB->Fill(eta,clusPerTrkB); + hclusPerTrkVsEtaF->Fill(eta,clusPerTrkF); + hclusPerTrkVsPt->Fill(pt,clusPerTrk); + hclusPerTrkVsls->Fill(lumiBlock,clusPerTrk); + + } +#endif // HISTOS + + + } // tracks + + +#ifdef HISTOS + // total layer histos + if(numberOfClusters>0) { + hclusPerLay1->Fill(float(numOfClustersPerLay1)); + hclusPerLay2->Fill(float(numOfClustersPerLay2)); + hclusPerLay3->Fill(float(numOfClustersPerLay3)); + + hclusPerDisk1->Fill(float(numOfClustersPerDisk1)); + hclusPerDisk2->Fill(float(numOfClustersPerDisk2)); + hclusPerDisk3->Fill(float(numOfClustersPerDisk3)); + hclusPerDisk4->Fill(float(numOfClustersPerDisk4)); + + //hdetsPerLay1->Fill(float(numberOfDetUnits1)); + //hdetsPerLay2->Fill(float(numberOfDetUnits2)); + //hdetsPerLay3->Fill(float(numberOfDetUnits3)); + //int tmp = numberOfDetUnits1+numberOfDetUnits2+numberOfDetUnits3; + //hpixPerLay1->Fill(float(numOfPixelsPerLay1)); + //hpixPerLay2->Fill(float(numOfPixelsPerLay2)); + //hpixPerLay3->Fill(float(numOfPixelsPerLay3)); + //htest7->Fill(float(tmp)); + hclusBpix->Fill(float(numberOfClusters)); + hpixBpix->Fill(float(numberOfPixels)); + + } + htracksGood->Fill(float(countNiceTracks)); + htracksGoodInPix->Fill(float(countPixTracks)); + htracks->Fill(float(trackNumber)); + + hbx->Fill(float(bx)); + hlumi->Fill(float(lumiBlock)); + + htracksls->Fill(float(lumiBlock),float(countPixTracks)); + hpvsls->Fill(float(lumiBlock),float(pvsTrue)); + if(instlumi>0.) { + float tmp = float(countPixTracks)/instlumi; + htrackslsn->Fill(float(lumiBlock),tmp); + tmp = float(pvsTrue)/instlumi; + hpvslsn->Fill(float(lumiBlock),tmp); + } + + +#ifdef VDM_STUDIES + + hclusls->Fill(float(lumiBlock),float(numberOfClusters)); // clusters fpix+bpix + //hpixls->Fill(float(lumiBlock),float(numberOfPixels)); // pixels fpix+bpix + + hclubx->Fill(float(bx),float(numberOfClusters)); // clusters fpix+bpix + //hpixbx->Fill(float(bx),float(numberOfPixels)); // pixels fpix+bpix + hpvbx->Fill(float(bx),float(pvsTrue)); // pvs + htrackbx->Fill(float(bx),float(countPixTracks)); // tracks + + hclusls1->Fill(float(lumiBlock),float(numOfClustersPerLay1)); // clusters bpix1 + //hpixls1->Fill( float(lumiBlock),float(numOfPixPerLay1)); // pixels bpix1 + hclusls2->Fill(float(lumiBlock),float(numOfClustersPerLay2)); // clusters bpix2 + //hpixls2->Fill( float(lumiBlock),float(numOfPixPerLay2)); // pixels bpix2 + hclusls3->Fill(float(lumiBlock),float(numOfClustersPerLay3)); // clusters bpix3 + //hpixls3->Fill( float(lumiBlock),float(numOfPixPerLay3)); // pixels bpix3 +#endif + + + +#endif // HISTOS + + // + countTracks += float(trackNumber); + countGoodTracks += float(countNiceTracks); + countTracksInPix += float(countPixTracks); + countPVs += float(pvsTrue); + countEvents++; + if(lumiBlock != lumiBlockOld) { + countLumi += intlumi; + lumiBlockOld = lumiBlock; + } + + + if(PRINT) cout<<" event with tracks = "< trajTrackCollectionHandle; + e.getByLabel(conf_.getParameter("trajectoryInput"),trajTrackCollectionHandle); + + TrajectoryStateCombiner tsoscomb; + + int NbrTracks = trajTrackCollectionHandle->size(); + std::cout << " track measurements " << trajTrackCollectionHandle->size() << std::endl; + + int trackNumber = 0; + int numberOfClusters = 0; + + for(TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(), itEnd = trajTrackCollectionHandle->end(); it!=itEnd;++it){ + + int pixelHits = 0; + int stripHits = 0; + const Track& track = *it->val; + const Trajectory& traj = *it->key; + + std::vector checkColl = traj.measurements(); + for(std::vector::const_iterator checkTraj = checkColl.begin(); + checkTraj != checkColl.end(); ++checkTraj) { + + if(! checkTraj->updatedState().isValid()) continue; + TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit(); + if(! testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker ) continue; + uint testSubDetID = (testhit->geographicalId().subdetId()); + if(testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap) pixelHits++; + else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB || + testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC) stripHits++; + + } + + + if (pixelHits == 0) continue; + + trackNumber++; + std::cout << " track " << trackNumber <<" has pixelhits "< tmColl = traj.measurements(); + for(std::vector::const_iterator itTraj = checkColl.begin(); + itTraj != checkColl.end(); ++itTraj) { + if(! itTraj->updatedState().isValid()) continue; + + TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() ); + TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit(); + if(! hit->isValid() || hit->geographicalId().det() != DetId::Tracker ) continue; + + const DetId & hit_detId = hit->geographicalId(); + uint IntSubDetID = (hit_detId.subdetId()); + + if(IntSubDetID == 0 ) continue; // Select ?? + if(IntSubDetID != PixelSubdetector::PixelBarrel) continue; // look only at bpix || IntSubDetID == PixelSubdetector::PixelEndcap) { + + +// const GeomDetUnit* detUnit = hit->detUnit(); +// if(detUnit) { +// const Surface& surface = hit->detUnit()->surface(); +// const TrackerGeometry& theTracker(*tkGeom_); +// const PixelGeomDetUnit* theGeomDet = dynamic_cast (theTracker.idToDet(hit_detId) ); +// const RectangularPixelTopology * topol = dynamic_cast(&(theGeomDet->specificTopology())); +// } + + + // get the enclosed persistent hit + const TrackingRecHit *persistentHit = hit->hit(); + // check if it's not null, and if it's a valid pixel hit + if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) { + // tell the C++ compiler that the hit is a pixel hit + const SiPixelRecHit* pixhit = dynamic_cast( hit->hit() ); + // get the edm::Ref to the cluster + edm::Ref, SiPixelCluster> const& clust = (*pixhit).cluster(); + // check if the ref is not null + if (clust.isNonnull()) { + + numberOfClusters++; + pixelHits++; + float charge = (clust->charge())/1000.0; // convert electrons to kilo-electrons + int size = clust->size(); + int size_x = clust->sizeX(); + int size_y = clust->sizeY(); + float row = clust->x(); + float col = clust->y(); + + //LocalPoint lp = topol->localPosition(MeasurementPoint(clust_.row,clust_.col)); + //float x = lp.x(); + //float y = lp.y(); + + int maxPixelCol = clust->maxPixelCol(); + int maxPixelRow = clust->maxPixelRow(); + int minPixelCol = clust->minPixelCol(); + int minPixelRow = clust->minPixelRow(); + + //int geoId = PixGeom->geographicalId().rawId(); + + // Replace with the topology methods + // edge method moved to topologi class + //int edgeHitX = (int) ( topol->isItEdgePixelInX( minPixelRow ) || topol->isItEdgePixelInX( maxPixelRow ) ); + //int edgeHitY = (int) ( topol->isItEdgePixelInY( minPixelCol ) || topol->isItEdgePixelInY( maxPixelCol ) ); + + // calculate alpha and beta from cluster position + //LocalTrajectoryParameters ltp = tsos.localParameters(); + //LocalVector localDir = ltp.momentum()/ltp.momentum().mag(); + + //float locx = localDir.x(); + //float locy = localDir.y(); + //float locz = localDir.z(); + //float loctheta = localDir.theta(); // currently unused + + //float alpha = atan2( locz, locx ); + //float beta = atan2( locz, locy ); + + //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)); + } // valid cluster + } // valid peristant hit + + } // loop over trajectory meas. + + if(PRINT) cout<<" Cluster for track "< +#include +#include + +// ROOT: +#include "TH1.h" +#include "TH2.h" +#include "TProfile.h" + +// CMS and user include files: +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +//#include +#include "FWCore/Framework/interface/MakerMacros.h" + +#include + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +//#include "DataFormats/METReco/interface/PFMET.h" +//#include "DataFormats/METReco/interface/PFMETFwd.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackExtra.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" + +#include + +// To convert detId to subdet/layer number: +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" + +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +//#include "DataFormats/GeometryVector/interface/GlobalPoint.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +//#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include +#include + +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" +#include +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" + +//#define SIM // use for single muon simulations + +// +// class declaration: +// +class myCounters{ +public: + static int neve; + static unsigned int prevrun; +}; + +int myCounters::neve = 0; +unsigned int myCounters::prevrun = 0; +// +// +// +class Triplet : public edm::EDAnalyzer { + +public: + explicit Triplet(const edm::ParameterSet&); + ~Triplet(); + +private: + virtual void beginJob() ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob(); + void triplets(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3, + double ptsig, double & dc,double & dz); + + + // ----------member data: + + TH1D *h000, *h001, *h002, *h003, *h004, *h005, *h006, *h007, *h008, *h009; + TH1D *h010, *h011, *h012, *h013, *h014, *h015, *h016, *h017, *h018, *h019; + TH1D *h020, *h021, *h022, *h023, *h024, *h025, *h026, *h027, *h028, *h029; + TH1D *h030, *h031, *h032, *h033, *h034, *h035, *h036, *h037, *h038, *h039; + TH1D *h040, *h041, *h042, *h043, *h044, *h045, *h046, *h047, *h048, *h049; + TH1D *h050, *h051, *h052, *h053, *h054, *h055, *h056, *h057, *h058, *h059; + TH1D *h060, *h061, *h062, *h063, *h064, *h065, *h067, *h068, *h069; + TH2D *h066; + TH1D *h070, *h071, *h072, *h073, *h074, *h075, *h076, *h077, *h078, *h079; + TH1D *h080, *h081, *h082, *h083, *h084, *h085, *h086, *h087, *h088, *h089; + TH1D *h090, *h091, *h092, *h093, *h094; + TH2D *h095, *h096, *h097, *h098, *h099; + + TH1D *h100, *h101, *h102, *h103, *h104, *h105, *h108; + TH2D *h106, *h107, *h109; + TH1D *h110, *h112, *h113, *h114, *h115, *h116, *h118, *h119; + TH2D *h111, *h117; + TH1D *h120, *h121, *h122, *h123, *h124, *h125, *h126, *h127, *h128, *h129; + TH1D *h130, *h131, *h132, *h133, *h134, *h135, *h136, *h137, *h138, *h139; + TH1D *h140, *h141, *h142, *h143, *h144, *h145, *h146, *h147, *h148, *h149; + TH1D *h150, *h151, *h152, *h153, *h154, *h155, *h156, *h157, *h158, *h159; + TH1D *h160, *h161, *h162, *h163, *h164, *h165, *h166, *h167, *h168, *h169; + TH1D *h170, *h171, *h172, *h173, *h174, *h175, *h176, *h177, *h178, *h179; + TH1D *h180, *h181, *h182, *h183, *h184, *h185, *h186, *h187, *h188, *h189; + TH1D *h190, *h191, *h192, *h193, *h194, *h195, *h196, *h197, *h198, *h199; + + TH1D *h200, *h201, *h202, *h203, *h204, *h205, *h208; + TH2D *h206, *h207, *h209; + + TH1D *h300, *h301, *h302, *h303, *h304, *h305, *h308; + TH2D *h306, *h307, *h309; + + TH1D *h400, *h401, *h402, *h403, *h404, *h405, *h406, *h407, *h408; + TProfile *h409; + TH1D *h410, *h411; + TProfile *h412, *h413, *h414, *h415, *h416, *h417, *h418, *h419; + TH1D *h420, *h421; + TProfile *h422, *h423, *h424, *h425, *h426, *h427, *h428, *h429; + TH1D *h430, *h431, *h432, *h435, *h436, *h437, *h438, *h439; + TProfile *h433, *h434; + TH1D *h440, *h441; + TProfile *h442, *h443, *h444, *h445, *h446, *h447, *h448, *h449; + TH1D *h450, *h451; + +}; + +// +// constants, enums and typedefs: +// + +// +// static data member definitions: +// + +// +// constructor: +// +Triplet::Triplet(const edm::ParameterSet& iConfig) +{ + std::cout << "Triplet constructed\n"; +} +// +// destructor: +// +Triplet::~Triplet() +{ + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) +} +// +// member functions: +// method called once each job just before starting event loop +// +void Triplet::beginJob() +{ + edm::Service fs; + + h000 = fs->make("h000", "beta star; beta star [cm]", 100, 0, 500 ); + h001 = fs->make("h001", "emittance; emittance x [cm]", 100, 0, 1e-5 ); + h002 = fs->make("h002", "beam width x; beam width x [um]", 100, 0, 200 ); + h003 = fs->make("h003", "beam width y; beam width y [um]", 100, 0, 200 ); + h004 = fs->make("h004", "beam spot sigma z; beam spot sigma z [cm]", 100, 0, 20 ); + h005 = fs->make("h005", "beam spot x; beam spot x [cm]", 100, -1, 1 ); + h006 = fs->make("h006", "beam spot y; beam spot y [cm]", 100, -1, 1 ); + h007 = fs->make("h007", "beam spot z; beam spot z [cm]", 100, -5, 5 ); + h008 = fs->make("h008", "beam slope dxdz; beam slope dxdz [mrad]", 100, -5, 5 ); + h009 = fs->make("h009", "beam slope dydz; beam slope dydz [mrad]", 100, -5, 5 ); + + h010 = fs->make("h010", "number of primary vertices; vertices; events", 31, -0.5, 30.5 ); + h011 = fs->make("h011", "invalid z-vertex;z [cm]", 100, -50, 50 ); + h012 = fs->make("h012", "fake z-vertex;z [cm]", 100, -50, 50 ); + h013 = fs->make("h013", "non-fake z-vertex;z [cm]", 100, -50, 50 ); + h014 = fs->make("h014", "vertex x;x [cm]", 100, -0.5, 0.5 ); + h015 = fs->make("h015", "vertex y;y [cm]", 100, -0.5, 0.5 ); + h016 = fs->make("h016", "tracks per vertex; tracks; vertices", 101, -0.5, 100.5 ); + h017 = fs->make("h017", "tracks per vertex; tracks; vertices", 100, 5, 505 ); + h018 = fs->make("h018", "z-vertex with refitted tracks;z [cm]", 100, -50, 50 ); + h019 = fs->make("h019", "z-vertex without refitted tracks;z [cm]", 100, -50, 50 ); + + h021 = fs->make("h021", "vertex sum pt; sum pt [GeV]", 100, 0, 100 ); + h022 = fs->make("h022", "vertex max sum pt; max sum pt [GeV]", 100, 0, 100 ); + + h023 = fs->make("h023", "best vertex x;x [cm]", 100, -0.25, 0.25 ); + h024 = fs->make("h024", "best vertex y;y [cm]", 100, -0.25, 0.25 ); + h025 = fs->make("h025", "best vertex z;z [cm]", 100, -25, 25 ); + + //h026 = fs->make("h026", "Sum Et; Sum Et [GeV]", 100, 0, 1000 ); + //h027 = fs->make("h027", "MET; MET [GeV]", 100, 0, 200 ); + + h028 = fs->make("h028", "sum track pt; sum track Pt [GeV]", 100, 0, 500 ); + h029 = fs->make("h029", "sum primary track charge; sum track charge", 41, -20.5, 20.5 ); + + h040 = fs->make("h040", "number of tracks; tracks", 101, -5, 1005 ); + h041 = fs->make("h041", "track charge; charge", 11, -5.5, 5.5 ); + h042 = fs->make("h042", "pt; pt [GeV]", 100, 0, 5 ); + h043 = fs->make("h043", "pt use logy, pt [GeV]", 100, 0, 100 ); + + h044 = fs->make("h044", "number of rec hits; hits; tracks", 51, -0.5, 50.5 ); + h045 = fs->make("h045", "valid tracker hits; tracker hits; tracks", 51, -0.5, 50.5 ); + h046 = fs->make("h046", "valid pixel barrel hits; valid pixel barrel hits; tracks", 11, -0.5, 10.5 ); + h047 = fs->make("h047", "tracker layers; tracker layers; tracks", 31, -0.5, 30.5 ); + h048 = fs->make("h048", "pixel barrel layers; pixel barrel layers; tracks", 11, -0.5, 10.5 ); + + h051 = fs->make("h051", "kap-kap; dkap; tracks", 100, -0.01, 0.01 ); + h052 = fs->make("h052", "phi-phi; dphi; tracks", 100, -0.1, 0.1 ); + h053 = fs->make("h053", "dca-dca; ddca; tracks", 100, -0.1, 0.1 ); + h054 = fs->make("h054", "dip-dip; ddip; tracks", 100, -0.1, 0.1 ); + h055 = fs->make("h055", "z0-z0; dz0; tracks", 100, -0.1, 0.1 ); + + h056 = fs->make("h056", "tracks", 100, 0.01142, 0.01143 ); + h049 = fs->make("h049", "tracks", 1000, -0.000001, 0.000001); + + h057 = fs->make("h057", "tscp ref x; x [cm]; hits", 100, -1, 1 ); + h058 = fs->make("h058", "tscp ref y; y [cm]; hits", 100, -1, 1 ); + h059 = fs->make("h059", "tscp ref z; z [cm]; hits", 100, -10, 10 ); + + h060 = fs->make("h060", "rec hit tracker subdet; subdet ID; tracks", 11, -0.5, 10.5 ); + h061 = fs->make("h061", "rec hits local x; x [cm]; hits", 120, -6, 6 ); + h062 = fs->make("h062", "rec hits local y; y [cm]; hits", 80, -4, 4 ); + + h063 = fs->make("h063", "rec hits global R; R [cm]; hits", 120, 0, 120 ); + h064 = fs->make("h064", "rec hits global phi; phi [deg]; hits", 180, -180, 180 ); + h065 = fs->make("h065", "rec hits global z; z [cm]; hits", 300, -300, 300 ); + + h066 = fs->make("h066", "rec hits barrel x-y; x [cm]; y [cm]", 240, -120, 120, 240, -120, 120 ); + + h100 = fs->make("h100", "hits on tracks PXB layer; PXB layer; hits", 6, -0.5, 5.5 ); + + h101 = fs->make("h101", "hits on tracks PXB1 ladder; ladder; hits", 22, -0.5, 21.5 ); + h102 = fs->make("h102", "hits on tracks PXB1 module; module; hits", 10, -0.5, 9.5 ); + h103 = fs->make("h103", "hits on tracks PXB1 R; R [cm]; hits", 150, 0, 15 ); + h104 = fs->make("h104", "hits on tracks PXB1 phi; phi [deg]; hits", 180, -180, 180 ); + h105 = fs->make("h105", "hits on tracks PXB1 z; z [cm]; hits", 600, -30, 30 ); + h106 = fs->make("h106", "hits on tracks PXB1 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30 ); + h107 = fs->make("h107", "hits local x-y PXB1; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3 ); + + h111 = fs->make("h111", "hits on tracks PXB1 x-y; x [cm]; y [cm]", 240, -6, 6, 240, -6, 6 ); + h112 = fs->make("h112", "residuals PXB1 dU; dU [um]; hits", 100, -250, 250 ); + h113 = fs->make("h113", "residuals PXB1 dZ; dZ [um]; hits", 100, -250, 250 ); + h114 = fs->make("h114", "residuals PXB1 dU; dU [um]; hits", 100, -250, 250 ); + h115 = fs->make("h115", "residuals PXB1 dZ; dZ [um]; hits", 100, -250, 250 ); + + h201 = fs->make("h201", "hits on tracks PXB2 ladder; ladder; hits", 34, -0.5, 33.5 ); + h202 = fs->make("h202", "hits on tracks PXB2 module; module; hits", 10, -0.5, 9.5 ); + h203 = fs->make("h203", "hits on tracks PXB2 R; R [cm]; hits", 150, 0, 15 ); + h204 = fs->make("h204", "hits on tracks PXB2 phi; phi [deg]; hits", 180, -180, 180 ); + h205 = fs->make("h205", "hits on tracks PXB2 z; z [cm]; hits", 600, -30, 30 ); + h206 = fs->make("h206", "hits on tracks PXB2 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30 ); + h207 = fs->make("h207", "hits local x-y PXB2; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3 ); + + h301 = fs->make("h301", "hits on tracks PXB3 ladder; ladder; hits", 46, -0.5, 45.5 ); + h302 = fs->make("h302", "hits on tracks PXB3 module; module; hits", 10, -0.5, 9.5 ); + h303 = fs->make("h303", "hits on tracks PXB3 R; R [cm]; hits", 150, 0, 15 ); + h304 = fs->make("h304", "hits on tracks PXB3 phi; phi [deg]; hits", 180, -180, 180 ); + h305 = fs->make("h305", "hits on tracks PXB3 z; z [cm]; hits", 600, -30, 30 ); + h306 = fs->make("h306", "hits on tracks PXB3 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30 ); + h307 = fs->make("h307", "hits local x-y PXB3; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3 ); + // + // triplets: + // + h401 = fs->make("h401", "triplets z2; z [cm]; hits", 600, -30, 30 ); + h402 = fs->make("h402", "uphi-phi; dphi; tracks", 100, -0.1, 0.1 ); + h403 = fs->make("h403", "udca-dca; ddca; tracks", 100, -0.1, 0.1 ); + h404 = fs->make("h404", "udip-dip; ddip; tracks", 100, -0.1, 0.1 ); + h405 = fs->make("h405", "uz0-z0; dz0; tracks", 100, -0.1, 0.1 ); + + h406 = fs->make("h406", "valid tracker hits; tracker hits; tracks", 51, -0.5, 50.5 ); + h407 = fs->make("h407", "valid pixel barrel hits; valid pixel barrel hits; tracks", 11, -0.5, 10.5 ); + h408 = fs->make("h408", "tracker layers; tracker layers; tracks", 31, -0.5, 30.5 ); + h409 = fs->make("h409", "angle of incidence; phi2 [deg]; phi inc 2 [deg]", 180, -180, 180, -90, 90 ); + + h410 = fs->make("h410", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150 ); + h411 = fs->make("h411", "residuals PXB2 dz2 ; dz2 [um]; hits", 100, -300, 300 ); + // + // mean resid profiles: + // + h412 = fs->make("h412", "PXB2 dxy vs phi; phi2 [deg]; [um]", 180, -180, 180, -99, 99 ); + h413 = fs->make("h413", "PXB2 dz vs phi; phi2 [deg]; [um]", 180, -180, 180, -99, 199 ); + + h414 = fs->make("h414", "PXB2 dxy vs z; z2 [cm]; [um]", 80, -20, 20, -99, 99 ); + h415 = fs->make("h415", "PXB2 dz vs z; z2 [cm]; [um]", 80, -20, 20, -199, 199 ); + + h416 = fs->make("h416", "PXB2 dxy vs pt; log(pt [GeV]); [um]", 20, 0, 2, -99, 99 ); + h417 = fs->make("h417", "PXB2 dz vs pt; log(pt [GeV]); [um]", 20, 0, 2, -199, 199 ); + + h418 = fs->make("h418", "PXB2 dxy vs pt +; log(pt [GeV]); [um]", 20, 0, 2, -99, 99 ); + h419 = fs->make("h419", "PXB2 dxy vs pt -; log(pt [GeV]); [um]", 20, 0, 2, -99, 99 ); + + h420 = fs->make("h420", "residuals PXB2 dca2, pt > 12; dca2 [um]; hits", 100, -150, 150 ); + h421 = fs->make("h421", "residuals PXB2 dz2, pt > 12; dz2 [um]; hits", 100, -300, 300 ); + // + // width profiles: + // + h422 = fs->make("h422", "PXB2 sxy vs phi; phi2 [deg]; sxy [um]", 360, -180, 180, 0, 99 ); + h423 = fs->make("h423", "PXB2 sz vs phi; phi2 [deg]; sz [um]", 360, -180, 180, 0, 199 ); + + h424 = fs->make("h424", "PXB2 sxy vs z; z2 [cm]; sxy [um]", 80, -20, 20, 0, 99 ); + h425 = fs->make("h425", "PXB2 sz vs z; z2 [cm]; sz [um]", 80, -20, 20, 0, 199 ); + + h426 = fs->make("h426", "PXB2 sxy vs pt; log(pt [GeV]); sxy [um]", 20, 0, 2, 0, 99 ); + h427 = fs->make("h427", "PXB2 sz vs pt; log(pt [GeV]); sz [um]", 20, 0, 2, 0, 199 ); + + h428 = fs->make("h428", "PXB2 sxy vs dip; dip [deg]; sxy [um]", 70, -70, 70, 0, 99 ); + h429 = fs->make("h429", "PXB2 sz vs dip; dip [deg]; sz [um]", 70, -70, 70, 0, 199 ); + + h430 = fs->make("h430", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150 ); + h431 = fs->make("h431", "residuals PXB2 dz2; dz2 [um]; hits", 100, -300, 300 ); + + h432 = fs->make("h432", "residuals PXB2 dz2, 18 < |dip| < 50; dz2 [um]; hits", 100, -300, 300 ); + h433 = fs->make("h433", "PXB2 sz vs pt, best dip; log(pt [GeV]); sz [um]", 20, 0, 2, 0, 199 ); + + h434 = fs->make("h434", "PXB2 sxy vs inc; phi inc 2 [deg]; sxy [um]", 40, -10, 10, 0, 99 ); + + h435 = fs->make("h435", "ukap-kap; dkap; tracks", 100, -0.01, 0.01 ); + h436 = fs->make("h436", "uphi-phi; dphi; tracks", 100, -0.1, 0.1 ); + h437 = fs->make("h437", "udca-dca; ddca; tracks", 100, -0.1, 0.1 ); + + h440 = fs->make("h440", "pixel track dcap, pt > 2; dcap [um]; hits", 100, -1000, 1000 ); + h441 = fs->make("h441", "pixel track dcap, pt > 4; dcap [um]; hits", 100, -1000, 1000 ); + + h442 = fs->make("h442", "pixel track dcap vs phi; phi2 [deg]; [um]", 180, -180, 180, -500, 499 ); + h443 = fs->make("h443", "pixel tracks dcap vs pt; log(pt [GeV]); [um]", 20, 0, 2, -500, 499 ); + + h444 = fs->make("h444", "pixel track sdcap vs phi; phi2 [deg]; sdcap [um]", 180, -180, 180, 0, 499 ); + h445 = fs->make("h445", "pixel tracks sdcap vs pt; log(pt [GeV]); sdcap [um]", 20, 0, 2, 0, 499 ); + + h450 = fs->make("h450", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150 ); + h451 = fs->make("h451", "residuals PXB2 dz2; dz2 [um]; hits", 100, -300, 300 ); + +} +// +//---------------------------------------------------------------------- +// method called for each event: +// +void Triplet::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){ + + using namespace std; + using namespace edm; + using namespace reco; + using namespace math; + + const double pi = 4*atan(1); + const double wt = 180/pi; + const double twopi = 2*pi; + const double pihalf = 2*atan(1); + //const double sqrtpihalf = sqrt(pihalf); + + myCounters::neve++; + + if( myCounters::prevrun != iEvent.run() ){ + + time_t unixZeit = iEvent.time().unixTime(); + + cout << "new run " << iEvent.run(); + cout << ", LumiBlock " << iEvent.luminosityBlock(); + cout << " taken " << ctime(&unixZeit); // ctime has endline + + myCounters::prevrun = iEvent.run(); + + }// new run + + int idbg = 0; + if( myCounters::neve < 1 ) idbg = 1; + // + //-------------------------------------------------------------------- + // beam spot: + // + edm::Handle rbs; + iEvent.getByLabel( "offlineBeamSpot", rbs ); + + XYZPoint bsP = XYZPoint(0,0,0); + //int ibs = 0; + + if( !rbs.failedToGet() && rbs.isValid() ){ + + //ibs = 1; + h000->Fill( rbs->betaStar() ); + h001->Fill( rbs->emittanceX() ); + h002->Fill( rbs->BeamWidthX()*1e4 ); + h003->Fill( rbs->BeamWidthY()*1e4 ); + h004->Fill( rbs->sigmaZ() ); + h005->Fill( rbs->x0() ); + h006->Fill( rbs->y0() ); + h007->Fill( rbs->z0() ); + h008->Fill( rbs->dxdz()*1e3 ); + h009->Fill( rbs->dydz()*1e3 ); + bsP = XYZPoint( rbs->x0(), rbs->y0(), rbs->z0() ); + + if( idbg ){ + cout << "beam spot x " << rbs->x0(); + cout << ", y " << rbs->y0(); + cout << ", z " << rbs->z0(); + cout << endl; + } + + }// bs valid + // + //-------------------------------------------------------------------- + // primary vertices: + // + Handle vertices; + iEvent.getByLabel( "offlinePrimaryVertices", vertices ); + + if( vertices.failedToGet() ) return; + if( !vertices.isValid() ) return; + + h010->Fill( vertices->size() ); + + // need vertex global point for tracks + // from #include "DataFormats/GeometryVector/interface/GlobalPoint.h" + // Global points are three-dimensional by default + // typedef Global3DPoint GlobalPoint; + // typedef Point3DBase< float, GlobalTag> Global3DPoint; + + XYZPoint vtxN = XYZPoint(0,0,0); + XYZPoint vtxP = XYZPoint(0,0,0); + + double bestNdof = 0; + double maxSumPt = 0; + Vertex bestPvx; + + for( VertexCollection::const_iterator iVertex = vertices->begin(); + iVertex != vertices->end(); ++iVertex ) { + + if( ! iVertex->isValid() ) h011->Fill( iVertex->z() ); + else{ + if( iVertex->isFake() ) + h012->Fill( iVertex->z() ); + + else{ + h013->Fill( iVertex->z() ); + h014->Fill( iVertex->x() ); + h015->Fill( iVertex->y() ); + h016->Fill( iVertex->ndof() ); + h017->Fill( iVertex->ndof() ); + + if( idbg ){ + cout << "vertex"; + cout << ": x " << iVertex->x(); + cout << ", y " << iVertex->y(); + cout << ", z " << iVertex->z(); + cout << ", ndof " << iVertex->ndof(); + cout << ", sumpt " << iVertex->p4().pt(); + cout << endl; + } + + if( iVertex->hasRefittedTracks() ) + h018->Fill( iVertex->z() ); // empty in Zmumu sample + else + h019->Fill( iVertex->z() ); // all here: why? + + if( iVertex->ndof() > bestNdof ) { + bestNdof = iVertex->ndof(); + vtxN = XYZPoint( iVertex->x(), iVertex->y(), iVertex->z() ); + } + + h021->Fill( iVertex->p4().pt() ); + + if( iVertex->p4().pt() > maxSumPt ) { + maxSumPt = iVertex->p4().pt(); + vtxP = XYZPoint( iVertex->x(), iVertex->y(), iVertex->z() ); + bestPvx = *iVertex; + } + }// non-fake + }//valid + } // loop over vertices + + h022->Fill( maxSumPt ); + +#ifndef SIM + if( maxSumPt < 1 ) return; +#endif + if( maxSumPt < 1 ) vtxP = vtxN; + + /* + if( ibs ) { + vtxP.SetX( bsP.x() ); // beam spot. should take tilt into account! + vtxP.SetY( bsP.y() ); // beam spot. should take tilt into account! + } + */ + + h023->Fill( vtxP.x() ); + h024->Fill( vtxP.y() ); + h025->Fill( vtxP.z() ); + + //double xbs = 0; + //double ybs = 0; + //if( ibs ) { + //xbs = bsP.x(); + //ybs = bsP.y(); + //} + //else { + //xbs = vtxP.x(); + //ybs = vtxP.y(); + //} + // + //-------------------------------------------------------------------- + // MET: + // +// edm::Handle< edm::View > pfMEThandle; +// iEvent.getByLabel("pfMet", pfMEThandle); + +// if( !pfMEThandle.failedToGet() && pfMEThandle.isValid()){ + +// h026->Fill( pfMEThandle->front().sumEt() ); +// h027->Fill( pfMEThandle->front().et() ); + +// } + // + //-------------------------------------------------------------------- + // tracks: + // + Handle tracks; + + iEvent.getByLabel( "generalTracks", tracks ); + + if( tracks.failedToGet() ) return; + if( !tracks.isValid() ) return; + + h040->Fill( tracks->size() ); + // + // get tracker geometry: + // + edm::ESHandle pDD; + iSetup.get().get(pDD); + + if( !pDD.isValid() ) { + cout << "Unable to find TrackerDigiGeometry. Return\n"; + return; + } + // + // loop over tracker detectors: + // + for( TrackerGeometry::DetContainer::const_iterator idet = pDD->dets().begin(); + idet != pDD->dets().end(); ++idet ) { + + DetId mydetId = (*idet)->geographicalId(); + uint32_t mysubDet = mydetId.subdetId(); + + if( idbg ){ + + cout << "Det " << mydetId.det(); + cout << ", subDet " << mydetId.subdetId(); + + if( mysubDet == PixelSubdetector::PixelBarrel ) { + + cout << ": PXB layer " << PXBDetId(mydetId).layer(); + cout << ", ladder " << PXBDetId(mydetId).ladder(); + cout << ", module " << PXBDetId(mydetId).module(); + cout << ", at R " << (*idet)->position().perp(); + cout << ", F " << (*idet)->position().barePhi()*wt; + cout << ", z " << (*idet)->position().z(); + cout << endl; + cout << "rot x"; + cout << "\t" << (*idet)->rotation().xx(); + cout << "\t" << (*idet)->rotation().xy(); + cout << "\t" << (*idet)->rotation().xz(); + cout << endl; + cout << "rot y"; + cout << "\t" << (*idet)->rotation().yx(); + cout << "\t" << (*idet)->rotation().yy(); + cout << "\t" << (*idet)->rotation().yz(); + cout << endl; + cout << "rot z"; + cout << "\t" << (*idet)->rotation().zx(); + cout << "\t" << (*idet)->rotation().zy(); + cout << "\t" << (*idet)->rotation().zz(); + cout << endl; + // + // normal vector: includes alignment (varies from module to module along z on one ladder) + // neighbouring ladders alternate with inward/outward orientation + // + cout << "normal"; + cout << ": x " << (*idet)->surface().normalVector().x(); + cout << ", y " << (*idet)->surface().normalVector().y(); + cout << ", z " << (*idet)->surface().normalVector().z(); + cout << ", f " << (*idet)->surface().normalVector().barePhi()*wt; + + }//PXB + + if( mysubDet == PixelSubdetector::PixelEndcap ) { + + cout << ": PXD side " << PXFDetId(mydetId).side(); + cout << ", disk " << PXFDetId(mydetId).disk(); + cout << ", blade " << PXFDetId(mydetId).blade(); + cout << ", panel " << PXFDetId(mydetId).panel(); + cout << ", module " << PXFDetId(mydetId).module(); + cout << ", at R " << (*idet)->position().perp(); + cout << ", F " << (*idet)->position().barePhi()*wt; + cout << ", z " << (*idet)->position().z(); + cout << endl; + cout << "rot x"; + cout << "\t" << (*idet)->rotation().xx(); + cout << "\t" << (*idet)->rotation().xy(); + cout << "\t" << (*idet)->rotation().xz(); + cout << endl; + cout << "rot y"; + cout << "\t" << (*idet)->rotation().yx(); + cout << "\t" << (*idet)->rotation().yy(); + cout << "\t" << (*idet)->rotation().yz(); + cout << endl; + cout << "rot z"; + cout << "\t" << (*idet)->rotation().zx(); + cout << "\t" << (*idet)->rotation().zy(); + cout << "\t" << (*idet)->rotation().zz(); + cout << endl; + cout << "normal"; + cout << ": x " << (*idet)->surface().normalVector().x(); + cout << ", y " << (*idet)->surface().normalVector().y(); + cout << ", z " << (*idet)->surface().normalVector().z(); + cout << ", f " << (*idet)->surface().normalVector().barePhi()*wt; + + }//PXD + + cout << endl; + + }//idbg + + }//idet + + // + // transient track builder, needs B-field from data base (global tag in .py) + // + edm::ESHandle theB; + + iSetup.get().get( "TransientTrackBuilder", theB ); + // + // transient rec hits: + // + ESHandle hitBuilder; + iSetup.get().get( "WithTrackAngle", hitBuilder ); + // + // + // + double sumpt = 0; + double sumq = 0; + int kk = -1; + Surface::GlobalPoint origin = Surface::GlobalPoint(0,0,0); + // + //---------------------------------------------------------------------------- + // Tracks: + // + for( TrackCollection::const_iterator iTrack = tracks->begin(); + iTrack != tracks->end(); ++iTrack ) { + + kk++; + + // cpt = cqRB = 0.3*R[m]*B[T] = 1.14*R[m] for B=3.8T + // D = 2R = 2*pt/1.14 + // calo: D = 1.3 m => pt = 0.74 GeV/c + + double pt = iTrack->pt(); + + if( pt < 0.75 ) continue;// curls up + //if( pt < 1.75 ) continue;// want sharper image + + if( abs( iTrack->dxy(vtxP) ) > 5*iTrack->dxyError() ) continue; // not prompt + + double logpt = log(pt) / log(10); + double charge = iTrack->charge(); + h041->Fill( charge ); + h042->Fill( pt ); + h043->Fill( pt ); + + if( idbg ) { + cout << "Track " << kk; + cout << ": pt " << iTrack->pt(); + cout << ", eta " << iTrack->eta(); + cout << ", phi " << iTrack->phi()*wt; + cout << ", dxyv " << iTrack->dxy(vtxP); + cout << ", dzv " << iTrack->dz(vtxP); + cout << endl; + } + + const reco::HitPattern& hp = iTrack->hitPattern(); + + h045->Fill( hp.numberOfValidTrackerHits() ); + h046->Fill( hp.numberOfValidPixelBarrelHits() ); + h047->Fill( hp.trackerLayersWithMeasurement() ); + h048->Fill( hp.pixelBarrelLayersWithMeasurement() ); + + double phi = iTrack->phi(); + double dca = iTrack->d0(); // w.r.t. origin + //double dca = -iTrack->dxy(); // dxy = -d0 + double dip = iTrack->lambda(); + double z0 = iTrack->dz(); + double tet = pihalf - dip; + //double eta = iTrack->eta(); + // + // transient track: + // + TransientTrack tTrack = theB->build(*iTrack); + + double kap = tTrack.initialFreeState().transverseCurvature(); + + TrajectoryStateClosestToPoint tscp = tTrack.trajectoryStateClosestToPoint( origin ); + + //cout<Fill( tscp.referencePoint().x() ); // 0.0 + h058->Fill( tscp.referencePoint().y() ); // 0.0 + h059->Fill( tscp.referencePoint().z() ); // 0.0 + kap = tscp.perigeeParameters().transverseCurvature(); + phi = tscp.perigeeParameters().phi(); + dca = tscp.perigeeParameters().transverseImpactParameter(); + tet = tscp.perigeeParameters().theta(); + z0 = tscp.perigeeParameters().longitudinalImpactParameter(); + dip = pihalf - tet; + + h051->Fill( kap - tTrack.initialFreeState().transverseCurvature() ); + h052->Fill( phi - iTrack->phi() ); + h053->Fill( dca - iTrack->d0() ); + h054->Fill( dip - iTrack->lambda() ); + h055->Fill( z0 - iTrack->dz() ); + + } + + double tmp = abs(kap * pt); + h056->Fill( tmp ); + double rho1 = pt/0.0114257; + if(charge>0) rho1 = -rho1; + double kap1 = 1./rho1; + double tmp1 = (kap1-kap); + h049->Fill( tmp1 ); + + //cout<= abs(dca) ) { + // + // sin(delta phi): + // + double sindp = ( 0.5*kap * (R1*R1 + dca*dca) - dca ) / (R1*erd); + fpos1 = phi + asin(sindp); // phi position at R1 + // + // sin(alpha): + // + double sina = R1*kap * sqrt( 1.0 - sindp*sindp ); + // + // s = alpha / kappa: + // + if( sina >= 1.0 ) + s1 = pi / kap; + else{ + if( sina <= -1.0 ) + s1 = -pi / kap; + else + s1 = asin(sina) / kap;//always positive + } + // + // Check direction: limit to half-turn + // + if( hkk * ( R1*R1 - dca*dca ) > erd ) s1 = pi/abs(kap) - s1; // always positive + + }// R1 > dca + + if( fpos1 > pi ) fpos1 -= twopi; + else if( fpos1 < -pi ) fpos1 += twopi; + + double zR1 = z0 + s1*tan(dip); // z at R1 + // + //-------------------------------------------------------------------------- + // loop over tracker detectors: + // + double xcrss[99]; + double ycrss[99]; + double zcrss[99]; + int ncrss = 0; + + for( TrackerGeometry::DetContainer::const_iterator idet = pDD->dets().begin(); + idet != pDD->dets().end(); ++idet ) { + + DetId mydetId = (*idet)->geographicalId(); + uint32_t mysubDet = mydetId.subdetId(); + + if( mysubDet != PixelSubdetector::PixelBarrel ) continue; + /* + cout << ": PXB layer " << PXBDetId(mydetId).layer(); + cout << ", ladder " << PXBDetId(mydetId).ladder(); + cout << ", module " << PXBDetId(mydetId).module(); + cout << ", at R1 " << (*idet)->position().perp(); + cout << ", F " << (*idet)->position().barePhi()*wt; + cout << ", z " << (*idet)->position().z(); + cout << endl; + */ + + if( PXBDetId(mydetId).layer() == 1 ) { + + double dz = zR1 - (*idet)->position().z(); + + if( abs(dz) > 4.0 ) continue; + + double df = fpos1 - (*idet)->position().barePhi(); + + if( df > pi ) df -= twopi; + else if( df < -pi ) df += twopi; + + if( abs(df)*wt > 36.0 ) continue; + // + // normal vector: includes alignment (varies from module to module along z on one ladder) + // neighbouring ladders alternate with inward/outward orientation + // + /* + cout << "normal"; + cout << ": x " << (*idet)->surface().normalVector().x(); + cout << ", y " << (*idet)->surface().normalVector().y(); + cout << ", z " << (*idet)->surface().normalVector().z(); + cout << ", f " << (*idet)->surface().normalVector().barePhi()*wt; + cout << endl; + */ + + double phiN = (*idet)->surface().normalVector().barePhi();//normal vector + + double phidet = phiN - pihalf;// orientation of sensor plane in x-y + double ux = cos(phidet);// vector in sensor plane + double uy = sin(phidet); + double x = (*idet)->position().x(); + double y = (*idet)->position().y(); + // + // intersect helix with line: FUNRXY (in FUNPHI) from V. Blobel + // factor f for intersect point (x + f*ux, y + f*uy) + // + double a = kap * ( ux*ux + uy*uy ) * 0.5; + double b = erd * ( ux*sf - uy*cf ) + kap * ( ux*x + uy*y ); + double c = drd + erd * ( x*sf - y*cf ) + kap * ( x*x + y*y ) * 0.5; + double dis = b*b - 4.0*a*c; + double f = 0; + + if( dis > 0 ) { + + dis = sqrt(dis); + double f1 = 0; + double f2 = 0; + + if( b < 0 ) { + f1 = 0.5 * ( dis - b ) / a; + f2 = 2.0 * c / ( dis - b ); + } + else{ + f1 = -0.5 * ( dis + b ) / a; + f2 = -2.0 * c / ( dis + b ); + } + + f = f1; + if( abs(f2) < abs(f1) ) f = f2; + + }//dis + + xcrss[ncrss] = x + f*ux; + ycrss[ncrss] = y + f*uy; + double r = sqrt( xcrss[ncrss]*xcrss[ncrss] + ycrss[ncrss]*ycrss[ncrss] ); + + double s = 0; + + if( r >= abs(dca) ) { + double sindp = ( 0.5*kap * ( r*r + dca*dca ) - dca ) / (r*erd); + double sina = r*kap * sqrt( 1.0 - sindp*sindp ); + if( sina >= 1.0 ) + s = pi / kap; + else{ + if( sina <= -1.0 ) + s = -pi / kap; + else + s = asin(sina) / kap; + } + if( hkk * ( r*r - dca*dca ) > erd ) s = pi/abs(kap) - s; + } + + zcrss[ncrss] = z0 + s*tan(dip); // z at r + + ncrss++; + + }//PXB1 + + }//idet + // + //-------------------------------------------------------------------------- + // rec hits from track extra: + // + if( iTrack->extra().isNonnull() &&iTrack->extra().isAvailable() ){ + + h044->Fill( tTrack.recHitsSize() ); // tTrack + + double x1 = 0; + double y1 = 0; + double z1 = 0; + double x2 = 0; + double y2 = 0; + double z2 = 0; + double x3 = 0; + double y3 = 0; + double z3 = 0; + int n1 = 0; + int n2 = 0; + int n3 = 0; + //double phiN2 = 0; + + for( trackingRecHit_iterator irecHit = iTrack->recHitsBegin(); + irecHit != iTrack->recHitsEnd(); ++irecHit){ + + if( (*irecHit)->isValid() ){ + + DetId detId = (*irecHit)->geographicalId(); + + // enum Detector { Tracker=1, Muon=2, Ecal=3, Hcal=4, Calo=5 }; + + if( detId.det() != 1 ){ + cout << "rec hit ID = " << detId.det() << " not in tracker!?!?\n"; + continue; + } + + uint32_t subDet = detId.subdetId(); + + // enum SubDetector{ PixelBarrel=1, PixelEndcap=2 }; + // enum SubDetector{ TIB=3, TID=4, TOB=5, TEC=6 }; + + h060->Fill( subDet ); + // + // build hit: (from what?) + // + TransientTrackingRecHit::RecHitPointer trecHit = hitBuilder->build( &*(*irecHit) ); + + double xloc = trecHit->localPosition().x();// 1st meas coord + double yloc = trecHit->localPosition().y();// 2nd meas coord or zero + //double zloc = trecHit->localPosition().z();// up, always zero + h061->Fill( xloc ); + h062->Fill( yloc ); + + double gX = trecHit->globalPosition().x(); + double gY = trecHit->globalPosition().y(); + double gZ = trecHit->globalPosition().z(); + double gF = atan2( gY, gX ); + double gR = sqrt( gX*gX + gY*gY ); + + h063->Fill( gR ); + h064->Fill( gF*wt ); + h065->Fill( gZ ); + + // GeomDet* igeomdet = trecHit->det(); + //double phiN = trecHit->det()->surface().normalVector().barePhi();//normal vector + + if( subDet == PixelSubdetector::PixelBarrel || + subDet == StripSubdetector::TIB || + subDet == StripSubdetector::TOB ) { // barrel + + h066->Fill( gX, gY ); + + }//barrel + + if( subDet == PixelSubdetector::PixelBarrel ) { + + int ilay = PXBDetId(detId).layer(); + int ilad = PXBDetId(detId).ladder(); + int imod = PXBDetId(detId).module(); + bool halfmod = 0; + + h100->Fill( ilay ); // 1,2,3 + + if( ilay == 1 ){ + + n1++; + x1 = gX; + y1 = gY; + z1 = gZ; + + h101->Fill( ilad );// 1..20 + h102->Fill( imod );// 1..8 + + h103->Fill( gR ); + h104->Fill( gF*wt ); + h105->Fill( gZ ); + + h106->Fill( gF*wt, gZ ); // phi-z of hit + + if( ilad == 5 ) halfmod = 1; + else if( ilad == 6 ) halfmod = 1; + else if( ilad == 15 ) halfmod = 1; + else if( ilad == 16 ) halfmod = 1; + + if( !halfmod ){ + h107->Fill( xloc, yloc ); // hit within one module + } + // + // my crossings: + // + for( int icrss = 0; icrss < ncrss; ++icrss ){ + + double fcrss = atan2( ycrss[icrss], xcrss[icrss] ); + double df = gF - fcrss; + if( df > pi ) df -= twopi; + else if( df < -pi ) df += twopi; + double du = gR*df; + double dz = gZ - zcrss[icrss]; + + if( abs(du) < 0.01 && abs(dz) < 0.05 ) h111->Fill( gX, gY ); + h112->Fill( du*1E4 ); + h113->Fill( dz*1E4 ); + + if( abs(dz) < 0.02 ) h114->Fill( du*1E4 ); + if( abs(du) < 0.01 ) h115->Fill( dz*1E4 ); + + }//crss + + }//PXB1 + + if( ilay == 2 ){ + + n2++; + x2 = gX; + y2 = gY; + z2 = gZ; + //phiN2 = phiN; + + h201->Fill( ilad );// 1..32 + h202->Fill( imod );//1..8 + + h203->Fill( gR ); + h204->Fill( gF*wt ); + h205->Fill( gZ ); + + h206->Fill( gF*wt, gZ ); // phi-z of hit + + if( ilad == 8 ) halfmod = 1; + else if( ilad == 9 ) halfmod = 1; + else if( ilad == 24 ) halfmod = 1; + else if( ilad == 25 ) halfmod = 1; + + if( !halfmod ){ + h207->Fill( xloc, yloc ); // hit within one module + } + + }//PXB2 + + if( ilay == 3 ){ + + n3++; + x3 = gX; + y3 = gY; + z3 = gZ; + + h301->Fill( ilad );//1..44 + h302->Fill( imod );//1..8 + + h303->Fill( gR ); + h304->Fill( gF*wt ); + h305->Fill( gZ ); + + h306->Fill( gF*wt, gZ ); // phi-z of hit + + if( ilad == 11 ) halfmod = 1; + if( ilad == 12 ) halfmod = 1; + if( ilad == 33 ) halfmod = 1; + if( ilad == 34 ) halfmod = 1; + + if( !halfmod ){ + h307->Fill( xloc, yloc ); // hit within one module + } + + }//PXB3 + + }//PXB + + }//valid + + }//loop rechits + // + // 1-2-3 triplet: + // + //if( hp.pixelBarrelLayersWithMeasurement() == 3 ){ + if( n1*n2*n3 > 0 ) { + + double dca2 = 0., dz2=0.; + //triplets(x1,y1,z1,x2,y2,z2,x3,y3,z3,kap,dca2,dz2); + double ptsig = pt; + if(charge<0.) ptsig = -pt; + triplets(x1,y1,z1,x2,y2,z2,x3,y3,z3,ptsig,dca2,dz2); + + if( pt > 4 ) { + h410->Fill( dca2*1E4 ); + h411->Fill( dz2*1E4 ); + } + if( pt > 12 ) { + h420->Fill( dca2*1E4 ); + h421->Fill( dz2*1E4 ); + if( hp.trackerLayersWithMeasurement() > 8 ) { + h430->Fill( dca2*1E4 ); + h431->Fill( dz2*1E4 ); + } + //if( phiinc*wt > -1 && phiinc*wt < 7 ){ + //h450->Fill( dca2*1E4 ); + //h451->Fill( dz2*1E4 ); + //} + } + + // + // residual profiles: alignment check + // + if( pt > 4 ) { + //h412->Fill( f2*wt, dca2*1E4 ); + //h413->Fill( f2*wt, dz2*1E4 ); + + h414->Fill( z2, dca2*1E4 ); + h415->Fill( z2, dz2*1E4 ); + } + + h416->Fill( logpt, dca2*1E4 ); + h417->Fill( logpt, dz2*1E4 ); + if( iTrack->charge() > 0 ) h418->Fill( logpt, dca2*1E4 ); + else h419->Fill( logpt, dca2*1E4 ); + + // profile of abs(dca) gives mean abs(dca): + // mean of abs(Gauss) = 0.7979 * RMS = 1/sqrt(pi/2) + // => rms = sqrt(pi/2) * mean of abs (sqrt(pi/2) = 1.2533) + // point resolution = 1/sqrt(3/2) * triplet middle residual width + // => sqrt(pi/2)*sqrt(2/3) = sqrt(pi/3) = 1.0233, almost one + + if( pt > 4 ) { + + //h422->Fill( f2*wt, abs(dca2)*1E4 ); + //h423->Fill( f2*wt, abs(dz2)*1E4 ); + + h424->Fill( z2, abs(dca2)*1E4 ); + h425->Fill( z2, abs(dz2)*1E4 ); + + //h428->Fill( udip*wt, abs(dca2)*1E4 ); + //h429->Fill( udip*wt, abs(dz2)*1E4 ); + //if( abs(dip)*wt > 18 && abs(dip)*wt < 50 ) h432->Fill( dz2*1E4 ); + + //h434->Fill( phiinc*wt, abs(dca2)*1E4 ); + + }//pt + + h426->Fill( logpt, abs(dca2)*1E4 ); + h427->Fill( logpt, abs(dz2)*1E4 ); + //if( abs(dip)*wt > 18 && abs(dip)*wt < 50 ) h433->Fill( logpt, abs(dz2)*1E4 ); + + }//3 PXB layers + + }//extra + + sumpt += iTrack->pt(); + sumq += iTrack->charge(); + + }// loop over tracks + + h028->Fill( sumpt ); + h029->Fill( sumq ); + +} + +void Triplet::triplets(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3, + double ptsig, double & dca2,double & dz2) { + using namespace std; + const double pi = 4*atan(1); + //const double wt = 180/pi; + const double twopi = 2*pi; + //const double pihalf = 2*atan(1); + //const double sqrtpihalf = sqrt(pihalf); + + double pt = abs(ptsig); + //double rho = pt/0.0114257; + double kap = 0.0114257/pt; + if(ptsig>0) kap = -kap; // kap i snegative for positive charge + + double rho = 1/kap; + double rinv = -kap; // Karimaki + + + //double f2 = atan2( y2, x2 );//position angle + + //h406->Fill( hp.numberOfValidTrackerHits() ); + //h407->Fill( hp.numberOfValidPixelBarrelHits() ); + //h408->Fill( hp.trackerLayersWithMeasurement() ); + + // Author: Johannes Gassner (15.11.1996) + // Make track from 2 space points and kappa (cmz98/ftn/csmktr.f) + // Definition of the Helix : + // + // x( t ) = X0 + KAPPA^-1 * SIN( PHI0 + t ) + // y( t ) = Y0 - KAPPA^-1 * COS( PHI0 + t ) t > 0 + // z( t ) = Z0 + KAPPA^-1 * TAN( DIP ) * t + + // Center of the helix in the xy-projection: + + // X0 = + ( DCA - KAPPA^-1 ) * SIN( PHI0 ) + // Y0 = - ( DCA - KAPPA^-1 ) * COS( PHI0 ) + // + // Point 1 must be in the inner layer, 3 in the outer: + // + double r1 = sqrt( x1*x1 + y1*y1 ); + double r3 = sqrt( x3*x3 + y3*y3 ); + + if( r3-r1 < 2.0 ) cout << "warn r1 = " << r1 << ", r3 = " << r3 << endl; + // + // Calculate the centre of the helix in xy-projection that + // transverses the two spacepoints. The points with the same + // distance from the two points are lying on a line. + // LAMBDA is the distance between the point in the middle of + // the two spacepoints and the centre of the helix. + // + // we already have kap and rho = 1/kap + // + + double lam = sqrt( -0.25 + + rho*rho / ( ( x1 - x3 )*( x1 - x3 ) + ( y1 - y3 )*( y1 - y3 ) ) ); + // + // There are two solutions, the sign of kap gives the information + // which of them is correct. + // + if( kap > 0 ) lam = -lam; + // + // ( X0, Y0 ) is the centre of the circle that describes the helix in xy-projection. + // + double x0 = 0.5*( x1 + x3 ) + lam * ( -y1 + y3 ); + double y0 = 0.5*( y1 + y3 ) + lam * ( x1 - x3 ); + // + // Calculate theta : + // + double num = ( y3 - y0 ) * ( x1 - x0 ) - ( x3 - x0 ) * ( y1 - y0 ); + double den = ( x1 - x0 ) * ( x3 - x0 ) + ( y1 - y0 ) * ( y3 - y0 ); + double tandip = kap * ( z3 - z1 ) / atan( num / den ); + //double udip = atan(tandip); + //double utet = pihalf - udip; + // + // To get phi0 in the right intervall one must differ two cases + // with positve and negative kap: + // + double uphi; + if( kap > 0 ) uphi = atan2( -x0, y0 ); + else uphi = atan2( x0, -y0 ); + // + // The distance of the closest approach DCA depends on the sign + // of kappa : + // + double udca; + if( kap > 0 ) udca = rho - sqrt( x0*x0 + y0*y0 ); + else udca = rho + sqrt( x0*x0 + y0*y0 ); + // + // angle from first hit to dca point: + // + double dphi = atan( ( ( x1 - x0 ) * y0 - ( y1 - y0 ) * x0 ) + / ( ( x1 - x0 ) * x0 + ( y1 - y0 ) * y0 ) ); + + double uz0 = z1 + tandip * dphi * rho; + + //h401->Fill( z2 ); + //h402->Fill( uphi - iTrack->phi() ); + //h403->Fill( udca - iTrack->d0() ); + //h404->Fill( udip - iTrack->lambda() ); + //h405->Fill( uz0 - iTrack->dz() ); + // + // interpolate to middle hit: + // cirmov + // we already have rinv = -kap + // + double cosphi = cos(uphi); + double sinphi = sin(uphi); + double dp = -x2*sinphi + y2*cosphi + udca; + double dl = -x2*cosphi - y2*sinphi; + double sa = 2*dp + rinv*(dp*dp+dl*dl); + dca2 = sa / ( 1 + sqrt(1 + rinv*sa) );// distance to hit 2 + + double ud = 1 + rinv*udca; + double phi2 = atan2( -rinv*x2 + ud*sinphi, rinv*y2 + ud*cosphi );//direction + + //double phiinc = phi2 - phiN2;//angle of incidence in rphi w.r.t. normal vector + // + // phiN alternates inward/outward + // reduce phiinc + //if( phiinc > pihalf ) phiinc -= pi; + //else if( phiinc < -pihalf ) phiinc += pi; + //h409->Fill( f2*wt, phiinc*wt ); + // + // arc length: + // + double xx = x2 + dca2 * sin(phi2); // point on track + double yy = y2 - dca2 * cos(phi2); + + double f0 = uphi;// + double kx = kap*xx; + double ky = kap*yy; + double kd = kap*udca; + // + // Solve track equation for s: + // + double dx = kx - (kd-1)*sin(f0); + double dy = ky + (kd-1)*cos(f0); + double ks = atan2( dx, -dy ) - f0;// turn angle + // + //--- Limit to half-turn: + // + if( ks > pi ) ks = ks - twopi; + else if( ks < -pi ) ks = ks + twopi; + + double s = ks*rho;// signed + double uz2 = uz0 + s*tandip; //track z at R2 + dz2 = z2 - uz2; + +} +//---------------------------------------------------------------------- +// method called just after ending the event loop: +// +void Triplet::endJob() { + + std::cout << "end of job after " << myCounters::neve << " events.\n"; + +} + +//define this as a plug-in +DEFINE_FWK_MODULE(Triplet); diff --git a/RecoLocalTracker/SiPixelClusterizer/test/readClusters_cfg.py b/RecoLocalTracker/SiPixelClusterizer/test/readClusters_cfg.py index d3de66780044e..5ad3efea0c8f0 100644 --- a/RecoLocalTracker/SiPixelClusterizer/test/readClusters_cfg.py +++ b/RecoLocalTracker/SiPixelClusterizer/test/readClusters_cfg.py @@ -32,7 +32,7 @@ process.hltPhysicsDeclared.L1GtReadoutRecordTag = 'gtDigis' process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(100) + input = cms.untracked.int32(-1) ) process.MessageLogger = cms.Service("MessageLogger", @@ -50,11 +50,18 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( # fill 3273 run 206940 - "/store/data/Run2012D/MinimumBias/RECO/PromptReco-v1/000/206/940/FA55823C-312C-E211-94AB-001D09F29533.root", +# "/store/data/Run2012D/MinimumBias/RECO/PromptReco-v1/000/206/940/FA55823C-312C-E211-94AB-001D09F29533.root", +# for MC +# 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100_pre10/clus/clus1.root' +# 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/clus/clus1.root' +# 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/rechits/rechits1.root' + 'file:../../../../../CMSSW_7_0_0_pre8/src/EventFilter/SiPixelRawToDigi/test/digis.root' + ) ) -process.source.lumisToProcess = cms.untracked.VLuminosityBlockRange('206940:0-206940:1027') +# for data +#process.source.lumisToProcess = cms.untracked.VLuminosityBlockRange('206940:0-206940:1027') process.TFileService = cms.Service("TFileService", fileName = cms.string('histo.root') @@ -63,25 +70,22 @@ process.load("Configuration.StandardSequences.Geometry_cff") process.load("Configuration.StandardSequences.MagneticField_38T_cff") -# what is this? -# process.load("Configuration.StandardSequences.Services_cff") -# what is this? -#process.load("SimTracker.Configuration.SimTracker_cff") -# needed for global transformation -# process.load("Configuration.StandardSequences.FakeConditions_cff") - process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: # 2012 -process.GlobalTag.globaltag = 'GR_P_V40::All' +#process.GlobalTag.globaltag = 'GR_P_V40::All' +# MC 2913 +process.GlobalTag.globaltag = 'MC_70_V1::All' process.analysis = cms.EDAnalyzer("ReadPixClusters", Verbosity = cms.untracked.bool(True), src = cms.InputTag("siPixelClusters"), ) -process.p = cms.Path(process.hltPhysicsDeclared*process.hltfilter*process.analysis) +# for data +#process.p = cms.Path(process.hltPhysicsDeclared*process.hltfilter*process.analysis) #process.p = cms.Path(process.hltPhysicsDeclared*process.analysis) -#process.p = cms.Path(process.analysis) +# for MC +process.p = cms.Path(process.analysis) # define an EndPath to analyze all other path results diff --git a/RecoLocalTracker/SiPixelClusterizer/test/testTracks_cfg.py b/RecoLocalTracker/SiPixelClusterizer/test/testTracks_cfg.py new file mode 100644 index 0000000000000..51c5c460d300d --- /dev/null +++ b/RecoLocalTracker/SiPixelClusterizer/test/testTracks_cfg.py @@ -0,0 +1,98 @@ +# +# Last update: new version for python +# +# +import FWCore.ParameterSet.Config as cms + +process = cms.Process("T") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(-1) +) + +process.MessageLogger = cms.Service("MessageLogger", +# debugModules = cms.untracked.vstring('TestPixTracks'), + destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), + cout = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') + threshold = cms.untracked.string('ERROR') + ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +) + +import HLTrigger.HLTfilters.hltHighLevel_cfi as hlt +# accept if 'path_1' succeeds +process.hltfilter = hlt.hltHighLevel.clone( +# Min-Bias +# HLTPaths = ['HLT_Physics_v*'], +# HLTPaths = ['HLT_Random_v*'], + HLTPaths = ['HLT_ZeroBias*'], +# HLTPaths = ['HLT_L1Tech54_ZeroBias*'], +# HLTPaths = ['HLT_L1Tech53_MB*'], +# Commissioning: +# HLTPaths = ['HLT_L1_Interbunch_BSC_v*'], +# HLTPaths = ['HLT_L1_PreCollisions_v*'], +# HLTPaths = ['HLT_BeamGas_BSC_v*'], +# HLTPaths = ['HLT_BeamGas_HF_v*'], +# HLTPaths = ['p*'], +# HLTPaths = ['path_?'], + andOr = True, # False = and, True=or + throw = False + ) + + +# to select PhysicsBit +process.load('HLTrigger.special.hltPhysicsDeclared_cfi') +process.hltPhysicsDeclared.L1GtReadoutRecordTag = 'gtDigis' + + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( +# 'file:../../../SimTracker/SiPixelDigitizer/test/tracks.root' + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/tracks/tracks1.root' + ) +) + + +process.TFileService = cms.Service("TFileService", + fileName = cms.string('histo_tracks.root') +) + +process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load("Configuration.StandardSequences.MagneticField_38T_cff") + +# needed for global transformation +# process.load("Configuration.StandardSequences.FakeConditions_cff") + +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: +process.GlobalTag.globaltag = "MC_70_V1::All" +# 2012 +#process.GlobalTag.globaltag = "GR_P_V40::All" +#process.GlobalTag.globaltag = "GR_P_V28::All" # A&B +# 2011 +# process.GlobalTag.globaltag = "GR_P_V20::All" +# process.GlobalTag.globaltag = "GR_R_311_V2::All" +# process.GlobalTag.globaltag = "GR_R_310_V2::All" +# for 2010 +# process.GlobalTag.globaltag = 'GR10_P_V5::All' +# process.GlobalTag.globaltag = 'GR10_P_V4::All' +# OK for 2009 LHC data +# process.GlobalTag.globaltag = 'CRAFT09_R_V4::All' + +process.d = cms.EDAnalyzer("TestWithTracks", + Verbosity = cms.untracked.bool(False), + src = cms.InputTag("generalTracks"), +# PrimaryVertexLabel = cms.untracked.InputTag("offlinePrimaryVertices"), +# trajectoryInput = cms.string("TrackRefitterP5") +# trajectoryInput = cms.string('cosmictrackfinderP5') +) + +#process.p = cms.Path(process.hltPhysicsDeclared*process.hltfilter*process.d) +# process.p = cms.Path(process.hltPhysicsDeclared*process.d) +process.p = cms.Path(process.d) + + + diff --git a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h index 64714832a20c5..74cf33c6be2c2 100644 --- a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h +++ b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h @@ -41,7 +41,6 @@ #include - class RectangularPixelTopology; class MagneticField; class PixelCPEBase : public PixelClusterParameterEstimator @@ -52,12 +51,16 @@ class PixelCPEBase : public PixelClusterParameterEstimator Param() : bz(-9e10f) {} float bz; // local Bz LocalVector drift; + float widthLAFraction; // Width-LA to Offset-LA }; public: PixelCPEBase(edm::ParameterSet const& conf, const MagneticField * mag = 0, - const SiPixelLorentzAngle * lorentzAngle = 0, const SiPixelCPEGenericErrorParm * genErrorParm = 0, - const SiPixelTemplateDBObject * templateDBobject = 0); + const SiPixelLorentzAngle * lorentzAngle = 0, + const SiPixelCPEGenericErrorParm * genErrorParm = 0, + const SiPixelTemplateDBObject * templateDBobject = 0, + const SiPixelLorentzAngle * lorentzAngleWidth = 0 + ); //-------------------------------------------------------------------------- @@ -74,6 +77,7 @@ class PixelCPEBase : public PixelClusterParameterEstimator const GeomDetUnit & det ) const { nRecHitsTotal_++ ; + //std::cout<<" in PixelCPEBase:localParameters(all) - "< cpe_; edm::ParameterSet pset_; + bool useLAWidthFromDB_; + bool useLAAlignmentOffsets_; }; diff --git a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h index 82ef9ec85f2ad..e02a644bb789d 100644 --- a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h +++ b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h @@ -78,6 +78,8 @@ class PixelCPETemplateReco : public PixelCPEBase mutable bool DoCosmics_; + mutable bool DoLorentz_; + mutable bool LoadTemplatesFromDB_; }; diff --git a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h index 82c5238af4f38..6e3abb210c8f7 100644 --- a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h +++ b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h @@ -15,6 +15,7 @@ class PixelCPETemplateRecoESProducer: public edm::ESProducer{ private: boost::shared_ptr cpe_; edm::ParameterSet pset_; + bool DoLorentz_; }; diff --git a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPEGenericESProducer.cc b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPEGenericESProducer.cc index 9b7cc87bf9916..46fd03c907948 100644 --- a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPEGenericESProducer.cc +++ b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPEGenericESProducer.cc @@ -20,8 +20,20 @@ using namespace edm; PixelCPEGenericESProducer::PixelCPEGenericESProducer(const edm::ParameterSet & p) { std::string myname = p.getParameter("ComponentName"); + // Use LA-width from DB. If both (upper and this) are false LA-width is calcuated from LA-offset + //useLAWidthFromDB_ = p.getParameter("useLAWidthFromDB"); + useLAWidthFromDB_ = p.existsAs("useLAWidthFromDB")? + p.getParameter("useLAWidthFromDB"):false; + // Use Alignment LA-offset + //useLAAlignmentOffsets_ = p.getParameter("useLAAlignmentOffsets"); + useLAAlignmentOffsets_ = p.existsAs("useLAAlignmentOffsets")? + p.getParameter("useLAAlignmentOffsets"):false; + pset_ = p; setWhatProduced(this,myname); + + //std::cout<<" ESProducer "< pDD; iRecord.getRecord().get( pDD ); + // Lorant angle for offsets ESHandle lorentzAngle; - iRecord.getRecord().get(lorentzAngle ); - - ESHandle genErrorParm; - iRecord.getRecord().get(genErrorParm); + if(useLAAlignmentOffsets_) // LA offsets from alignment + iRecord.getRecord().get("fromAlignment",lorentzAngle ); + else // standard LA, from calibration, label="" + iRecord.getRecord().get(lorentzAngle ); + + // add the new la width object + ESHandle lorentzAngleWidth; + const SiPixelLorentzAngle * lorentzAngleWidthProduct = 0; + if(useLAWidthFromDB_) { // use the width LA + iRecord.getRecord().get("forWidth",lorentzAngleWidth ); + lorentzAngleWidthProduct = lorentzAngleWidth.product(); + } else { lorentzAngleWidthProduct = NULL;} // do not use it + //std::cout<<" la width "< genErrorParm; + iRecord.getRecord().get(genErrorParm); + + // errors come from this, replace by a lighter object + ESHandle templateDBobject; + iRecord.getRecord().get(templateDBobject); + - ESHandle templateDBobject; - iRecord.getRecord().get(templateDBobject); + cpe_ = boost::shared_ptr(new PixelCPEGeneric(pset_,magfield.product(),lorentzAngle.product(),genErrorParm.product(),templateDBobject.product(),lorentzAngleWidthProduct) ); - cpe_ = boost::shared_ptr(new PixelCPEGeneric(pset_,magfield.product(),lorentzAngle.product(),genErrorParm.product(),templateDBobject.product()) ); - //ToDo? Replace blah.product() with ESHandle + //ToDo? Replace blah.product() with ESHandle return cpe_; } diff --git a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc index 9319a36e10e2c..bc1312291f0e9 100644 --- a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc +++ b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc @@ -20,8 +20,15 @@ using namespace edm; PixelCPETemplateRecoESProducer::PixelCPETemplateRecoESProducer(const edm::ParameterSet & p) { std::string myname = p.getParameter("ComponentName"); + + //DoLorentz_ = p.getParameter("DoLorentz"); // True when LA from alignment is used + DoLorentz_ = p.existsAs("DoLorentz")?p.getParameter("DoLorentz"):false; + pset_ = p; setWhatProduced(this,myname); + + //std::cout<<" from ES Producer Templates "< pDD; iRecord.getRecord().get( pDD ); - edm::ESHandle lorentzAngle; - iRecord.getRecord().get(lorentzAngle); - - ESHandle templateDBobject; - iRecord.getRecord().get(templateDBobject); - - cpe_ = boost::shared_ptr(new PixelCPETemplateReco(pset_,magfield.product(),lorentzAngle.product(),templateDBobject.product() )); + edm::ESHandle lorentzAngle; + const SiPixelLorentzAngle * lorentzAngleProduct = 0; + if(DoLorentz_) { // LA correction from alignment + iRecord.getRecord().get("fromAlignment",lorentzAngle); + lorentzAngleProduct = lorentzAngle.product(); + } else { // Normal, deafult LA actually is NOT needed + //iRecord.getRecord().get(lorentzAngle); + lorentzAngleProduct=NULL; // null is ok becuse LA is not use by templates in this mode + } + + ESHandle templateDBobject; + iRecord.getRecord().get(templateDBobject); + + // cpe_ = boost::shared_ptr(new PixelCPETemplateReco(pset_,magfield.product(),lorentzAngle.product(),templateDBobject.product() )); + cpe_ = boost::shared_ptr(new PixelCPETemplateReco(pset_,magfield.product(),lorentzAngleProduct,templateDBobject.product() )); return cpe_; } diff --git a/RecoLocalTracker/SiPixelRecHits/python/PixelCPEGeneric_cfi.py b/RecoLocalTracker/SiPixelRecHits/python/PixelCPEGeneric_cfi.py index 61b2962406569..afa6ba500f905 100644 --- a/RecoLocalTracker/SiPixelRecHits/python/PixelCPEGeneric_cfi.py +++ b/RecoLocalTracker/SiPixelRecHits/python/PixelCPEGeneric_cfi.py @@ -3,7 +3,6 @@ PixelCPEGenericESProducer = cms.ESProducer("PixelCPEGenericESProducer", ComponentName = cms.string('PixelCPEGeneric'), - TanLorentzAnglePerTesla = cms.double(0.106), Alpha2Order = cms.bool(True), PixelErrorParametrization = cms.string('NOTcmsim'), @@ -40,8 +39,21 @@ LoadTemplatesFromDB = cms.bool(True), # petar, for clusterProbability() from TTRHs - ClusterProbComputationFlag = cms.int32(0) - + ClusterProbComputationFlag = cms.int32(0), + + # new parameters added in 1/14, dk + # LA defined by hand, FOR TESTING ONLY, not for production + # 0.0 means that the offset is taken from DB + #lAOffset = cms..double(0.0), + #lAWidthBPix = cms.double(0.0), + #lAWidthFPix = cms.double(0.0), + + # Flag to select the source of LA-Width + # Normal = True, use LA from DB, not ready yet + useLAWidthFromDB = cms.bool(False), + # if lAWith=0 and useLAWidthFromDB=false than width is calculated from lAOffset. + # Use the LA-Offsets from Alignment instead of our calibration + useLAAlignmentOffsets = cms.bool(False), ) diff --git a/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py b/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py index 7c91b1cf6b8a1..a3e543472930c 100644 --- a/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py +++ b/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py @@ -2,7 +2,6 @@ templates = cms.ESProducer("PixelCPETemplateRecoESProducer", ComponentName = cms.string('PixelCPETemplateReco'), - #TanLorentzAnglePerTesla = cms.double(0.106), speed = cms.int32(-2), #PixelErrorParametrization = cms.string('NOTcmsim'), Alpha2Order = cms.bool(True), @@ -10,9 +9,12 @@ # petar, for clusterProbability() from TTRHs ClusterProbComputationFlag = cms.int32(0), - # gavril - DoCosmics = cms.bool(False), - + DoCosmics = cms.bool(False), + # The flag to regulate if the LA offset is taken from Alignment + # Will be True in the future for offline RECO. kevin + DoLorentz = cms.bool(False), + LoadTemplatesFromDB = cms.bool(True) + ) diff --git a/RecoLocalTracker/SiPixelRecHits/python/SiPixelRecHits_cfi.py b/RecoLocalTracker/SiPixelRecHits/python/SiPixelRecHits_cfi.py index ae12a26b45dfe..ba16ed203f569 100644 --- a/RecoLocalTracker/SiPixelRecHits/python/SiPixelRecHits_cfi.py +++ b/RecoLocalTracker/SiPixelRecHits/python/SiPixelRecHits_cfi.py @@ -2,12 +2,8 @@ siPixelRecHits = cms.EDProducer("SiPixelRecHitConverter", src = cms.InputTag("siPixelClusters"), - # untracked string ClusterCollLabel = "siPixelClusters" CPE = cms.string('PixelCPEGeneric'), VerboseLevel = cms.untracked.int32(0), - #TanLorentzAnglePerTesla = cms.double(0.106), - #Alpha2Order = cms.bool(True), - #speed = cms.int32(0) ) diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEBase.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEBase.cc index 8bee8db3fd7bd..8a383d11c8366 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEBase.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEBase.cc @@ -18,7 +18,6 @@ #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h" -//#define TPDEBUG #define CORRECT_FOR_BIG_PIXELS // MessageLogger @@ -46,18 +45,24 @@ namespace { // A fairly boring constructor. All quantities are DetUnit-dependent, and // will be initialized in setTheDet(). //----------------------------------------------------------------------------- -PixelCPEBase::PixelCPEBase(edm::ParameterSet const & conf, const MagneticField *mag, const SiPixelLorentzAngle * lorentzAngle, - const SiPixelCPEGenericErrorParm * genErrorParm, const SiPixelTemplateDBObject * templateDBobject) +PixelCPEBase::PixelCPEBase(edm::ParameterSet const & conf, const MagneticField *mag, + const SiPixelLorentzAngle * lorentzAngle, + const SiPixelCPEGenericErrorParm * genErrorParm, + const SiPixelTemplateDBObject * templateDBobject, + const SiPixelLorentzAngle * lorentzAngleWidth) : theDet(nullptr), theTopol(nullptr), theRecTopol(nullptr), theParam(nullptr), nRecHitsTotal_(0), nRecHitsUsedEdge_(0), probabilityX_(0.0), probabilityY_(0.0), probabilityQ_(0.0), qBin_(0), isOnEdge_(false), hasBadPixels_(false), spansTwoROCs_(false), hasFilledProb_(false), + useLAAlignmentOffsets_(false), useLAOffsetFromConfig_(false), + useLAWidthFromConfig_(false), useLAWidthFromDB_(false), loc_trk_pred_(0.0, 0.0, 0.0, 0.0) { //--- Lorentz angle tangent per Tesla lorentzAngle_ = lorentzAngle; + lorentzAngleWidth_ = lorentzAngleWidth; //--- Algorithm's verbosity theVerboseLevel = @@ -83,6 +88,14 @@ PixelCPEBase::PixelCPEBase(edm::ParameterSet const & conf, const MagneticField * // clusterProbComputationFlag_ = (unsigned int) conf.getParameter("ClusterProbComputationFlag"); + + // For safety initilaize the parameters which are used by generic algo only to 0 + lAOffset_ = 0.0; + lAWidthBPix_ = 0.0; + lAWidthFPix_ = 0.0; + + LogDebug("PixelCPEBase") <<" LA constants - " + <( &det ); @@ -107,19 +123,23 @@ PixelCPEBase::setTheDet( const GeomDetUnit & det, const SiPixelCluster & cluster //--- theDet->type() returns a GeomDetType, which implements subDetector() thePart = theDet->type().subDetector(); -#ifdef EDM_ML_DEBUG + //cout<<" in PixelCPEBase:setTheDet - in det "<surface(), but @@ -155,36 +175,31 @@ PixelCPEBase::setTheDet( const GeomDetUnit & det, const SiPixelCluster & cluster theSign = isFlipped() ? -1 : 1; + widthLAFraction_=1.0; // preset the LA fraction to 1. - // will cache if not yest there (need some of the above) + // will cache if not yet there (need some of the above) theParam = ¶m(); - // this "has wrong sign..." driftDirection_ = (*theParam).drift; - - - //--- The Lorentz shift. - theLShiftX = lorentzShiftX(); - - theLShiftY = lorentzShiftY(); + widthLAFraction_ = (*theParam).widthLAFraction; + + //cout<<" in PixelCPEBase:setTheDet - "<index(); if (i>=int(m_Params.size())) m_Params.resize(i+1); // should never happen! @@ -392,44 +410,11 @@ PixelCPEBase::Param const & PixelCPEBase::param() const { LocalVector Bfield = theDet->surface().toLocal(magfield_->inTesla(theDet->surface().position())); p.drift = driftDirection(Bfield ); p.bz = Bfield.z(); + p.widthLAFraction = widthLAFraction_; } return p; } - -//----------------------------------------------------------------------------- -// HALF OF the Lorentz shift (so for the full shift multiply by 2), and -// in the units of pitch. (So note these are neither local nor measurement -// units!) -//----------------------------------------------------------------------------- -float PixelCPEBase::lorentzShiftX() const -{ - LocalVector dir = getDrift(); - - // max shift in cm - float xdrift = dir.x()/dir.z() * theThickness; - // express the shift in units of pitch, - // divide by 2 to get the average correction - float lshift = xdrift / (thePitchX*2.); - - return lshift; - - -} - -float PixelCPEBase::lorentzShiftY() const -{ - - LocalVector dir = getDrift(); - - float ydrift = dir.y()/dir.z() * theThickness; - float lshift = ydrift / (thePitchY * 2.f); - return lshift; - - -} - - //----------------------------------------------------------------------------- // Drift direction. // Works OK for barrel and forward. @@ -437,8 +422,6 @@ float PixelCPEBase::lorentzShiftY() const // used in the digitizer (SiPixelDigitizerAlgorithm.cc). // Assumption: setTheDet() has been called already. // -// Petar (2/23/07): uhm, actually, there is a bug in the sign for both X and Y! -// (The signs have been fixed in SiPixelDigitizer, but not in here.) //----------------------------------------------------------------------------- LocalVector PixelCPEBase::driftDirection( GlobalVector bfield ) const { @@ -451,13 +434,43 @@ PixelCPEBase::driftDirection( GlobalVector bfield ) const { LocalVector PixelCPEBase::driftDirection( LocalVector Bfield ) const { - - - auto langle = lorentzAngle_->getLorentzAngle(theDet->geographicalId().rawId()); + const bool LocalPrint = false; + + //auto langle = lorentzAngle_->getLorentzAngle(theDet->geographicalId().rawId()); + // Use LA from DB or from config + float langle = 0.; + if( !useLAOffsetFromConfig_ ) { // get it from DB + if(lorentzAngle_ != NULL) { // a real LA object + langle = lorentzAngle_->getLorentzAngle(theDet->geographicalId().rawId()); + } else { // no LA, unused + //cout<<" LA object is NULL, assume LA = 0"<getLorentzAngle(theDet->geographicalId().rawId()); + if(langleWidth!=0.0) widthLAFraction_ = std::abs(langleWidth/langle); + else widthLAFraction_ = 1.0; + if(LocalPrint) cout<<" Will use LA Width from DB "< 0) @@ -59,17 +63,43 @@ PixelCPEGeneric::PixelCPEGeneric(edm::ParameterSet const & conf, DoCosmics_ = conf.getParameter("DoCosmics"); LoadTemplatesFromDB_ = conf.getParameter("LoadTemplatesFromDB"); + // This LA related parameters are only relevant for the Generic algo + + //lAOffset_ = conf.getUntrackedParameter("lAOffset",0.0); + lAOffset_ = conf.existsAs("lAOffset")? + conf.getParameter("lAOffset"):0.0; + //lAWidthBPix_ = conf.getUntrackedParameter("lAWidthBPix",0.0); + lAWidthBPix_ = conf.existsAs("lAWidthBPix")? + conf.getParameter("lAWidthBPix"):0.0; + //lAWidthFPix_ = conf.getUntrackedParameter("lAWidthFPix",0.0); + lAWidthFPix_ = conf.existsAs("lAWidthFPix")? + conf.getParameter("lAWidthFPix"):0.0; + + // Use LA-offset from config, for testing only + if(lAOffset_>0.0) useLAOffsetFromConfig_ = true; + // Use LA-width from config, split into fpix & bpix, for testing only + if(lAWidthBPix_>0.0 || lAWidthFPix_>0.0) useLAWidthFromConfig_ = true; + + // Use LA-width from DB. If both (upper and this) are false LA-width is calcuated from LA-offset + //useLAWidthFromDB_ = conf.getParameter("useLAWidthFromDB"); + useLAWidthFromDB_ = conf.existsAs("useLAWidthFromDB")? + conf.getParameter("useLAWidthFromDB"):false; + + // Use Alignment LA-offset + useLAAlignmentOffsets_ = conf.existsAs("useLAAlignmentOffsets")? + conf.getParameter("useLAAlignmentOffsets"):false; + + if ( !UseErrorsFromTemplates_ && ( TruncatePixelCharge_ || IrradiationBiasCorrection_ || DoCosmics_ || - LoadTemplatesFromDB_ ) ) - { - throw cms::Exception("PixelCPEGeneric::PixelCPEGeneric: ") - << "\nERROR: UseErrorsFromTemplates_ is set to False in PixelCPEGeneric_cfi.py. " - << " In this case it does not make sense to set any of the following to True: " - << " TruncatePixelCharge_, IrradiationBiasCorrection_, DoCosmics_, LoadTemplatesFromDB_ !!!" - << "\n\n"; - } + LoadTemplatesFromDB_ ) ) { + throw cms::Exception("PixelCPEGeneric::PixelCPEGeneric: ") + << "\nERROR: UseErrorsFromTemplates_ is set to False in PixelCPEGeneric_cfi.py. " + << " In this case it does not make sense to set any of the following to True: " + << " TruncatePixelCharge_, IrradiationBiasCorrection_, DoCosmics_, LoadTemplatesFromDB_ !!!" + << "\n\n"; + } if ( UseErrorsFromTemplates_ ) { @@ -100,8 +130,8 @@ PixelCPEGeneric::PixelCPEGeneric(edm::ParameterSet const & conf, //cout << "(int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl; //cout << endl; - //yes, these should be config parameters! - //default case... + // Default case for rechit errors in case other, more correct, errors are not vailable + // This are constants. Maybe there is a more efficienct way to store them. xerr_barrel_l1_= {0.00115, 0.00120, 0.00088}; xerr_barrel_l1_def_=0.01030; yerr_barrel_l1_= {0.00375,0.00230,0.00250,0.00250,0.00230,0.00230,0.00210,0.00210,0.00240}; @@ -155,7 +185,12 @@ PixelCPEGeneric::PixelCPEGeneric(edm::ParameterSet const & conf, LocalPoint PixelCPEGeneric::localPosition(const SiPixelCluster& cluster) const { + + //cout<<" in PixelCPEGeneric:localPosition - "<getTemplateID(theDet->geographicalId().rawId()); @@ -283,11 +318,12 @@ PixelCPEGeneric::localPosition(const SiPixelCluster& cluster) const cout << "\t >>> Generic:: processing X" << endl; #endif + float chargeWidth = (lorentzShiftInCmX_ * widthLAFraction_); float xPos = generic_position_formula( cluster.sizeX(), Q_f_X, Q_l_X, local_URcorn_LLpix.x(), local_LLcorn_URpix.x(), - 0.5*lorentzShiftInCmX_, // 0.5 * lorentz shift in + chargeWidth, // lorentz shift in cm cotalpha_, thePitchX, theRecTopol->isItBigPixelInX( cluster.minPixelRow() ), @@ -297,33 +333,42 @@ PixelCPEGeneric::localPosition(const SiPixelCluster& cluster) const the_size_cutX); // cut for eff charge width &&& - #ifdef EDM_ML_DEBUG + // apply the lorentz offset correction + xPos = xPos + (0.5 * lorentzShiftInCmX_); + +#ifdef EDM_ML_DEBUG if (theVerboseLevel > 20) cout << "\t >>> Generic:: processing Y" << endl; #endif + chargeWidth = (lorentzShiftInCmY_ * widthLAFraction_); float yPos = generic_position_formula( cluster.sizeY(), Q_f_Y, Q_l_Y, local_URcorn_LLpix.y(), local_LLcorn_URpix.y(), - 0.5*lorentzShiftInCmY_, // 0.5 * lorentz shift in cm + chargeWidth, // lorentz shift in cm cotbeta_, - thePitchY, // 0.5 * lorentz shift (may be 0) + thePitchY, theRecTopol->isItBigPixelInY( cluster.minPixelCol() ), theRecTopol->isItBigPixelInY( cluster.maxPixelCol() ), the_eff_charge_cut_lowY, the_eff_charge_cut_highY, the_size_cutY); // cut for eff charge width &&& - - // Apply irradiation corrections. + + // apply the lorentz offset correction + yPos = yPos + (0.5 * lorentzShiftInCmY_); + + // Apply irradiation corrections. NOT USED FOR NOW if ( IrradiationBiasCorrection_ ) { if ( cluster.sizeX() == 1 ) { + // ggiurgiu@jhu.edu, 02/03/09 : for size = 1, the Lorentz shift is already accounted by the irradiation correction + xPos = xPos - (0.5 * lorentzShiftInCmX_); + // Find if pixel is double (big). - bool bigInX = theRecTopol->isItBigPixelInX( cluster.maxPixelRow() ); - + bool bigInX = theRecTopol->isItBigPixelInX( cluster.maxPixelRow() ); if ( !bigInX ) { //cout << "Apply correction dx1 = " << dx1 << " to xPos = " << xPos << endl; @@ -344,9 +389,11 @@ PixelCPEGeneric::localPosition(const SiPixelCluster& cluster) const if ( cluster.sizeY() == 1 ) { + // ggiurgiu@jhu.edu, 02/03/09 : for size = 1, the Lorentz shift is already accounted by the irradiation correction + yPos = yPos - (0.5 * lorentzShiftInCmY_); + // Find if pixel is double (big). - bool bigInY = theRecTopol->isItBigPixelInY( cluster.maxPixelCol() ); - + bool bigInY = theRecTopol->isItBigPixelInY( cluster.maxPixelCol() ); if ( !bigInY ) { //cout << "Apply correction dy1 = " << dy1 << " to yPos = " << yPos << endl; @@ -366,6 +413,8 @@ PixelCPEGeneric::localPosition(const SiPixelCluster& cluster) const } // if ( IrradiationBiasCorrection_ ) + //cout<<" in PixelCPEGeneric:localPosition - pos = "<containsBigPixelInX( minPixelRow, maxPixelRow ); bool bigInY = theRecTopol->containsBigPixelInY( minPixelCol, maxPixelCol ); - if unlikely( isUpgrade_ ||(!with_track_angle && DoCosmics_) ) + if unlikely( isUpgrade_ ||(!with_track_angle && DoCosmics_) ) { //cout << "Track angles are not known and we are processing cosmics." << endl; //cout << "Default angle estimation which assumes track from PV (0,0,0) does not work." << endl; diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPETemplateReco.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPETemplateReco.cc index fd88b736f7217..ec2f6975a6f68 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/PixelCPETemplateReco.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/PixelCPETemplateReco.cc @@ -46,7 +46,7 @@ const int cluster_matrix_size_y = 21; PixelCPETemplateReco::PixelCPETemplateReco(edm::ParameterSet const & conf, const MagneticField * mag, const SiPixelLorentzAngle * lorentzAngle, const SiPixelTemplateDBObject * templateDBobject) - : PixelCPEBase(conf, mag, lorentzAngle, 0, templateDBobject) + : PixelCPEBase(conf, mag, lorentzAngle, 0, templateDBobject, 0) { //cout << endl; //cout << "Constructing PixelCPETemplateReco::PixelCPETemplateReco(...)................................................." << endl; @@ -57,7 +57,11 @@ PixelCPETemplateReco::PixelCPETemplateReco(edm::ParameterSet const & conf, //-- Use Magnetic field at (0,0,0) to select a template ID [Morris, 6/25/08] (temporary until we implement DB access) DoCosmics_ = conf.getParameter("DoCosmics"); + //DoLorentz_ = conf.getParameter("DoLorentz"); // True when LA from alignment is used + DoLorentz_ = conf.existsAs("DoLorentz")?conf.getParameter("DoLorentz"):false; + LoadTemplatesFromDB_ = conf.getParameter("LoadTemplatesFromDB"); + //cout << " PixelCPETemplateReco : (int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl; //cout << "field_magnitude = " << field_magnitude << endl; @@ -119,6 +123,11 @@ PixelCPETemplateReco::localPosition(const SiPixelCluster& cluster) const else fpix = true; // yes, it's forward + // Compute the Lorentz shifts for this detector element for the Alignment Group + if ( DoLorentz_ ) computeLorentzShifts(); + + + int ID = -9999; if ( LoadTemplatesFromDB_ ) @@ -444,6 +453,23 @@ PixelCPETemplateReco::localPosition(const SiPixelCluster& cluster) const // go back to the module coordinate system templXrec_ += lp.x(); templYrec_ += lp.y(); + + // Compute the Alignment Group Corrections [template ID should already be selected from call to reco procedure] + if ( DoLorentz_ ) { + // Donly if the lotentzshift has meaningfull numbers + if( lorentzShiftInCmX_!= 0.0 || lorentzShiftInCmY_!= 0.0 ) { + // the LA width/shift returned by templates use (+) + // the LA width/shift produced by PixelCPEBase for positive LA is (-) + // correct this by iserting (-) + float templateLorwidthCmX = -micronsToCm*templ_.lorxwidth(); + float templateLorwidthCmY = -micronsToCm*templ_.lorywidth(); + // now, correctly, we can use the difference of shifts + templXrec_ += 0.5*(lorentzShiftInCmX_ - templateLorwidthCmX); + templYrec_ += 0.5*(lorentzShiftInCmY_ - templateLorwidthCmY); + //cout << "Templates: la lorentz offset = " <<(0.5*(lorentzShiftInCmX_-templateLorwidthCmX))<< endl; //dk + } //else {cout<<" LA is 0, disable offset corrections "< + + + + + + +# + + +# +# + + +# +# + + +# + + + +# for tracks +# +# + + + + + + + + diff --git a/RecoLocalTracker/SiPixelRecHits/test/CPEAccessTester.cc b/RecoLocalTracker/SiPixelRecHits/test/CPEAccessTester.cc index 611a19030e83d..99e7829ed3a9a 100644 --- a/RecoLocalTracker/SiPixelRecHits/test/CPEAccessTester.cc +++ b/RecoLocalTracker/SiPixelRecHits/test/CPEAccessTester.cc @@ -2,21 +2,23 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/EDAnalyzer.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" + #include "FWCore/Framework/interface/Event.h" -#include "DataFormats/Common/interface/Handle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/InputTag.h" -#include "FWCore/Framework/interface/ESHandle.h" + +#include "DataFormats/Common/interface/Handle.h" #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h" - #include #include @@ -45,4 +47,5 @@ class CPEAccessTester : public edm::EDAnalyzer { private: edm::ParameterSet conf_; }; - +//define this as a plug-in +DEFINE_FWK_MODULE(CPEAccessTester); diff --git a/RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.cc b/RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.cc index add3dd5d7a0be..217502b1e26a2 100644 --- a/RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.cc +++ b/RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.cc @@ -10,11 +10,11 @@ #include #include -#include "RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.h" - -#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" + //#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" @@ -26,6 +26,9 @@ #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "Geometry/CommonTopologies/interface/PixelTopology.h" +#include "RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" + #include "DataFormats/DetId/interface/DetId.h" #include "DataFormats/SiPixelDetId/interface/PXBDetId.h" #include "DataFormats/SiPixelDetId/interface/PXFDetId.h" @@ -174,116 +177,6 @@ void ReadPixelRecHit::beginJob() { hdetsPerLay2F = fs->make( "hdetsPerLay2F", "Full dets per layer l2", 257, -0.5, 256.5); - -// hpixid = new TH1F( "hpixid", "Pix det id", 10, 0., 10.); -// hpixsubid = new TH1F( "hpixsubid", "Pix Barrel id", 10, 0., 10.); -// hlayerid = new TH1F( "hlayerid", "Pix layer id", 10, 0., 10.); -// hladder1id = new TH1F( "hladder1id", "Ladder L1 id", 50, 0., 50.); -// hladder2id = new TH1F( "hladder2id", "Ladder L2 id", 50, 0., 50.); -// hladder3id = new TH1F( "hladder3id", "Ladder L3 id", 50, 0., 50.); -// hz1id = new TH1F( "hz1id", "Z-index id L1", 10, 0., 10.); -// hz2id = new TH1F( "hz2id", "Z-index id L2", 10, 0., 10.); -// hz3id = new TH1F( "hz3id", "Z-index id L3", 10, 0., 10.); - -// hrecHitsPerDet1 = new TH1F( "hrecHitsPerDet1", "RecHits per det l1", -// 200, -0.5, 199.5); -// hrecHitsPerDet2 = new TH1F( "hrecHitsPerDet2", "RecHits per det l2", -// 200, -0.5, 199.5); -// hrecHitsPerDet3 = new TH1F( "hrecHitsPerDet3", "RecHits per det l3", -// 200, -0.5, 199.5); -// hrecHitsPerLay1 = new TH1F( "hrecHitsPerLay1", "RecHits per layer l1", -// 2000, -0.5, 1999.5); -// hrecHitsPerLay2 = new TH1F( "hrecHitsPerLay2", "RecHits per layer l2", -// 2000, -0.5, 1999.5); - - -// hrecHitsPerLay3 = new TH1F( "hrecHitsPerLay3", "RecHits per layer l3", -// 2000, -0.5, 1999.5); -// hdetsPerLay1 = new TH1F( "hdetsPerLay1", "Full dets per layer l1", -// 161, -0.5, 160.5); -// hdetsPerLay3 = new TH1F( "hdetsPerLay3", "Full dets per layer l3", -// 353, -0.5, 352.5); -// hdetsPerLay2 = new TH1F( "hdetsPerLay2", "Full dets per layer l2", -// 257, -0.5, 256.5); - -// hcharge1 = new TH1F( "hcharge1", "Clu charge l1", 200, 0.,200.); //in ke -// hcharge2 = new TH1F( "hcharge2", "Clu charge l2", 200, 0.,200.); -// hcharge3 = new TH1F( "hcharge3", "Clu charge l3", 200, 0.,200.); -// hadcCharge1 = new TH1F( "hadcCharge1", "pix charge l1", 50, 0.,50.); //in ke -// hadcCharge2 = new TH1F( "hadcCharge2", "pix charge l2", 50, 0.,50.); //in ke -// hadcCharge3 = new TH1F( "hadcCharge3", "pix charge l3", 50, 0.,50.); //in ke -// hadcCharge1big = new TH1F( "hadcCharge1big", "big pix charge l1", 50, 0.,50.); //in ke - -// hxpos1 = new TH1F( "hxpos1", "Layer 1 cols", 700,-3.5,3.5); -// hxpos2 = new TH1F( "hxpos2", "Layer 2 cols", 700,-3.5,3.5); -// hxpos3 = new TH1F( "hxpos3", "Layer 3 cols", 700,-3.5,3.5); - -// hypos1 = new TH1F( "hypos1", "Layer 1 rows", 200,-1.,1.); -// hypos2 = new TH1F( "hypos2", "Layer 2 rows", 200,-1.,1.); -// hypos3 = new TH1F( "hypos3", "layer 3 rows", 200,-1.,1.); - -// hsize1 = new TH1F( "hsize1", "layer 1 clu size",100,-0.5,99.5); -// hsize2 = new TH1F( "hsize2", "layer 2 clu size",100,-0.5,99.5); -// hsize3 = new TH1F( "hsize3", "layer 3 clu size",100,-0.5,99.5); -// hsizex1 = new TH1F( "hsizex1", "lay1 clu size in x", -// 10,-0.5,9.5); -// hsizex2 = new TH1F( "hsizex2", "lay2 clu size in x", -// 10,-0.5,9.5); -// hsizex3 = new TH1F( "hsizex3", "lay3 clu size in x", -// 10,-0.5,9.5); -// hsizey1 = new TH1F( "hsizey1", "lay1 clu size in y", -// 20,-0.5,19.5); -// hsizey2 = new TH1F( "hsizey2", "lay2 clu size in y", -// 20,-0.5,19.5); -// hsizey3 = new TH1F( "hsizey3", "lay3 clu size in y", -// 20,-0.5,19.5); - -// hdetr = new TH1F("hdetr","det r",150,0.,15.); -// hdetz = new TH1F("hdetz","det z",520,-26.,26.); - -// // Forward edcaps -// hdetrF = new TH1F("hdetrF","Fdet r",150,5.,20.); -// hdetzF = new TH1F("hdetzF","Fdet z",600,-60.,60.); - -// hdisk = new TH1F( "hdisk", "FPix disk id", 10, 0., 10.); -// hblade = new TH1F( "hblade", "FPix blade id", 30, 0., 30.); -// hmodule = new TH1F( "hmodule", "FPix plaq. id", 10, 0., 10.); -// hpanel = new TH1F( "hpanel", "FPix panel id", 10, 0., 10.); -// hside = new TH1F( "hside", "FPix size id", 10, 0., 10.); - -// hcharge1F = new TH1F( "hcharge1F", "Clu charge 21", 200, 0.,200.); //in ke -// hcharge2F = new TH1F( "hcharge2F", "Clu charge 22", 200, 0.,200.); -// hxpos1F = new TH1F( "hxpos1F", "Disk 1 cols", 700,-3.5,3.5); -// hxpos2F = new TH1F( "hxpos2F", "Disk 2 cols", 700,-3.5,3.5); -// hypos1F = new TH1F( "hypos1F", "Disk 1 rows", 200,-1.,1.); -// hypos2F = new TH1F( "hypos2F", "Disk 2 rows", 200,-1.,1.); -// hsize1F = new TH1F( "hsize1F", "Disk 1 clu size",100,-0.5,99.5); -// hsize2F = new TH1F( "hsize2F", "Disk 2 clu size",100,-0.5,99.5); -// hsizex1F = new TH1F( "hsizex1F", "d1 clu size in x", -// 10,-0.5,9.5); -// hsizex2F = new TH1F( "hsizex2F", "d2 clu size in x", -// 10,-0.5,9.5); -// hsizey1F = new TH1F( "hsizey1F", "d1 clu size in y", -// 20,-0.5,19.5); -// hsizey2F = new TH1F( "hsizey2F", "d2 clu size in y", -// 20,-0.5,19.5); -// hadcCharge1F = new TH1F( "hadcCharge1F", "pix charge d1", 50, 0.,50.); //in ke -// hadcCharge2F = new TH1F( "hadcCharge2F", "pix charge d2", 50, 0.,50.); //in ke - -// hrecHitsPerDet1F = new TH1F( "hrecHitsPerDet1F", "RecHits per det l1", -// 200, -0.5, 199.5); -// hrecHitsPerDet2F = new TH1F( "hrecHitsPerDet2F", "RecHits per det l2", -// 200, -0.5, 199.5); -// hrecHitsPerLay1F = new TH1F( "hrecHitsPerLay1F", "RecHits per layer l1", -// 2000, -0.5, 1999.5); -// hrecHitsPerLay2F = new TH1F( "hrecHitsPerLay2F", "RecHits per layer l2", -// 2000, -0.5, 1999.5); -// hdetsPerLay1F = new TH1F( "hdetsPerLay1F", "Full dets per layer l1", -// 161, -0.5, 160.5); -// hdetsPerLay2F = new TH1F( "hdetsPerLay2F", "Full dets per layer l2", -// 257, -0.5, 256.5); - - cout<<" book histos "<ls(); - //hFile->pwd(); - //hFile->Write(); - //hFile->Close(); -#endif - } //--------------------------------------------------------------------- // Functions that gets called by framework every event void ReadPixelRecHit::analyze(const edm::Event& e, const edm::EventSetup& es) { using namespace edm; - const bool localPrint = false; - //const bool localPrint = true; + //const bool localPrint = false; + const bool localPrint = true; // Get event setup (to get global transformation) edm::ESHandle geom; @@ -367,24 +253,24 @@ void ReadPixelRecHit::analyze(const edm::Event& e, double detZ = theGeomDet->surface().position().z(); double detR = theGeomDet->surface().position().perp(); - const BoundPlane& plane = theGeomDet->surface(); //for transf. - double detThick = theGeomDet->specificSurface().bounds().thickness(); + //const BoundPlane& plane = theGeomDet->surface(); //for transf. unused + //double detThick = theGeomDet->specificSurface().bounds().thickness(); unused //const RectangularPixelTopology * topol = //dynamic_cast(&(theGeomDet->specificTopology())); const PixelTopology * topol = &(theGeomDet->specificTopology()); - int cols = theGeomDet->specificTopology().ncolumns(); - int rows = theGeomDet->specificTopology().nrows(); + //int cols = theGeomDet->specificTopology().ncolumns(); UNUSED + //int rows = theGeomDet->specificTopology().nrows(); unsigned int layer=0, disk=0, ladder=0, zindex=0, blade=0, panel=0; if(subid==1) { // Subdet it, pix barrel=1 ++numberOfDetUnits; PXBDetId pdetId = PXBDetId(detId); - unsigned int detTypeP=pdetId.det(); - unsigned int subidP=pdetId.subdetId(); + //unsigned int detTypeP=pdetId.det(); unused + //unsigned int subidP=pdetId.subdetId(); unused // Barell layer = 1,2,3 layer=pdetId.layer(); // Barrel ladder id 1-20,32,44. @@ -459,16 +345,20 @@ void ReadPixelRecHit::analyze(const edm::Event& e, //----Loop over rechits for this detId SiPixelRecHitCollection::DetSet::const_iterator pixeliter=detset.begin(); SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end(); - for(;pixeliter!=rechitRangeIteratorEnd;++pixeliter){//loop on the rechit - if(print) cout <<" Position " << pixeliter->localPosition() << endl; + for(;pixeliter!=rechitRangeIteratorEnd;++pixeliter) { //loop on the rechit + + if(print) cout <<" No Position " << endl; numOfRecHits++; - LocalPoint lp = pixeliter->localPosition(); - LocalError le = pixeliter->localPositionError(); - float xRecHit = lp.x(); - float yRecHit = lp.y(); - if(localPrint) cout<<" RecHit: "<localPosition(); + //LocalError le = pixeliter->localPositionError(); UNUSED + //float xRecHit = lp.x(); + //float yRecHit = lp.y(); + //if(localPrint) cout<<" RecHit: "<measurementPosition(xRecHit,yRecHit); //GlobalPoint GP = PixGeom->surface().toGlobal(Local3DPoint(lp)); @@ -505,24 +395,25 @@ void ReadPixelRecHit::analyze(const edm::Event& e, //if(localPrint) cout<<" Pixels in this cluster "< > chanMap; // Channel map // Look at pixels in this cluster. ADC is calibrated, in electrons - for (int i = 0; i < pixelsVec.size(); ++i) { + for (unsigned int i = 0; i < pixelsVec.size(); ++i) { float pixx = pixelsVec[i].x; // index as a float so = i+0.5 float pixy = pixelsVec[i].y; float adc = ((pixelsVec[i].adc)/1000); // in kelec. - int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy)); + // OLD way + //int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy)); //if(RectangularPixelTopology::isItBigPixelInX(int(pixx))) bigInX=true; //if(RectangularPixelTopology::isItBigPixelInY(int(pixy))) bigInY=true; - bool bigInX = (PixelTopology::isItBigPixelInX(int(pixx))); - bool bigInY = (PixelTopology::isItBigPixelInY(int(pixy))); + //bool bigInX = (PixelTopology::isItBigPixelInX(int(pixx))); + //bool bigInY = (PixelTopology::isItBigPixelInY(int(pixy))); bool edgeInX = topol->isItEdgePixelInX(int(pixx)); bool edgeInY = topol->isItEdgePixelInY(int(pixy)); if(localPrint) - cout<Fill(adc); - if(bigInX || bigInY) hadcCharge1big->Fill(adc); + //if(bigInX || bigInY) hadcCharge1big->Fill(adc); } else if(layer==2) { hadcCharge2->Fill(adc); } else if(layer==3) { @@ -548,8 +439,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, #ifdef DO_HISTO if(layer==1) { hcharge1->Fill(ch); - hxpos1->Fill(yRecHit); - hypos1->Fill(xRecHit); + //hxpos1->Fill(yRecHit); + //hypos1->Fill(xRecHit); hsize1->Fill(float(size)); hsizex1->Fill(float(sizeX)); hsizey1->Fill(float(sizeY)); @@ -558,8 +449,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, } else if(layer==2) { // layer 2 hcharge2->Fill(ch); - hxpos2->Fill(yRecHit); - hypos2->Fill(xRecHit); + //hxpos2->Fill(yRecHit); + //hypos2->Fill(xRecHit); hsize2->Fill(float(size)); hsizex2->Fill(float(sizeX)); hsizey2->Fill(float(sizeY)); @@ -568,8 +459,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, } else if(layer==3) { // Layer 3 hcharge3->Fill(ch); - hxpos3->Fill(yRecHit); - hypos3->Fill(xRecHit); + //hxpos3->Fill(yRecHit); + //hypos3->Fill(xRecHit); hsize3->Fill(float(size)); hsizex3->Fill(float(sizeX)); hsizey3->Fill(float(sizeY)); @@ -578,8 +469,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, } else if(disk==1) { hcharge1F->Fill(ch); - hxpos1F->Fill(yRecHit); - hypos1F->Fill(xRecHit); + //hxpos1F->Fill(yRecHit); + //hypos1F->Fill(xRecHit); hsize1F->Fill(float(size)); hsizex1F->Fill(float(sizeX)); hsizey1F->Fill(float(sizeY)); @@ -589,8 +480,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, } else if(disk==2) { // disk 2 hcharge2F->Fill(ch); - hxpos2F->Fill(yRecHit); - hypos2F->Fill(xRecHit); + //hxpos2F->Fill(yRecHit); + //hypos2F->Fill(xRecHit); hsize2F->Fill(float(size)); hsizex2F->Fill(float(sizeX)); hsizey2F->Fill(float(sizeY)); @@ -634,5 +525,8 @@ void ReadPixelRecHit::analyze(const edm::Event& e, hdetsPerLay2F ->Fill(float(numberOfDetUnits2F)); hdetsPerLay3 ->Fill(float(numberOfDetUnits3)); #endif + } +//define this as a plug-in +DEFINE_FWK_MODULE(ReadPixelRecHit); diff --git a/RecoLocalTracker/SiPixelRecHits/test/SealModules.cc b/RecoLocalTracker/SiPixelRecHits/test/SealModules.cc deleted file mode 100644 index d59078891ae09..0000000000000 --- a/RecoLocalTracker/SiPixelRecHits/test/SealModules.cc +++ /dev/null @@ -1,12 +0,0 @@ - -//#include "PluginManager/ModuleDef.h" - -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.h" -#include "RecoLocalTracker/SiPixelRecHits/test/CPEAccessTester.cc" - - - -DEFINE_FWK_MODULE(ReadPixelRecHit); -DEFINE_FWK_MODULE(CPEAccessTester); diff --git a/RecoLocalTracker/SiPixelRecHits/test/SiPixelDBTest_cfg.py b/RecoLocalTracker/SiPixelRecHits/test/SiPixelDBTest_cfg.py index 446916cab6542..5e05a1286c058 100644 --- a/RecoLocalTracker/SiPixelRecHits/test/SiPixelDBTest_cfg.py +++ b/RecoLocalTracker/SiPixelRecHits/test/SiPixelDBTest_cfg.py @@ -1,10 +1,19 @@ import FWCore.ParameterSet.Config as cms process = cms.Process("SiPixelDBTestt") + +# does not run - fake conditions do not exist process.load("Configuration.StandardSequences.FakeConditions_cff") -useFakeSource = True -#useFakeSource = False +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: +# 2012 +#process.GlobalTag.globaltag = 'GR_P_V40::All' +# MC 2913 +#process.GlobalTag.globaltag = 'MC_70_V1::All' + + +#useFakeSource = True +useFakeSource = False useCPEGeneric = True #useCPEGeneric = False @@ -31,7 +40,8 @@ if useCPEGeneric: process.load("CalibTracker.SiPixelESProducers.SiPixelFakeTemplateDBObjectESSource_cfi") process.PoolDBESSource = cms.ESSource("PoolDBESSource", - process.CondDBSetup,loadAll = cms.bool(True), + process.CondDBSetup, + loadAll = cms.bool(True), toGet = cms.VPSet(cms.PSet( record = cms.string('SiPixelCPEGenericErrorParmRcd'), tag = cms.string('SiPixelCPEGenericErrorParm') @@ -56,12 +66,10 @@ process.source = cms.Source("PoolSource", - debugFlag = cms.untracked.bool(True), - debugVebosity = cms.untracked.uint32(1), - - fileNames = cms.untracked.vstring( - '/store/relval/CMSSW_2_1_9/RelValSingleMuPt10/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v2/0000/2A00EECC-A185-DD11-93A9-000423D9517C.root' +# '/store/relval/CMSSW_2_1_9/RelValSingleMuPt10/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v2/0000/2A00EECC-A185-DD11-93A9-000423D9517C.root' + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/rechits/rechits1.root' + ) ) diff --git a/RecoLocalTracker/SiPixelRecHits/test/readRecHits_cfg.py b/RecoLocalTracker/SiPixelRecHits/test/readRecHits_cfg.py index 6de508131e53e..0c39be8b076a1 100644 --- a/RecoLocalTracker/SiPixelRecHits/test/readRecHits_cfg.py +++ b/RecoLocalTracker/SiPixelRecHits/test/readRecHits_cfg.py @@ -24,7 +24,9 @@ process.source = cms.Source("PoolSource", # fileNames = cms.untracked.vstring('file:/scratch/dkotlins/digis.root') - fileNames = cms.untracked.vstring('file:/scratch/dkotlins/promptrecoCosmics_1.root') + fileNames = cms.untracked.vstring( + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/rechits/rechits1.root' + ) ) @@ -33,34 +35,33 @@ fileName = cms.string('histo.root') ) - process.load("Configuration.StandardSequences.Geometry_cff") process.load("Configuration.StandardSequences.MagneticField_38T_cff") - -# what is this? # process.load("Configuration.StandardSequences.Services_cff") -# what is this? -#process.load("SimTracker.Configuration.SimTracker_cff") - -# needed for global transformation -process.load("Configuration.StandardSequences.FakeConditions_cff") - -# Initialize magnetic field -# include "MagneticField/Engine/data/volumeBasedMagneticField.cfi" -# Tracker SimGeometryXML -# include "Geometry/TrackerSimData/data/trackerSimGeometryXML.cfi" -# Tracker Geometry Builder -# include "Geometry/TrackerGeometryBuilder/data/trackerGeometry.cfi" -# Tracker Numbering Builder -# include "Geometry/TrackerNumberingBuilder/data/trackerNumberingGeometry.cfi" +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: +# 2012 +#process.GlobalTag.globaltag = 'GR_P_V40::All' +# MC 2913 +process.GlobalTag.globaltag = 'MC_70_V1::All' +# read rechits process.analysis = cms.EDAnalyzer("ReadPixelRecHit", Verbosity = cms.untracked.bool(True), src = cms.InputTag("siPixelRecHits"), ) +# test the DB object, works +process.load("RecoLocalTracker.SiPixelRecHits.PixelCPEESProducers_cff") +#process.load("RecoLocalTracker.SiPixelRecHits.SiPixelRecHits_cfi") +#process.load("CalibTracker.SiPixelESProducers.SiPixelFakeTemplateDBObjectESSource_cfi" +#process.load("CalibTracker.SiPixelESProducers.SiPixelFakeCPEGenericErrorParmESSource_cfi" +process.test = cms.EDAnalyzer("CPEAccessTester", +# PixelCPE = cms.string('PixelCPEGeneric'), + PixelCPE = cms.string('PixelCPETemplateReco'), +) process.p = cms.Path(process.analysis) +#process.p = cms.Path(process.test) diff --git a/SimTracker/SiPixelDigitizer/test/Digitze_Pixels_And_Strips.py b/SimTracker/SiPixelDigitizer/test/Digitze_Pixels_And_Strips.py new file mode 100644 index 0000000000000..ca6ac5fc28fb5 --- /dev/null +++ b/SimTracker/SiPixelDigitizer/test/Digitze_Pixels_And_Strips.py @@ -0,0 +1,173 @@ +############################################################################## + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("Test") + +process.load("FWCore.MessageLogger.MessageLogger_cfi") +# process.load("Configuration.StandardSequences.Geometry_cff") +process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") + +# process.load("SimTracker.Configuration.SimTracker_cff") +process.load("SimG4Core.Configuration.SimG4Core_cff") + +# for strips +process.load("CalibTracker.SiStripESProducers.SiStripGainSimESProducer_cfi") + +# process.load("SimGeneral.MixingModule.mixNoPU_cfi") + +from SimGeneral.MixingModule.aliases_cfi import * +from SimGeneral.MixingModule.mixObjects_cfi import * +# from SimGeneral.MixingModule.digitizers_cfi import * +from SimGeneral.MixingModule.pixelDigitizer_cfi import * +from SimGeneral.MixingModule.stripDigitizer_cfi import * +from SimGeneral.MixingModule.trackingTruthProducer_cfi import * + + +process.mix = cms.EDProducer("MixingModule", +# digitizers = cms.PSet(theDigitizers), +# digitizers = cms.PSet( +# mergedtruth = cms.PSet( +# trackingParticles +# ) +# ), + + digitizers = cms.PSet( + pixel = cms.PSet( + pixelDigitizer + ), + strip = cms.PSet( + stripDigitizer + ), + ), + +#theDigitizersValid = cms.PSet( +# pixel = cms.PSet( +# pixelDigitizer +# ), +# strip = cms.PSet( +# stripDigitizer +# ), +# ecal = cms.PSet( +# ecalDigitizer +# ), +# hcal = cms.PSet( +# hcalDigitizer +# ), +# castor = cms.PSet( +# castorDigitizer +# ), +# mergedtruth = cms.PSet( +# trackingParticles +# ) +#), + + + LabelPlayback = cms.string(' '), + maxBunch = cms.int32(3), + minBunch = cms.int32(-5), ## in terms of 25 ns + + bunchspace = cms.int32(25), + mixProdStep1 = cms.bool(False), + mixProdStep2 = cms.bool(False), + + playback = cms.untracked.bool(False), + useCurrentProcessOnly = cms.bool(False), + mixObjects = cms.PSet( + mixTracks = cms.PSet( + mixSimTracks + ), + mixVertices = cms.PSet( + mixSimVertices + ), + mixSH = cms.PSet( +# mixPixSimHits +# mixPixSimHits = cms.PSet( + input = cms.VInputTag(cms.InputTag("g4SimHits","TrackerHitsPixelBarrelHighTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelBarrelLowTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapHighTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapLowTof"), + cms.InputTag("g4SimHits","TrackerHitsTECHighTof"), + cms.InputTag("g4SimHits","TrackerHitsTECLowTof"), + cms.InputTag("g4SimHits","TrackerHitsTIBHighTof"), + cms.InputTag("g4SimHits","TrackerHitsTIBLowTof"), + cms.InputTag("g4SimHits","TrackerHitsTIDHighTof"), + cms.InputTag("g4SimHits","TrackerHitsTIDLowTof"), + cms.InputTag("g4SimHits","TrackerHitsTOBHighTof"), + cms.InputTag("g4SimHits","TrackerHitsTOBLowTof") + ), + type = cms.string('PSimHit'), + subdets = cms.vstring( + 'TrackerHitsPixelBarrelHighTof', + 'TrackerHitsPixelBarrelLowTof', + 'TrackerHitsPixelEndcapHighTof', + 'TrackerHitsPixelEndcapLowTof', + 'TrackerHitsTECHighTof', + 'TrackerHitsTECLowTof', + 'TrackerHitsTIBHighTof', + 'TrackerHitsTIBLowTof', + 'TrackerHitsTIDHighTof', + 'TrackerHitsTIDLowTof', + 'TrackerHitsTOBHighTof', + 'TrackerHitsTOBLowTof' + ), + crossingFrames = cms.untracked.vstring(), +# 'MuonCSCHits', +# 'MuonDTHits', +# 'MuonRPCHits'), +#) + ), + mixHepMC = cms.PSet( + mixHepMCProducts + ) + ) +) + +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", +# simMuonCSCDigis = cms.PSet( +# initialSeed = cms.untracked.uint32(1234567), +# engineName = cms.untracked.string('TRandom3') +# ), + mix = cms.PSet( + initialSeed = cms.untracked.uint32(1234567), + engineName = cms.untracked.string('TRandom3') + ) +) + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(-1) +) + +process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( + 'file:/afs/cern.ch/work/d/dkotlins/public//MC/mu/pt100/simhits/simHits1.root' + ) +) + +# Choose the global tag here: +# for v7.0 +#process.GlobalTag.globaltag = 'MC_70_V1::All' +process.GlobalTag.globaltag = 'START70_V1::All' +#process.GlobalTag.globaltag = 'DESIGN70_V1::All' + +process.o1 = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *','keep *_*_*_Test'), + fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/digis_trk/digis1.root') +# fileName = cms.untracked.string('file:dummy.root') +) + +process.g4SimHits.Generator.HepMCProductLabel = 'source' + +# modify digitizer parameters +#process.mix.digitizers.pixel.ThresholdInElectrons_BPix = 3500.0 + +#This process is to run the digitizer: +process.p1 = cms.Path(process.mix) + +process.outpath = cms.EndPath(process.o1) + + diff --git a/SimTracker/SiPixelDigitizer/test/PixelSimHitsTest.cc b/SimTracker/SiPixelDigitizer/test/PixelSimHitsTest.cc index 6c509e3bd118d..3b12af03e117a 100644 --- a/SimTracker/SiPixelDigitizer/test/PixelSimHitsTest.cc +++ b/SimTracker/SiPixelDigitizer/test/PixelSimHitsTest.cc @@ -108,11 +108,9 @@ class PixelSimHitsTest : public edm::EDAnalyzer { TH1F *hdetr, *hdetz, *hdetphi1, *hdetphi2, *hdetphi3; TH1F *hglobr1,*hglobr2,*hglobr3,*hglobz1, *hglobz2, *hglobz3; TH1F *hcolsB, *hrowsB, *hcolsF, *hrowsF; - TH1F *hglox1; TH2F *htest, *htest2, *htest3, *htest4, *htest5; - - TProfile *hp1, *hp2, *hp3, *hp4, *hp5; + //TProfile *hp1, *hp2, *hp3, *hp4, *hp5; #ifdef CHECK_GEOM float modulePositionZ[3][44][8]; @@ -172,12 +170,12 @@ void PixelSimHitsTest::beginJob() { hpixid = fs->make( "hpixid", "Pix det id", 10, 0., 10.); hpixsubid = fs->make( "hpixsubid", "Pix Barrel id", 10, 0., 10.); hlayerid = fs->make( "hlayerid", "Pix layer id", 10, 0., 10.); - hladder1id = fs->make( "hladder1id", "Ladder L1 id", 100, -0.5, 49.5); - hladder2id = fs->make( "hladder2id", "Ladder L2 id", 100, -0.5, 49.5); - hladder3id = fs->make( "hladder3id", "Ladder L3 id", 100, -0.5, 49.5); - hz1id = fs->make( "hz1id", "Z-index id L1", 10, 0., 10.); - hz2id = fs->make( "hz2id", "Z-index id L2", 10, 0., 10.); - hz3id = fs->make( "hz3id", "Z-index id L3", 10, 0., 10.); + hladder1id = fs->make( "hladder1id", "Ladder L1 id", 102, -25.5, 25.5); + hladder2id = fs->make( "hladder2id", "Ladder L2 id", 102, -25.5, 25.5); + hladder3id = fs->make( "hladder3id", "Ladder L3 id", 102, -25.5, 25.5); + hz1id = fs->make( "hz1id", "Z-index id L1", 10, -5., 5.); + hz2id = fs->make( "hz2id", "Z-index id L2", 10, -5., 5.); + hz3id = fs->make( "hz3id", "Z-index id L3", 10, -5., 5.); hthick1 = fs->make( "hthick1", "Det 1 Thinckess", 400, 0.,0.04); hthick2 = fs->make( "hthick2", "Det 2 Thinckess", 400, 0.,0.04); @@ -221,9 +219,9 @@ void PixelSimHitsTest::beginJob() { heloss2mu = fs->make( "heloss2mu", "Eloss mu l2", 100, 0., max_charge); heloss3mu = fs->make( "heloss3mu", "Eloss mu l3", 100, 0., max_charge); - htheta1 = fs->make( "htheta1", "Theta det1",350,0.0,3.5); - htheta2 = fs->make( "htheta2", "Theta det2",350,0.0,3.5); - htheta3 = fs->make( "htheta3", "Theta det3",350,0.0,3.5); + htheta1 = fs->make( "htheta1", "Theta l1",350,0.0,3.5); + htheta2 = fs->make( "htheta2", "Theta l2",350,0.0,3.5); + htheta3 = fs->make( "htheta3", "Theta l3",350,0.0,3.5); hphi1 = fs->make("hphi1","phi l1",1400,-3.5,3.5); hphi2 = fs->make("hphi2","phi l2",1400,-3.5,3.5); hphi3 = fs->make("hphi3","phi l3",1400,-3.5,3.5); @@ -251,18 +249,15 @@ void PixelSimHitsTest::beginJob() { hglobr3 = fs->make("hglobr3","global r3",150,0.,15.); hglobz3 = fs->make("hglobz3","global z3",540,-27.,27.); - hglox1 = fs->make("hglox1","global x in l1",200,-1.,1.); - - htest = fs->make("htest"," ",108,-27.,27.,35,-3.5,3.5); - htest2 = fs->make("htest2"," ",108,-27.,27.,60,0.,600.); + // layer 1 only + htest = fs->make("htest"," ",108,-27.,27.,35,-3.5,3.5); // global z versus local y + htest2 = fs->make("htest2"," ",108,-27.,27.,60,0.,600.); // global z versus eloss htest3 = fs->make("htest3"," ",240,-12.,12.,240,-12.,12.); // x-y plane //htest4 = fs->make("htest4"," ",80,-4.,4.,100,-5.,5.); - hp1 = fs->make("hp1"," ",50,0.,50.); // default option - hp2 = fs->make("hp2"," ",50,0.,50.," "); // option set to " " + //hp1 = fs->make("hp1"," ",50,0.,50.); // default option + //hp2 = fs->make("hp2"," ",50,0.,50.," "); // option set to " " //hp3 = fs->make("hp3"," ",50,0.,50.,-100.,100.); // - //hp4 = fs->make("hp4"," ",50,0.,50.); - // hp5 = fs->make("hp5"," ",50,0.,50.); #ifdef CHECK_GEOM // To get the module position @@ -360,8 +355,6 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, int cols = theGeomDet->specificTopology().ncolumns(); int rows = theGeomDet->specificTopology().nrows(); - hcolsB->Fill(float(cols)); - hrowsB->Fill(float(rows)); if(DEBUG) cout<<"det z/r "<Fill(float(cols)); + hrowsF->Fill(float(rows)); } else if(mode_ == "bpix") { // Barrel @@ -433,14 +429,17 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, // <Fill(float(cols)); + hrowsB->Fill(float(rows)); + } // end fb-bar #ifdef CHECK_GEOM // To get the module position - modulePositionR[layer-1][ladder-1][zindex-1] = detR; - modulePositionZ[layer-1][ladder-1][zindex-1] = detZ; - modulePositionPhi[layer-1][ladder-1][zindex-1] = detPhi; + modulePositionR[layer-1][ladderC-1][zindex-1] = detR; + modulePositionZ[layer-1][ladderC-1][zindex-1] = detZ; + modulePositionPhi[layer-1][ladderC-1][zindex-1] = detPhi; #endif // SimHit information @@ -469,14 +468,14 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, float ypos = (y+y2)/2.; float zpos = (z+z2)/2.; + if(pt>0.1) {goodHits++;} if(PRINT) { - if(pt>0.1) {cout<<" simhit: id "<0.1) cout<<" simhit: "; + else if(pid==11) cout<<" delta: "; + else cout<<" low pt (sec?): "; + cout<<" id "<Fill(eloss); else heloss1mu->Fill(eloss); hladder1id->Fill(float(ladder)); - hz1id->Fill(float(zindex)); + hz1id->Fill(float(module)); hthick1->Fill(dz); hlength1->Fill(y); hwidth1->Fill(x); @@ -559,7 +558,6 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, // hwidth1h->Fill(x); // if(pid==13 && p>1.) { // select primary muons with mom above 1. // hphi1h->Fill(phi); - // hglox1->Fill(gloX); // hglobr1h->Fill(gloR); // } // } else { @@ -576,10 +574,8 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, //if(pid!=11 && moduleDirectionUp) hladder1idUp->Fill(float(ladder)); - //if(ladder==6) htest4->Fill(xpos,gloX); - - hp1->Fill(float(ladder),detR,1); - hp2->Fill(float(ladder),detPhi); + //hp1->Fill(float(ladder),detR,1); + //hp2->Fill(float(ladder),detPhi); hdetphi1->Fill(detPhi); } else if(layer==2) { @@ -590,7 +586,7 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, if(abs(pid)==11) heloss2e->Fill(eloss); else heloss2mu->Fill(eloss); hladder2id->Fill(float(ladder)); - hz2id->Fill(float(zindex)); + hz2id->Fill(float(module)); hthick2->Fill(dz); hlength2->Fill(y); hwidth2->Fill(x); @@ -619,7 +615,7 @@ void PixelSimHitsTest::analyze(const edm::Event& iEvent, else heloss3mu->Fill(eloss); hladder3id->Fill(float(ladder)); - hz3id->Fill(float(zindex)); + hz3id->Fill(float(module)); hthick3->Fill(dz); hlength3->Fill(y); hwidth3->Fill(x); diff --git a/SimTracker/SiPixelDigitizer/test/gen.py b/SimTracker/SiPixelDigitizer/test/gen.py new file mode 100644 index 0000000000000..1d385446733fa --- /dev/null +++ b/SimTracker/SiPixelDigitizer/test/gen.py @@ -0,0 +1,201 @@ +# +import FWCore.ParameterSet.Config as cms + +process = cms.Process("GenTest") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(20000) +) + +# Input source +process.source = cms.Source("EmptySource") + +process.options = cms.untracked.PSet( + +) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') + +process.load('SimGeneral.MixingModule.mixNoPU_cfi') + +process.load('Configuration.StandardSequences.Generator_cff') + +process.load('GeneratorInterface.Core.genFilterSummary_cff') + +process.load('Configuration.StandardSequences.SimIdeal_cff') +# process.load('Configuration/StandardSequences/Sim_cff') + +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", + generator = cms.PSet( +# 1st file +# initialSeed = cms.untracked.uint32(100000000), +# 2nd file +# initialSeed = cms.untracked.uint32(123456789), +# 3rd file +# initialSeed = cms.untracked.uint32(200000000), +# 4th file +# initialSeed = cms.untracked.uint32(300000000), +# 5th file +# initialSeed = cms.untracked.uint32(400000000), +# 6th file + initialSeed = cms.untracked.uint32(500000000), + + engineName = cms.untracked.string('HepJamesRandom') + ) +) + + +########################################### +# list of possible events to be generated # +########################################### +# +# single muons +process.generator = cms.EDProducer("FlatRandomPtGunProducer", + PGunParameters = cms.PSet( + MaxPt = cms.double(100.1), + MinPt = cms.double(99.9), + PartID = cms.vint32(13), + MaxEta = cms.double(2.5), + MaxPhi = cms.double(3.14159265359), + MinEta = cms.double(-2.5), + MinPhi = cms.double(-3.14159265359) ## in radians + + ), + Verbosity = cms.untracked.int32(0), ## set to 1 (or greater) for printouts + + psethack = cms.string('single muon pt 100'), + AddAntiParticle = cms.bool(False), # True makes to muons per event +# firstRun = cms.untracked.uint32(1) + firstRun = cms.untracked.uint32(2) +# firstRun = cms.untracked.uint32(3) +# firstRun = cms.untracked.uint32(4) +) + +## # choose particle +#process.generator.PGunParameters.PartID[0] = 13 +## # example: for 4 muons to test with vertex +## #process.generator.PGunParameters.PartID = cms.untracked.vint32(13,-13,13,-13) +## # example: for opposite sign back-to-back dimuon pairs set to True +## # define limits for Pt +#process.generator.PGunParameters.MinPt = 40.0 +#process.generator.PGunParameters.MaxPt = 50.0 +## # define limits for Pseudorapidity +#process.generator.PGunParameters.MinEta = -3 +#process.generator.PGunParameters.MaxEta = 3 +#process.source.firstRun = cms.untracked.uint32(10) +#process.source.firstEvent = cms.untracked.uint32(9001) +#process.source.firstLuminosityBlock = cms.untracked.uint32(10) + + +#process.MessageLogger = cms.Service("MessageLogger", +# debugModules = cms.untracked.vstring('PixelDigisTest'), +# destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), +# cout = cms.untracked.PSet( +# threshold = cms.untracked.string('ERROR') +# ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +#) + +#process.load("Geometry.TrackerSimData.trackerSimGeometryXML_cfi") +#process.load("Geometry.CMSCommonData.cmsSimIdealGeometryXML_cfi") + +#process.load("Configuration.StandardSequences.Geometry_cff") +#process.load("Configuration.StandardSequences.GeometryIdeal_cff") +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.GeometrySimDB_cff') + + +# does using an empty PixelSkimmedGeometry.txt file speeds up job with lots more channels? +#process.siPixelFakeGainOfflineESSource = cms.ESSource("SiPixelFakeGainOfflineESSource", +# file = cms.FileInPath('SLHCUpgradeSimulations/Geometry/data/longbarrel/PixelSkimmedGeometry_empty.txt') +#) +#process.es_prefer_fake_gain = cms.ESPrefer("SiPixelFakeGainOfflineESSource","siPixelFakeGainOfflineESSource") + +#process.siPixelFakeLorentzAngleESSource = cms.ESSource("SiPixelFakeLorentzAngleESSource", +# file = cms.FileInPath('SLHCUpgradeSimulations/Geometry/data/longbarrel/PixelSkimmedGeometry.txt') +#) +#process.es_prefer_fake_lorentz = cms.ESPrefer("SiPixelFakeLorentzAngleESSource","siPixelFakeLorentzAngleESSource") + + +############################## +# magnetic field in solenoid # +############################## +# +process.load('Configuration.StandardSequences.MagneticField_38T_cff') + +# Parametrized magnetic field (new mapping, 4.0 and 3.8T) +process.VolumeBasedMagneticFieldESProducer.useParametrizedTrackerField = True + +######################### +# event vertex smearing # +######################### +# +#process.load("Configuration.StandardSequences.VtxSmearedGauss_cff") +#process.load("Configuration.StandardSequences.VtxSmearedBetafuncEarlyCollision_cff") +process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic8TeVCollision_cfi') + + +########### +# what is this? +# process.load("Configuration.StandardSequences.Services_cff") + +# what is this? +#process.load("SimTracker.Configuration.SimTracker_cff") + +# needed for global transformation +# this crashes +# process.load("Configuration.StandardSequences.FakeConditions_cff") + +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")# Choose the global tag here: +#process.GlobalTag.globaltag = 'MC_53_V15::All' +#process.GlobalTag.globaltag = 'DESIGN53_V15::All' +#process.GlobalTag.globaltag = 'START53_V15::All' +# ideal +process.GlobalTag.globaltag = 'MC_70_V1::All' +# realistiv alignment and calibrations +#process.GlobalTag.globaltag = 'START70_V1::All' + +#process.load("Geometry.TrackerGeometryBuilder.trackerGeometry_cfi") +#process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi") +#process.load("Configuration.StandardSequences.MagneticField_cff") +# include "MagneticField/Engine/data/volumeBasedMagneticField.cfi" + + +############################### +# global output of simulation # +############################### +# +process.o1 = cms.OutputModule( + "PoolOutputModule", +# definition of branches to keep or drop + outputCommands = cms.untracked.vstring('keep *','drop PCaloHits_*_*_*','drop *_*_MuonCSCHits_*', + 'drop *_*_MuonDTHits_*','drop *_*_MuonRPCHits_*', + 'drop *_*_MuonPLTHits_*','drop *_*_TotemHitsRP_*', + 'drop *_*_TotemHitsT1_*','drop *_*_TotemHitsT2Gem_*', + ), + +# definition of output file (full path) + fileName = cms.untracked.string('/afs/cern.ch/user/d/dkotlins/work/MC/mu/pt100/simhits/simHits6.root') +) + +# +process.outpath = cms.EndPath(process.o1) + +#process.simulation_step = cms.Path(process.psim) +# process.digitisation_step = cms.Path(process.pdigi) + +#process.p = cms.Path(process.generator*process.genParticles) +process.p = cms.Path(process.generator*process.genParticles*process.psim) + +process.schedule = cms.Schedule(process.p,process.outpath) diff --git a/SimTracker/SiPixelDigitizer/test/runDigisToTracks.py b/SimTracker/SiPixelDigitizer/test/runDigisToTracks.py new file mode 100644 index 0000000000000..b6e5caff3a651 --- /dev/null +++ b/SimTracker/SiPixelDigitizer/test/runDigisToTracks.py @@ -0,0 +1,177 @@ +# produce pixel cluster & rechits from digia +# works directly or through raw +# +# +############################################################################## + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("TrackTest") + +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +# clusterizer +process.load("RecoLocalTracker.Configuration.RecoLocalTracker_cff") + +# for raw +#process.load("EventFilter.SiPixelRawToDigi.SiPixelDigiToRaw_cfi") +#process.load("EventFilter.SiPixelRawToDigi.SiPixelRawToDigi_cfi") +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') + +# needed for pixel RecHits (templates?) +process.load("Configuration.StandardSequences.Reconstruction_cff") +#process.load("RecoTracker.Configuration.RecoTracker_cff") + +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1) +) + +process.MessageLogger = cms.Service("MessageLogger", +# debugModules = cms.untracked.vstring('SiPixelClusterizer'), + debugModules = cms.untracked.vstring('SiPixelRecHits'), + destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), + cout = cms.untracked.PSet( +# threshold = cms.untracked.string('INFO') + threshold = cms.untracked.string('ERROR') +# threshold = cms.untracked.string('WARNING') +# threshold = cms.untracked.string('DEBUG') + ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +) +# get the files from DBS: +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/digis_trk/digis1.root' + ) +) + +# Choose the global tag here: +process.GlobalTag.globaltag = 'MC_70_V1::All' + +#process.PoolDBESSource = cms.ESSource("PoolDBESSource", +# BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), +# DBParameters = cms.PSet( +# messageLevel = cms.untracked.int32(0), +# authenticationPath = cms.untracked.string('') +# ), +# timetype = cms.string('runnumber'), +# toGet = cms.VPSet(cms.PSet( +# record = cms.string('SiPixelQualityRcd'), +# tag = cms.string('SiPixelBadModule_test') +# )), +# connect = cms.string('sqlite_file:test.db') +#) +# +# To use a test DB instead of the official pixel object DB tag: +#process.customDead = cms.ESSource("PoolDBESSource", process.CondDBSetup, connect = cms.string('sqlite_file:/afs/cern.ch/user/v/vesna/Digitizer/dead_20100901.db'), toGet = cms.VPSet(cms.PSet(record = cms.string('SiPixelQualityRcd'), tag = cms.string('dead_20100901')))) +#process.es_prefer_customDead = cms.ESPrefer("PoolDBESSource","customDead") + + +process.o1 = cms.OutputModule("PoolOutputModule", +# outputCommands = cms.untracked.vstring('drop *','keep *_*_*_TrackTest'), + fileName = cms.untracked.string('file:tracks.root'), +# fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/tracks/tracks1.root'), +# splitLevel = cms.untracked.int32(0), +# eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + outputCommands = process.RECOSIMEventContent.outputCommands, + dataset = cms.untracked.PSet( + filterName = cms.untracked.string(''), + dataTier = cms.untracked.string('RECO') + ) + +) + +#process.Timing = cms.Service("Timing") +#process.SimpleMemoryCheck = cms.Service("SimpleMemoryCheck") + +# My +# modify clusterie parameters +#process.siPixelClusters.ClusterThreshold = 4000.0 + +# direct clusterization (no raw step) +# label of digis +#process.siPixelClusters.src = 'mix' + +# plus pixel clusters (OK) +#process.p1 = cms.Path(process.siPixelClusters) + +# plus pixel rechits (OK) +#process.p1 = cms.Path(process.pixeltrackerlocalreco) + +# clusterize through raw (OK) +# for digi to raw +process.siPixelRawData.InputLabel = 'mix' +process.SiStripDigiToRaw.InputModuleLabel = 'mix' +# for Raw2digi for simulations +process.siPixelDigis.InputLabel = 'siPixelRawData' +process.siStripDigis.ProductLabel = 'SiStripDigiToRaw' +# for digi to clu +process.siPixelClusters.src = 'siPixelDigis' + +# pixel only +#process.p1 = cms.Path(process.siPixelRawData) +#process.p1 = cms.Path(process.siPixelRawData*process.siPixelDigis) +#process.p1 = cms.Path(process.siPixelRawData*process.siPixelDigis*process.pixeltrackerlocalreco) + +# with strips ok +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw) +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis) + +# runs ok +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis*process.trackerlocalreco) + +# runs ok +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot) + +# runs ok +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.MeasurementTrackerEvent) + +# runs ok +#process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.MeasurementTrackerEvent*process.recopixelvertexing) + +# copy the sequence below from +# RecoTracker/IterativeTracking/python/iterativeTk_cff.py +process.myTracking = cms.Sequence(process.InitialStep* + process.LowPtTripletStep* + process.PixelPairStep* + process.DetachedTripletStep* + process.MixedTripletStep* + process.PixelLessStep* + process.TobTecStep* + process.earlyGeneralTracks* + #process.muonSeededStep* + process.preDuplicateMergingGeneralTracks* + process.generalTracksSequence* + process.ConvStep* + process.conversionStepTracks + ) + +# run full tracking +# trackingGlobalReco does not work, needs EarlyMuons for muon seeding. +# ckftracks & iterTracking does not work as well (same problem). +process.p1 = cms.Path(process.siPixelRawData*process.SiStripDigiToRaw*process.siPixelDigis*process.siStripDigis*process.trackerlocalreco*process.offlineBeamSpot*process.recopixelvertexing*process.MeasurementTrackerEvent*process.myTracking*process.vertexreco) + + +# unused +# Path and EndPath definitions +#process.raw2digi_step = cms.Path(process.RawToDigi) +#process.reconstruction_step = cms.Path(process.reconstruction) +#process.endjob_step = cms.EndPath(process.endOfProcess) +#process.RECOSIMoutput_step = cms.EndPath(process.RECOSIMoutput) +# Schedule definition +#process.schedule = cms.Schedule(process.raw2digi_step,process.reconstruction_step,process.endjob_step,process.RECOSIMoutput_step) + +process.outpath = cms.EndPath(process.o1) diff --git a/SimTracker/SiPixelDigitizer/test/runPixelRH.py b/SimTracker/SiPixelDigitizer/test/runPixelRH.py new file mode 100644 index 0000000000000..e3cc27c551cad --- /dev/null +++ b/SimTracker/SiPixelDigitizer/test/runPixelRH.py @@ -0,0 +1,116 @@ +# produce pixel cluster & rechits from digia +# works directly or through raw +# +# +############################################################################## + +import FWCore.ParameterSet.Config as cms + +process = cms.Process("ClusTest") + +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.load("Configuration.Geometry.GeometryIdeal_cff") +process.load("Configuration.StandardSequences.MagneticField_38T_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.load("Configuration.StandardSequences.Services_cff") + +# clusterizer +process.load("RecoLocalTracker.Configuration.RecoLocalTracker_cff") + +# for raw +#process.load("EventFilter.SiPixelRawToDigi.SiPixelDigiToRaw_cfi") +#process.load("EventFilter.SiPixelRawToDigi.SiPixelRawToDigi_cfi") +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') + + +# needed for pixel RecHits (templates?) +process.load("Configuration.StandardSequences.Reconstruction_cff") + +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(-1) +) + +process.MessageLogger = cms.Service("MessageLogger", + debugModules = cms.untracked.vstring('SiPixelClusterizer'), + destinations = cms.untracked.vstring('cout'), +# destinations = cms.untracked.vstring("log","cout"), + cout = cms.untracked.PSet( +# threshold = cms.untracked.string('INFO') +# threshold = cms.untracked.string('ERROR') + threshold = cms.untracked.string('WARNING') + ) +# log = cms.untracked.PSet( +# threshold = cms.untracked.string('DEBUG') +# ) +) +# get the files from DBS: +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/digis/digis1.root' + ) +) + +# Choose the global tag here: +process.GlobalTag.globaltag = 'MC_70_V1::All' + +#process.PoolDBESSource = cms.ESSource("PoolDBESSource", +# BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), +# DBParameters = cms.PSet( +# messageLevel = cms.untracked.int32(0), +# authenticationPath = cms.untracked.string('') +# ), +# timetype = cms.string('runnumber'), +# toGet = cms.VPSet(cms.PSet( +# record = cms.string('SiPixelQualityRcd'), +# tag = cms.string('SiPixelBadModule_test') +# )), +# connect = cms.string('sqlite_file:test.db') +#) +# +# To use a test DB instead of the official pixel object DB tag: +#process.customDead = cms.ESSource("PoolDBESSource", process.CondDBSetup, connect = cms.string('sqlite_file:/afs/cern.ch/user/v/vesna/Digitizer/dead_20100901.db'), toGet = cms.VPSet(cms.PSet(record = cms.string('SiPixelQualityRcd'), tag = cms.string('dead_20100901')))) +#process.es_prefer_customDead = cms.ESPrefer("PoolDBESSource","customDead") + + +process.o1 = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *','keep *_*_*_ClusTest'), +# fileName = cms.untracked.string('file:clus.root') + fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100_pre10/clus/clus1.root') +# fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/clus/clus1.root') +# fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/rechits/rechits1.root') +) + +#process.Timing = cms.Service("Timing") +#process.SimpleMemoryCheck = cms.Service("SimpleMemoryCheck") + +# My +# modify clusterie parameters +#process.siPixelClusters.ClusterThreshold = 4000.0 + +# DIRECT +# direct clusterization (no raw step) +# label of digis +process.siPixelClusters.src = 'mix' + +# plus pixel clusters (OK) +#process.p1 = cms.Path(process.siPixelClusters) +# plus pixel rechits (OK) +process.p1 = cms.Path(process.pixeltrackerlocalreco) + +# RAW +# clusterize through raw (OK) +# for Raw2digi for simulations +#process.siPixelRawData.InputLabel = 'mix' +#process.siPixelDigis.InputLabel = 'siPixelRawData' +# process.siStripDigis.ProductLabel = 'SiStripDigiToRaw' +#process.siPixelClusters.src = 'siPixelDigis' + +#process.p1 = cms.Path(process.siPixelRawData) +#process.p1 = cms.Path(process.siPixelRawData*process.siPixelDigis) +#process.p1 = cms.Path(process.siPixelRawData*process.siPixelDigis*process.pixeltrackerlocalreco) + +process.outpath = cms.EndPath(process.o1) diff --git a/SimTracker/SiPixelDigitizer/test/testPixelDigitizer.py b/SimTracker/SiPixelDigitizer/test/testPixelDigitizer.py index c86161e48a872..8713e91fd7311 100644 --- a/SimTracker/SiPixelDigitizer/test/testPixelDigitizer.py +++ b/SimTracker/SiPixelDigitizer/test/testPixelDigitizer.py @@ -1,202 +1,168 @@ -# python script Pixel Digitizer testing: -# author: V. Cuplov -# -# Process available: -# 1) simulation of hits in the Tracker -# 2) digitization + validation -# 3) reconstruction or clusterization + validation -# 4) tracks + validation -# ############################################################################## import FWCore.ParameterSet.Config as cms -process = cms.Process("DigiTest") +process = cms.Process("Test") process.load("FWCore.MessageLogger.MessageLogger_cfi") -process.load("Configuration.StandardSequences.Geometry_cff") +# process.load("Configuration.StandardSequences.Geometry_cff") +process.load("Configuration.Geometry.GeometryIdeal_cff") process.load("Configuration.StandardSequences.MagneticField_38T_cff") process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") process.load("Configuration.StandardSequences.Services_cff") process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") -process.load("Validation.TrackerHits.trackerHitsValidation_cff") -process.load("SimGeneral.MixingModule.mixNoPU_cfi") -process.load("SimTracker.Configuration.SimTracker_cff") -process.load("Validation.TrackerDigis.trackerDigisValidation_cff") + +# from v7 +#process.load("SimGeneral.MixingModule.pixelDigitizer_cfi") +# process.load("SimTracker.Configuration.SimTracker_cff") process.load("SimG4Core.Configuration.SimG4Core_cff") -process.load("Configuration.StandardSequences.Reconstruction_cff") -## For 31X use: -process.load("RecoTracker.TrackProducer.TrackRefitters_cff") -## for older releases -#process.load("RecoTracker.TrackProducer.RefitterWithMaterial_cff") +# process.load("SimGeneral.MixingModule.mixNoPU_cfi") +from SimGeneral.MixingModule.aliases_cfi import * +from SimGeneral.MixingModule.mixObjects_cfi import * +# from SimGeneral.MixingModule.digitizers_cfi import * +from SimGeneral.MixingModule.pixelDigitizer_cfi import * +from SimGeneral.MixingModule.stripDigitizer_cfi import * +from SimGeneral.MixingModule.trackingTruthProducer_cfi import * + +process.mix = cms.EDProducer("MixingModule", +# digitizers = cms.PSet(theDigitizers), +# digitizers = cms.PSet( +# mergedtruth = cms.PSet( +# trackingParticles +# ) +# ), + + digitizers = cms.PSet( + pixel = cms.PSet( + pixelDigitizer + ), +# strip = cms.PSet( +# stripDigitizer +# ), + ), + +#theDigitizersValid = cms.PSet( +# pixel = cms.PSet( +# pixelDigitizer +# ), +# strip = cms.PSet( +# stripDigitizer +# ), +# ecal = cms.PSet( +# ecalDigitizer +# ), +# hcal = cms.PSet( +# hcalDigitizer +# ), +# castor = cms.PSet( +# castorDigitizer +# ), +# mergedtruth = cms.PSet( +# trackingParticles +# ) +#), + + LabelPlayback = cms.string(' '), + maxBunch = cms.int32(3), + minBunch = cms.int32(-5), ## in terms of 25 ns + + bunchspace = cms.int32(25), + mixProdStep1 = cms.bool(False), + mixProdStep2 = cms.bool(False), + + playback = cms.untracked.bool(False), + useCurrentProcessOnly = cms.bool(False), + mixObjects = cms.PSet( + mixTracks = cms.PSet( + mixSimTracks + ), + mixVertices = cms.PSet( + mixSimVertices + ), + mixSH = cms.PSet( +# mixPixSimHits +# mixPixSimHits = cms.PSet( + input = cms.VInputTag(cms.InputTag("g4SimHits","TrackerHitsPixelBarrelHighTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelBarrelLowTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapHighTof"), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapLowTof"), +# cms.InputTag("g4SimHits","TrackerHitsTECHighTof"), +# cms.InputTag("g4SimHits","TrackerHitsTECLowTof"), +# cms.InputTag("g4SimHits","TrackerHitsTIBHighTof"), +# cms.InputTag("g4SimHits","TrackerHitsTIBLowTof"), +# cms.InputTag("g4SimHits","TrackerHitsTIDHighTof"), +# cms.InputTag("g4SimHits","TrackerHitsTIDLowTof"), +# cms.InputTag("g4SimHits","TrackerHitsTOBHighTof"), +# cms.InputTag("g4SimHits","TrackerHitsTOBLowTof") + ), + type = cms.string('PSimHit'), + subdets = cms.vstring( + 'TrackerHitsPixelBarrelHighTof', + 'TrackerHitsPixelBarrelLowTof', + 'TrackerHitsPixelEndcapHighTof', + 'TrackerHitsPixelEndcapLowTof', +# 'TrackerHitsTECHighTof', +# 'TrackerHitsTECLowTof', +# 'TrackerHitsTIBHighTof', +# 'TrackerHitsTIBLowTof', +# 'TrackerHitsTIDHighTof', +# 'TrackerHitsTIDLowTof', +# 'TrackerHitsTOBHighTof', +# 'TrackerHitsTOBLowTof' + ), + crossingFrames = cms.untracked.vstring(), +# 'MuonCSCHits', +# 'MuonDTHits', +# 'MuonRPCHits'), +#) + ), + mixHepMC = cms.PSet( + mixHepMCProducts + ) + ) +) -process.load("Validation.TrackerRecHits.trackerRecHitsValidation_cff") -process.load("SimGeneral.TrackingAnalysis.trackingParticles_cfi") -process.load("Validation.TrackingMCTruth.trackingTruthValidation_cfi") -process.load("Validation.RecoTrack.TrackValidation_cff") -process.load("Validation.RecoTrack.SiTrackingRecHitsValid_cff") -process.load("RecoLocalTracker.Configuration.RecoLocalTracker_cff") -process.load("EventFilter.SiPixelRawToDigi.SiPixelDigiToRaw_cfi") -process.load("EventFilter.SiPixelRawToDigi.SiPixelRawToDigi_cfi") +process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", +# simMuonCSCDigis = cms.PSet( +# initialSeed = cms.untracked.uint32(1234567), +# engineName = cms.untracked.string('TRandom3') +# ), + mix = cms.PSet( + initialSeed = cms.untracked.uint32(1234567), + engineName = cms.untracked.string('TRandom3') + ) +) process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) -##process.MessageLogger = cms.Service("MessageLogger", -## debugModules = cms.untracked.vstring('SiPixelDigitizer'), -## destinations = cms.untracked.vstring('cout'), -### destinations = cms.untracked.vstring("log","cout"), -## cout = cms.untracked.PSet( -### threshold = cms.untracked.string('INFO') -### threshold = cms.untracked.string('ERROR') -## threshold = cms.untracked.string('WARNING') -## ) -### log = cms.untracked.PSet( -### threshold = cms.untracked.string('DEBUG') -### ) -##) -### get the files from DBS: -###process.source = cms.Source("PoolSource", -### fileNames = cms.untracked.vstring('/store/relval/CMSSW_3_1_0_pre9/RelValSingleMuPt10/GEN-SIM-RECO/IDEAL_31X_v1/0007/B41EF45B-D14E-DE11-BC68-001D09F24600.root','/store/relval/CMSSW_3_1_0_pre9/RelValSingleMuPt10/GEN-SIM-RECO/IDEAL_31X_v1/0007/682FF80C-524F-DE11-9159-001D09F2543D.root','/store/relval/CMSSW_3_1_0_pre9/RelValSingleMuPt10/GEN-SIM-RECO/IDEAL_31X_v1/0007/0E9B84FD-D24E-DE11-94D9-001617C3B70E.root') -###) - -readFiles = cms.untracked.vstring() -secFiles = cms.untracked.vstring() -process.source = cms.Source("PoolSource", fileNames = readFiles, secondaryFileNames = secFiles) - - -readFiles.extend( [ -###Muons ls -ltr /scratch/jbutt/MuPt100/* | awk '{print ",\x27 file:" $9 "\x27 "}' -# 'file:/scratch/jbutt/MuPt100/3C1F82E7-6605-E111-B612-002354EF3BE2.root' -# ,'file:/scratch/jbutt/MuPt100/82C9F14B-B402-E111-87ED-002618943821.root' -###QCD ls -ltr /scratch/jbutt/QCD80120/* | awk '{print ",\x27 file:" $9 "\x27 "}' - ' file:/scratch/jbutt/QCD80120/16B15CEC-6605-E111-91B4-002618FDA204.root' - ,' file:/scratch/jbutt/QCD80120/423FA9C8-E606-E111-937B-002354EF3BE1.root' - ,' file:/scratch/jbutt/QCD80120/C2ACFE6C-2603-E111-B58F-001A9281173E.root' - ,' file:/scratch/jbutt/QCD80120/8C340B79-8002-E111-89EE-002354EF3BE0.root' - ,' file:/scratch/jbutt/QCD80120/7441D424-1702-E111-A600-0018F3D09696.root' - ,' file:/scratch/jbutt/QCD80120/B08F23A4-6901-E111-AA0C-002618943898.root' -####TTbar ls -ltr /scratch/jbutt/TTbar/* | awk '{print ",\x27 file:" $9 "\x27 "}' -# ' file:/scratch/jbutt/TTbar/1676ABF5-9004-E111-9B9B-002618FDA208.root' -#,' file:/scratch/jbutt/TTbar/5C73216C-2603-E111-AB90-002618FDA279.root' -#,' file:/scratch/jbutt/TTbar/F44EC1F2-6605-E111-85FF-003048678FDE.root' -#,' file:/scratch/jbutt/TTbar/C24DF25A-B302-E111-8AC4-001A92971BCA.root' -#,' file:/scratch/jbutt/TTbar/16DA979B-1A02-E111-9C54-0018F3D095F2.root' -#,' file:/scratch/jbutt/TTbar/34C0859D-1702-E111-835F-00304867BEDE.root' -#,' file:/scratch/jbutt/TTbar/2E90C1B8-6501-E111-8C4F-001A92971B08.root' -#,' file:/scratch/jbutt/TTbar/B4E41354-6601-E111-8D1B-00304867BED8.root' - ] ); - -secFiles.extend( [ -####Muons -# 'file:/scratch/jbutt/MuPt100/FE1E1645-B602-E111-9729-003048679168.root' -# ,'file:/scratch/jbutt/MuPt100/C6C414E9-6605-E111-854B-002618FDA21D.root' -# ,'file:/scratch/jbutt/MuPt100/B2649223-1702-E111-8399-001A9281174A.root' -# ,'file:/scratch/jbutt/MuPt100/3CFBFCFE-7E01-E111-B24B-0018F3D095FA.root' -###QCD - ' file:/scratch/jbutt/QCD80120/1228DEDC-2403-E111-914A-002618943896.root' - ,' file:/scratch/jbutt/QCD80120/182CF366-2603-E111-9F72-003048D42D92.root' - ,' file:/scratch/jbutt/QCD80120/72EE3DEF-6605-E111-8301-002618FDA248.root' - ,' file:/scratch/jbutt/QCD80120/DC1E3FC3-E606-E111-9C3E-00304867BEE4.root' - ,' file:/scratch/jbutt/QCD80120/EE35547C-8002-E111-A832-001A92810AAA.root' - ,' file:/scratch/jbutt/QCD80120/94B72EB2-1602-E111-8B8A-001A92810AE4.root' - ,' file:/scratch/jbutt/QCD80120/9C0615B6-1602-E111-B893-001A928116E0.root' - ,' file:/scratch/jbutt/QCD80120/E45BE72A-1902-E111-9BDE-00261894389E.root' - ,' file:/scratch/jbutt/QCD80120/4C14BA8B-6601-E111-8D01-001A92971BB8.root' - ,' file:/scratch/jbutt/QCD80120/625F1E38-6801-E111-A244-0030486792AC.root' - ,' file:/scratch/jbutt/QCD80120/849BD437-6801-E111-B8F1-002618943950.root' -###TTbar -# ' file:/scratch/jbutt/TTbar/4C4CEEEC-9004-E111-81CB-0026189438E4.root' -#,' file:/scratch/jbutt/TTbar/6851646C-2603-E111-A121-001BFCDBD11E.root' -#,' file:/scratch/jbutt/TTbar/8A6EC4EE-6605-E111-BAEB-0026189438A0.root' -#,' file:/scratch/jbutt/TTbar/86C08ADC-9E02-E111-BF3E-002618943916.root' -#,' file:/scratch/jbutt/TTbar/306D3C24-1802-E111-BB8C-002618943935.root' -#,' file:/scratch/jbutt/TTbar/628DBD21-1A02-E111-A5AE-0018F3D096F8.root' -#,' file:/scratch/jbutt/TTbar/D44E792A-1802-E111-BC00-001A92971B74.root' -#,' file:/scratch/jbutt/TTbar/E05D2324-1802-E111-8B0A-00248C0BE014.root' -#,' file:/scratch/jbutt/TTbar/12938537-6401-E111-A1BB-0018F3D096EE.root' -#,' file:/scratch/jbutt/TTbar/2A2223CF-6401-E111-8DAA-001A928116FA.root' -#,' file:/scratch/jbutt/TTbar/568D9338-6201-E111-9748-001A92810AEA.root' -#,' file:/scratch/jbutt/TTbar/6C1A3035-6701-E111-A24A-002618943807.root' -#,' file:/scratch/jbutt/TTbar/6CF830DB-6E01-E111-8E98-001A92810AD2.root' -#,' file:/scratch/jbutt/TTbar/B2972B10-6401-E111-AE9C-0018F3D09634.root' -#,' file:/scratch/jbutt/TTbar/C0A60C2F-6501-E111-932F-001A92810ADE.root' -#,' file:/scratch/jbutt/TTbar/CC654D60-6401-E111-8A27-00261894385A.root' -#,' file:/scratch/jbutt/TTbar/D0FCC289-6B01-E111-8AEB-001BFCDBD15E.root' - ] ) - - - -#process.PoolDBESSource = cms.ESSource("PoolDBESSource", -# BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), -# DBParameters = cms.PSet( -# messageLevel = cms.untracked.int32(0), -# authenticationPath = cms.untracked.string('') -# ), -# timetype = cms.string('runnumber'), -# toGet = cms.VPSet(cms.PSet( -# record = cms.string('SiPixelQualityRcd'), -# tag = cms.string('SiPixelBadModule_test') -# )), -# connect = cms.string('sqlite_file:test.db') -#) - +process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( + 'file:/afs/cern.ch/work/d/dkotlins/public//MC/mu/pt100/simhits/simHits1.root' + ) +) # Choose the global tag here: -#process.GlobalTag.globaltag = 'MC_38Y_V8::All' -process.GlobalTag.globaltag = 'MC_50_V0::All' - -process.load("CondCore.DBCommon.CondDBSetup_cfi") - -# To use a test DB instead of the official pixel object DB tag: -#process.customDead = cms.ESSource("PoolDBESSource", process.CondDBSetup, connect = cms.string('sqlite_file:/afs/cern.ch/user/v/vesna/Digitizer/dead_20100901.db'), toGet = cms.VPSet(cms.PSet(record = cms.string('SiPixelQualityRcd'), tag = cms.string('dead_20100901')))) -#process.es_prefer_customDead = cms.ESPrefer("PoolDBESSource","customDead") - +# for v7.0 +process.GlobalTag.globaltag = 'MC_70_V1::All' process.o1 = cms.OutputModule("PoolOutputModule", - outputCommands = cms.untracked.vstring('drop *','keep *_*_*_DigiTest'), -# fileName = cms.untracked.string('rfio:/castor/cern.ch/user/v/vesna/testDigis.root') - fileName = cms.untracked.string('file:dummy.root') + outputCommands = cms.untracked.vstring('drop *','keep *_*_*_Test'), + fileName = cms.untracked.string('file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/digis/digis1.root') +# fileName = cms.untracked.string('file:dummy.root') ) -#process.Timing = cms.Service("Timing") - -#process.SimpleMemoryCheck = cms.Service("SimpleMemoryCheck") - -# Simulation of hits in the Tracker: -#process.simhits = cms.Sequence(process.g4SimHits*process.trackerHitsValidation) - -# Digitization: -process.digis = cms.Sequence(process.simSiPixelDigis*process.pixelDigisValid) - -# Reconstruction or Clusterization: -#process.rechits = cms.Sequence(process.pixeltrackerlocalreco*process.pixRecHitsValid) -process.rechits = cms.Sequence(process.pixeltrackerlocalreco*process.siStripMatchedRecHits*process.pixRecHitsValid) - -#process.simSiPixelDigis.ChargeVCALSmearing=True -process.tracks = cms.Sequence(process.offlineBeamSpot*process.recopixelvertexing*process.trackingParticles*process.trackingTruthValid*process.ckftracks*process.trackerRecHitsValidation) -process.trackinghits = cms.Sequence(process.TrackRefitter*process.trackingRecHitsValid) - -# To use when you want to look at Residuals and Pulls: -#process.p1 = cms.Path(process.mix*process.digis*process.siPixelRawData*process.siPixelDigis*process.rechits*process.tracks*process.trackinghits) - -#process.p1 = cms.Path(process.mix*process.digis*process.siPixelRawData*process.siPixelDigis*process.pixeltrackerlocalreco) - -#This process to get events to be used to make DQM occupancy maps (cmsRun DQM_Pixel_digi.py after you ran testPixelDigitizer.py): -process.p1 = cms.Path(process.mix*process.simSiPixelDigis*process.siPixelRawData*process.siPixelDigis) - -# Look at cluster charge: -#process.p1 = cms.Path(process.mix*process.digis*process.siPixelRawData*process.siPixelDigis*process.pixeltrackerlocalreco) +process.g4SimHits.Generator.HepMCProductLabel = 'source' -#This process is to run the digis validation: -#process.p1 = cms.Path(process.mix*process.simSiPixelDigis*process.pixelDigisValid) +# modify digitizer parameters +#process.simSiPixelDigis.ThresholdInElectrons_BPix = 3500.0 +process.mix.digitizers.pixel.ThresholdInElectrons_BPix = 3500.0 -#This process is to run the digitizer: -#process.p1 = cms.Path(process.mix*process.simSiPixelDigis) +#This process is to run the digitizer, pixel gitizer is now clled by the mix module +process.p1 = cms.Path(process.mix) process.outpath = cms.EndPath(process.o1) -process.g4SimHits.Generator.HepMCProductLabel = 'source' + diff --git a/SimTracker/SiPixelDigitizer/test/testPxdigi_cfg.py b/SimTracker/SiPixelDigitizer/test/testPxdigi_cfg.py index 9ac8010399411..26fde277cd077 100644 --- a/SimTracker/SiPixelDigitizer/test/testPxdigi_cfg.py +++ b/SimTracker/SiPixelDigitizer/test/testPxdigi_cfg.py @@ -21,8 +21,8 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( - 'file:/afs/cern.ch/work/d/dkotlins/public/digis.root' -# 'file:/afs/cern.ch/user/d/dkotlins/scratch0/data/digis21.root' + 'file:/afs/cern.ch/work/d/dkotlins/public/MC/mu/pt100/digis/digis1.root' +# 'file:dummy_100.root' ) ) @@ -50,7 +50,9 @@ #process.GlobalTag.globaltag = 'MC_31X_V9::All' #process.GlobalTag.globaltag = 'CRAFT09_R_V4::All' # 2012 -process.GlobalTag.globaltag = 'GR_P_V40::All' +# process.GlobalTag.globaltag = 'GR_P_V40::All' +# 2013 MC +process.GlobalTag.globaltag = 'MC_70_V1::All' #process.load("Geometry.TrackerGeometryBuilder.trackerGeometry_cfi") #process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi") @@ -58,8 +60,11 @@ # include "MagneticField/Engine/data/volumeBasedMagneticField.cfi" process.analysis = cms.EDAnalyzer("PixelDigisTest", - Verbosity = cms.untracked.bool(True), - src = cms.InputTag("siPixelDigis"), + Verbosity = cms.untracked.bool(False), +# sim in V7 + src = cms.InputTag("mix"), +# old default +# src = cms.InputTag("siPixelDigis"), ) process.p = cms.Path(process.analysis) diff --git a/SimTracker/SiPixelDigitizer/test/testPxsims.py b/SimTracker/SiPixelDigitizer/test/testPxsims.py index 896f1e8c6ec1b..58a52f7079a20 100644 --- a/SimTracker/SiPixelDigitizer/test/testPxsims.py +++ b/SimTracker/SiPixelDigitizer/test/testPxsims.py @@ -23,8 +23,8 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( # '/store/user/kotlinski/mu100/simhits/simHits.root', - 'file:simHits.root' -# 'file:/scratch/kotlinski/sim/mu/pt100/mu100_10k_1.root' +# 'file:simHits.root' + 'file:/afs/cern.ch/work/d/dkotlins/public//MC/mu/pt100/simhits/simHits1.root' ) ) @@ -50,9 +50,9 @@ src = cms.string("g4SimHits"), # list = cms.string("TrackerHitsPixelBarrelLowTof"), # list = cms.string("TrackerHitsPixelBarrelHighTof"), -# list = cms.string("TrackerHitsPixelEndcapLowTof"), - list = cms.string("TrackerHitsPixelEndcapHighTof"), - Verbosity = cms.untracked.bool(True), + list = cms.string("TrackerHitsPixelEndcapLowTof"), +# list = cms.string("TrackerHitsPixelEndcapHighTof"), + Verbosity = cms.untracked.bool(False), # mode = cms.untracked.string("bpix"), mode = cms.untracked.string("fpix"), )