Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix bpix efficiency 76x #11487

Merged
merged 2 commits into from
Oct 7, 2015
Merged

fix bpix efficiency 76x #11487

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0451806
fix bpix efficiency 76x
fioriNTU Sep 25, 2015
ebcb5cf
Bpix efficiency fix plus Pixel cluster probability histos 76x
fioriNTU Sep 29, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 18 additions & 0 deletions DQM/SiPixelMonitorTrack/interface/SiPixelHitEfficiencySource.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
#include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
#include "DataFormats/SiPixelDetId/interface/PixelBarrelName.h"
#include "DataFormats/SiPixelDetId/interface/PixelBarrelNameUpgrade.h"
Expand All @@ -48,6 +49,7 @@ class SiPixelHitEfficiencySource : public DQMEDAnalyzer {
virtual void dqmBeginRun(const edm::Run& r, edm::EventSetup const& iSetup);
virtual void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
virtual void analyze(const edm::Event&, const edm::EventSetup&);
virtual void fillClusterProbability(int , int, bool, double );

private:
edm::ParameterSet pSet_;
Expand Down Expand Up @@ -86,6 +88,22 @@ class SiPixelHitEfficiencySource : public DQMEDAnalyzer {

bool isUpgrade;

//MEs for cluster probability
MonitorElement* meClusterProbabilityL1_Plus_;
MonitorElement* meClusterProbabilityL1_Minus_;

MonitorElement* meClusterProbabilityL2_Plus_;
MonitorElement* meClusterProbabilityL2_Minus_;

MonitorElement* meClusterProbabilityL3_Plus_;
MonitorElement* meClusterProbabilityL3_Minus_;

MonitorElement* meClusterProbabilityD1_Plus_;
MonitorElement* meClusterProbabilityD1_Minus_;

MonitorElement* meClusterProbabilityD2_Plus_;
MonitorElement* meClusterProbabilityD2_Minus_;

};

#endif
219 changes: 152 additions & 67 deletions DQM/SiPixelMonitorTrack/src/SiPixelHitEfficiencySource.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@

#include <iostream>
#include <map>
#include <math.h>
#include <string>
#include <vector>
#include <utility>
Expand Down Expand Up @@ -80,7 +81,7 @@ SiPixelHitEfficiencySource::SiPixelHitEfficiencySource(const edm::ParameterSet&
diskOn( pSet.getUntrackedParameter<bool>("diskOn",false) ),
isUpgrade( pSet.getUntrackedParameter<bool>("isUpgrade",false) )
//updateEfficiencies( pSet.getUntrackedParameter<bool>("updateEfficiencies",false) )
{
{
pSet_ = pSet;
debug_ = pSet_.getUntrackedParameter<bool>("debug", false);
applyEdgeCut_ = pSet_.getUntrackedParameter<bool>("applyEdgeCut");
Expand Down Expand Up @@ -111,17 +112,52 @@ SiPixelHitEfficiencySource::~SiPixelHitEfficiencySource() {
}
}

void SiPixelHitEfficiencySource::fillClusterProbability(int layer, int disk, bool plus, double probability){

//barrel
if (layer!=0){
if (layer==1){
if (plus) meClusterProbabilityL1_Plus_->Fill(probability);
else meClusterProbabilityL1_Minus_->Fill(probability);
}

else if (layer==2){
if (plus) meClusterProbabilityL2_Plus_->Fill(probability);
else meClusterProbabilityL2_Minus_->Fill(probability);
}

else if (layer==3){
if (plus) meClusterProbabilityL3_Plus_->Fill(probability);
else meClusterProbabilityL3_Minus_->Fill(probability);
}
}
//Endcap
if (disk!=0){
if (disk==1){
if (plus) meClusterProbabilityD1_Plus_->Fill(probability);
else meClusterProbabilityD1_Minus_->Fill(probability);
}
if (disk==2){
if (plus) meClusterProbabilityD2_Plus_->Fill(probability);
else meClusterProbabilityD2_Minus_->Fill(probability);
}
}
}





void SiPixelHitEfficiencySource::dqmBeginRun(const edm::Run& r, edm::EventSetup const& iSetup) {
LogInfo("PixelDQM") << "SiPixelHitEfficiencySource beginRun()" << endl;

if(firstRun){
// retrieve TrackerGeometry for pixel dets
// retrieve TrackerGeometry for pixel dets

nvalid=0;
nmissing=0;
nvalid=0;
nmissing=0;

firstRun = false;
firstRun = false;
}

edm::ESHandle<TrackerGeometry> TG;
Expand All @@ -144,7 +180,7 @@ void SiPixelHitEfficiencySource::dqmBeginRun(const edm::Run& r, edm::EventSetup
}
}
LogInfo("PixelDQM") << "SiPixelStructure size is " << theSiPixelStructure.size() << endl;

}

void SiPixelHitEfficiencySource::bookHistograms(DQMStore::IBooker & iBooker, edm::Run const & iRun, edm::EventSetup const & iSetup){
Expand Down Expand Up @@ -184,6 +220,40 @@ void SiPixelHitEfficiencySource::bookHistograms(DQMStore::IBooker & iBooker, edm
}
}

//book cluster probability histos for Barrel and Endcap
iBooker.setCurrentFolder("Pixel/Barrel");

meClusterProbabilityL1_Plus_ = iBooker.book1D("ClusterProbabilityXY_Layer1_Plus","ClusterProbabilityXY_Layer1_Plus",250,-5,0.1);
meClusterProbabilityL1_Plus_->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityL1_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer1_Minus","ClusterProbabilityXY_Layer1_Minus",250,-5,0.1);
meClusterProbabilityL1_Minus_->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityL2_Plus_ = iBooker.book1D("ClusterProbabilityXY_Layer2_Plus","ClusterProbabilityXY_Layer2_Plus",250,-5,0.1);
meClusterProbabilityL2_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityL2_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer2_Minus","ClusterProbabilityXY_Layer2_Minus",250,-5,0.1);
meClusterProbabilityL2_Minus_ ->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityL3_Plus_ = iBooker.book1D("ClusterProbabilityXY_Layer3_Plus","ClusterProbabilityXY_Layer3_Plus",250,-5,0.1);
meClusterProbabilityL3_Plus_->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityL3_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer3_Minus","ClusterProbabilityXY_Layer3_Minus",250,-5,0.1);
meClusterProbabilityL3_Minus_->setAxisTitle("Log(ClusterProbability)",1);

iBooker.setCurrentFolder("Pixel/Endcap");

meClusterProbabilityD1_Plus_ = iBooker.book1D("ClusterProbabilityXY_Disk1_Plus","ClusterProbabilityXY_Disk1_Plus",250,-5,0.1);
meClusterProbabilityD1_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityD1_Minus_ = iBooker.book1D("ClusterProbabilityXY_Disk1_Minus","ClusterProbabilityXY_Disk1_Minus",250,-5,0.1);
meClusterProbabilityD1_Minus_->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityD2_Plus_ = iBooker.book1D("ClusterProbabilityXY_Disk2_Plus","ClusterProbabilityXY_Disk2_Plus",250,-5,0.1);
meClusterProbabilityD2_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);

meClusterProbabilityD2_Minus_ = iBooker.book1D("ClusterProbabilityXY_Disk2_Minus","ClusterProbabilityXY_Disk2_Minus",250,-5,0.1);
meClusterProbabilityD2_Minus_->setAxisTitle("Log(ClusterProbability)",1);

}

Expand Down Expand Up @@ -368,64 +438,62 @@ void SiPixelHitEfficiencySource::analyze(const edm::Event& iEvent, const edm::Ev
continue;
}

TrajectoryStateOnSurface predTrajState=expTrajMeasurements[iexp].updatedState();
float x=predTrajState.globalPosition().x();
float y=predTrajState.globalPosition().y();
float z=predTrajState.globalPosition().z();
float dxyz=sqrt((chkx-x)*(chkx-x)+(chky-y)*(chky-y)+(chkz-z)*(chkz-z));

if (dxyz<=glmatch) {
glmatch=dxyz;
imatch=iexp;
imatches.push_back(int(imatch));
}

} // found the propagated traj best matching the hit in data

float lxmatch = 9999.0;
float lymatch = 9999.0;
if(!expTrajMeasurements.empty()){
if (glmatch<9999.) { // if there is any propagated trajectory for this hit
const DetId & matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();
TrajectoryStateOnSurface predTrajState=expTrajMeasurements[iexp].updatedState();
float x=predTrajState.globalPosition().x();
float y=predTrajState.globalPosition().y();
float z=predTrajState.globalPosition().z();
float dxyz=sqrt((chkx-x)*(chkx-x)+(chky-y)*(chky-y)+(chkz-z)*(chkz-z));

if (dxyz<=glmatch) {
glmatch=dxyz;
imatch=iexp;
imatches.push_back(int(imatch));
}

} // found the propagated traj best matching the hit in data

float lxmatch = 9999.0;
float lymatch = 9999.0;
if(!expTrajMeasurements.empty()){
if (glmatch<9999.) { // if there is any propagated trajectory for this hit
const DetId & matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();

int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
int dladder = abs(matchhit_ladder-hit_ladder);
if (dladder > 10) dladder = 20 - dladder;
LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
lxmatch=fabs(lp.x() - chklp.x());
lymatch=fabs(lp.y() - chklp.y());
}
if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {
int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
int dladder = abs(matchhit_ladder-hit_ladder);
if (dladder > 10) dladder = 20 - dladder;
LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
lxmatch=fabs(lp.x() - chklp.x());
lymatch=fabs(lp.y() - chklp.y());
}
if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {

if (testhit->getType()!=TrackingRecHit::missing || keepOriginalMissingHit_) {
expTrajMeasurements.erase(expTrajMeasurements.begin()+imatch);
if (testhit->getType()!=TrackingRecHit::missing || keepOriginalMissingHit_) {
expTrajMeasurements.erase(expTrajMeasurements.begin()+imatch);
}

}

} //expected trajectory measurment not empty
}

}

} //expected trajectory measurment not empty
}
}//loop on trajectory measurments tmeasColl

//if an extrapolated hit was found but not matched to an exisitng L1 hit then push the hit back into the collection
//now keep the first one that is left
if(!expTrajMeasurements.empty()){
for (size_t f=0; f<expTrajMeasurements.size(); f++) {
TrajectoryMeasurement AddHit=expTrajMeasurements[f];
if (AddHit.recHit()->getType()==TrackingRecHit::missing){
tmeasColl.push_back(AddHit);
isBpixtrack = true;
//now keep the first one that is left
if(!expTrajMeasurements.empty()){
for (size_t f=0; f<expTrajMeasurements.size(); f++) {
TrajectoryMeasurement AddHit=expTrajMeasurements[f];
if (AddHit.recHit()->getType()==TrackingRecHit::missing){
tmeasColl.push_back(AddHit);
isBpixtrack = true;

}
}

}



}
}

if(isBpixtrack || isFpixtrack){
if(trackref->pt()<0.6 ||
nStripHits<11 ||
nStripHits<8 ||
fabs(trackref->dxy(bestVtx->position()))>0.05 ||
fabs(trackref->dz(bestVtx->position()))>0.5) continue;

Expand All @@ -443,7 +511,7 @@ float y=predTrajState.globalPosition().y();
continue;
else {

// //residual
// //residual
const DetId & hit_detId = hit->geographicalId();
//uint IntRawDetID = (hit_detId.rawId());
uint IntSubDetID = (hit_detId.subdetId());
Expand All @@ -456,43 +524,59 @@ float y=predTrajState.globalPosition().y();
continue;

int disk=0; int layer=0; int panel=0; int module=0; bool isHalfModule=false;

const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit->hit());

if(IntSubDetID==PixelSubdetector::PixelBarrel){ // it's a BPIX hit
layer = PixelBarrelName(hit_detId,pTT,isUpgrade).layerName();
isHalfModule = PixelBarrelName(hit_detId,pTT,isUpgrade).isHalfModule();

if (hit->isValid()){ //fill the cluster probability in barrel
bool plus=true;
if ((PixelBarrelName(hit_detId,pTT,isUpgrade).shell()== PixelBarrelName::Shell::mO) || (PixelBarrelName(hit_detId,pTT,isUpgrade).shell()==PixelBarrelName::Shell::mI)) plus=false;
double clusterProbability= pixhit->clusterProbability(0);
if (clusterProbability!=0) fillClusterProbability(layer,0,plus,log10(clusterProbability));
}

}else if(IntSubDetID==PixelSubdetector::PixelEndcap){ // it's an FPIX hit
disk = PixelEndcapName(hit_detId,pTT,isUpgrade).diskName();
panel = PixelEndcapName(hit_detId,pTT,isUpgrade).pannelName();
module = PixelEndcapName(hit_detId,pTT,isUpgrade).plaquetteName();

if (hit->isValid()){
bool plus=true;
if ((PixelEndcapName(hit_detId,pTT,isUpgrade).halfCylinder()== PixelEndcapName::HalfCylinder::mO) || (PixelEndcapName(hit_detId,pTT,isUpgrade).halfCylinder()==PixelEndcapName::HalfCylinder::mI)) plus=false;
double clusterProbability= pixhit->clusterProbability(0);
if (clusterProbability!=0) fillClusterProbability(0,disk,plus,log10(clusterProbability));
}
}

if(layer==1){
if(fabs(trackref->dxy(bestVtx->position()))>0.01 ||
fabs(trackref->dz(bestVtx->position()))>0.1) continue;
if(!(L2hits>0&&L3hits>0&&L4hits>0) && !(L2hits>0&&D1hits>0&&D2hits) && !(D1hits>0&&D2hits>0&&D3hits>0)) continue;
if(!(L2hits>0&&L3hits>0) && !(L2hits>0&&D1hits>0) && !(D1hits>0&&D2hits>0)) continue;
}else if(layer==2){
if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
fabs(trackref->dz(bestVtx->position()))>0.1) continue;
if(!(L1hits>0&&L3hits>0&&L4hits>0) && !(L1hits>0&&L3hits>0&&D1hits>0) && !(L1hits>0&&D1hits>0&&D2hits>0)) continue;
if(!(L1hits>0&&L3hits>0) && !(L1hits>0&&D1hits>0)) continue;
}else if(layer==3){
if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
fabs(trackref->dz(bestVtx->position()))>0.1) continue;
if(!(L1hits>0&&L2hits>0&&L4hits>0) && !(L1hits>0&&L2hits>0&&D1hits>0)) continue;
if(!(L1hits>0&&L2hits>0)) continue;
}else if(layer==4){
if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
fabs(trackref->dz(bestVtx->position()))>0.1) continue;
if(!(L1hits>0&&L2hits>0&&L3hits>0)) continue;
}else if(disk==1){
if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
fabs(trackref->dz(bestVtx->position()))>0.5) continue;
if(!(L1hits>0&&L2hits>0&&D2hits>0) && !(L1hits>0&&D2hits>0&&D3hits>0) && !(L2hits>0&&D2hits>0&&D3hits>0)) continue;
if(!(L1hits>0&&D2hits>0) && !(L2hits>0&&D2hits>0)) continue;
}else if(disk==2){
if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
fabs(trackref->dz(bestVtx->position()))>0.5) continue;
if(!(L1hits>0&&L2hits>0&&D1hits>0) && !(L1hits>0&&D1hits>0&&D3hits>0) && !(L2hits>0&&D1hits>0&&D3hits>0)) continue;
if(!(L1hits>0&&D1hits>0)) continue;
}else if(disk==3){
if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
fabs(trackref->dz(bestVtx->position()))>0.5) continue;
if(!(L1hits>0&&D1hits>0&&D2hits>0) && !(L2hits>0&&D1hits>0&&D2hits>0)) continue;
}

//check whether hit is valid or missing using track algo flag
Expand All @@ -514,6 +598,7 @@ float y=predTrajState.globalPosition().y();
if(fabs(lx)>0.55 || fabs(ly)>3.0) continue;

bool passedFiducial=true;

// Module fiducials:
if(IntSubDetID==PixelSubdetector::PixelBarrel && fabs(ly)>=3.1) passedFiducial=false;
if(IntSubDetID==PixelSubdetector::PixelEndcap &&
Expand Down Expand Up @@ -623,18 +708,18 @@ float y=predTrajState.globalPosition().y();
float d_cl[2]; d_cl[0]=d_cl[1]=-9999.;
if(dx_cl[0]!=-9999. && dy_cl[0]!=-9999.) d_cl[0]=sqrt(dx_cl[0]*dx_cl[0]+dy_cl[0]*dy_cl[0]);
if(dx_cl[1]!=-9999. && dy_cl[1]!=-9999.) d_cl[1]=sqrt(dx_cl[1]*dx_cl[1]+dy_cl[1]*dy_cl[1]);
if(isHitMissing && (d_cl[0]<0.05 || d_cl[1]<0.05)){ isHitMissing=0; isHitValid=1; }

if(isHitMissing && (d_cl[0]<0.05 || d_cl[1]<0.05)){ isHitMissing=0; isHitValid=1; }

if(debug_){
std::cout << "Ready to add hit in histogram:\n";
//std::cout << "detid: "<<hit_detId<<std::endl;
std::cout << "isHitValid: "<<isHitValid<<std::endl;
std::cout << "isHitMissing: "<<isHitMissing<<std::endl;
//std::cout << "passedEdgeCut: "<<passedFiducial<<std::endl;
}


}

if (nStripHits<11) continue; //Efficiency plots are filled with hits on tracks that have at least 11 Strip hits

if(pxd!=theSiPixelStructure.end() && isHitValid && passedFiducial)
++nvalid;
if(pxd!=theSiPixelStructure.end() && isHitMissing && passedFiducial)
Expand Down